ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc
(Generate patch)

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.2 by jengbou, Thu Sep 3 19:29:38 2009 UTC vs.
Revision 1.3 by jengbou, Tue Sep 8 07:15:38 2009 UTC

# Line 59 | Line 59 | Implementation:
59   using namespace std;
60  
61   PVStudy::PVStudy(const edm::ParameterSet& iConfig):
62 <  nTrkMin_(10),nTrkMax_(150)
62 >  nTrkMin_(10),nTrkMax_(100)
63   {
64    simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
65    tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
# Line 67 | Line 67 | PVStudy::PVStudy(const edm::ParameterSet
67    verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
68    realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
69    analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
70 +  outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
71  
72    //now do what ever initialization is needed
73 <  pvanainfo.clear();
73 >  pvinfo.clear();
74  
75 +  file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
76 +  ftree_ = new TTree("mytree","mytree");
77 +  ftree_->AutoSave();
78 +  ftree_->Branch("residual",&fres,"fres[3]/D");
79 +  ftree_->Branch("error",&ferror,"ferror[3]/D");
80 +  ftree_->Branch("nTrk",&fntrk,"fntrk/I");
81 +  
82    edm::Service<TFileService> fs;
83   //   TFileDirectory subDir = fs->mkdir( "Summary" );
84   //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
# Line 92 | Line 100 | PVStudy::PVStudy(const edm::ParameterSet
100      if (analyze_) {
101        h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
102        h_resx_nTrk->SetMarkerStyle(21);
103 +      h_resx_nTrk->SetMarkerColor(4);
104 +      h_resx_nTrk->GetXaxis()->SetTitle("Num of tracks");
105        h_resx_nTrk->GetYaxis()->SetTitle("#mum");
106        h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
107        h_resy_nTrk->SetMarkerStyle(21);
108 +      h_resy_nTrk->SetMarkerColor(4);
109 +      h_resy_nTrk->GetXaxis()->SetTitle("Num of tracks");
110        h_resy_nTrk->GetYaxis()->SetTitle("#mum");
111        h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
112        h_resz_nTrk->SetMarkerStyle(21);
113 +      h_resz_nTrk->SetMarkerColor(4);
114 +      h_resz_nTrk->GetXaxis()->SetTitle("Num of tracks");
115        h_resz_nTrk->GetYaxis()->SetTitle("#mum");
116        h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
117        h_pullx_nTrk->SetMarkerStyle(21);
118 +      h_pullx_nTrk->SetMarkerColor(4);
119        h_pullx_nTrk->SetBit(TH1::kCanRebin);
120 <      h_pullx_nTrk->GetYaxis()->SetTitle("cm");
120 >      h_pullx_nTrk->GetXaxis()->SetTitle("Num of tracks");
121        h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
122        h_pully_nTrk->SetMarkerStyle(21);
123 +      h_pully_nTrk->SetMarkerColor(4);
124        h_pully_nTrk->SetBit(TH1::kCanRebin);
125 <      h_pully_nTrk->GetYaxis()->SetTitle("cm");
125 >      h_pully_nTrk->GetXaxis()->SetTitle("Num of tracks");
126        h_pullz_nTrk   = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
127        h_pullz_nTrk->SetMarkerStyle(21);
128 +      h_pullz_nTrk->SetMarkerColor(4);
129        h_pullz_nTrk->SetBit(TH1::kCanRebin);
130 <      h_pullz_nTrk->GetYaxis()->SetTitle("cm");
130 >      h_pullz_nTrk->GetXaxis()->SetTitle("Num of tracks");
131      }
132      for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
133        stringstream ss;
134        ss << ntrk;
135        string suffix = ss.str();
136 <      h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",400,-0.02,0.02);
136 >      h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
137        h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
138 <      h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",400,-0.02,0.02);
138 >      h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
139        h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
140 <      h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",400,-0.02,0.02);
140 >      h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
141        h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
142 <      h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",200,-10.,10.);
143 <      h["pullx_"+suffix]->GetXaxis()->SetTitle("cm");
144 <      h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",200,-10.,10.);
128 <      h["pully_"+suffix]->GetXaxis()->SetTitle("cm");
129 <      h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",200,-10.,10.);
130 <      h["pullz_"+suffix]->GetXaxis()->SetTitle("cm");
142 >      h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
143 >      h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
144 >      h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
145        suffix.clear();
146      }
147    }
# Line 305 | Line 319 | void
319   PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
320   {
321    using namespace edm;
308 –
322    using reco::TrackCollection;
323 <  
323 >
324 >  ftree_->SetBranchAddress("residual",&fres);
325 >  ftree_->SetBranchAddress("error",&ferror);
326 >  ftree_->SetBranchAddress("nTrk",&fntrk);
327 >
328    Handle<TrackCollection> tracks;
329    iEvent.getByLabel(tracksLabel_,tracks);
330  
# Line 379 | Line 396 | PVStudy::analyze(const edm::Event& iEven
396            if(verbose_) {
397              cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
398            }
399 +          fres[0] = vsim->recVtx->x()-vsim->x;
400 +          fres[1] = vsim->recVtx->y()-vsim->y;
401 +          fres[2] = vsim->recVtx->z()-vsim->z;
402 +          ferror[0] = vsim->recVtx->xError();
403 +          ferror[1] = vsim->recVtx->yError();
404 +          ferror[2] = vsim->recVtx->zError();
405 +          fntrk = nrectrk;
406 +          ftree_->Fill();
407 +
408            h_deltax->Fill( vsim->recVtx->x()-vsim->x );
409            h_deltay->Fill( vsim->recVtx->y()-vsim->y );
410            h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
411            h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
412            h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
413            h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
414 <          pvanainfo.push_back(PVStudy::PVAnaInfo( vsim->recVtx->x()-vsim->x,
415 <                                                  vsim->recVtx->y()-vsim->y,
416 <                                                  vsim->recVtx->z()-vsim->z,
417 <                                                  (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
418 <                                                  (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
419 <                                                  (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError(),
420 <                                                  nrectrk )
421 <                              );
414 >          pvinfo.push_back(PVStudy::PVInfo(res(vsim->recVtx->x()-vsim->x,
415 >                                               vsim->recVtx->y()-vsim->y,
416 >                                               vsim->recVtx->z()-vsim->z),
417 >                                           error(vsim->recVtx->xError(),
418 >                                                 vsim->recVtx->yError(),
419 >                                                 vsim->recVtx->zError()),
420 >                                           nrectrk));
421 >          // Fill histo according to its track multiplicity
422 >          fillHisto(res(vsim->recVtx->x()-vsim->x,
423 >                        vsim->recVtx->y()-vsim->y,
424 >                        vsim->recVtx->z()-vsim->z),
425 >                    pull((vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
426 >                         (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
427 >                         (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError()),
428 >                    nrectrk);
429          }
430          else{  // no rec vertex found for this simvertex
431            if(verbose_) {
# Line 456 | Line 489 | PVStudy::endJob() {
489        if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
490        PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
491        if (analyze_) {
492 <        if ( ipv.resx > 0 ) h_resx_nTrk->Fill(ntrk,ipv.resx/sqt2, 1.);
493 <        if ( ipv.resy > 0 ) h_resy_nTrk->Fill(ntrk,ipv.resy/sqt2, 1.);
494 <        if ( ipv.resz > 0 ) h_resz_nTrk->Fill(ntrk,ipv.resz/sqt2, 1.);
495 <        if ( ipv.pullx > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pullx, 1.);
496 <        if ( ipv.pully > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pully, 1.);
497 <        if ( ipv.pullz > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pullz, 1.);
492 > //      if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
493 > //      if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
494 > //      if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
495 >        if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x(), 1.);
496 >        if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y(), 1.);
497 >        if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z(), 1.);
498 >        if ( ipv.pull_.x() > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pull_.x(), 1.);
499 >        if ( ipv.pull_.y() > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pull_.y(), 1.);
500 >        if ( ipv.pull_.z() > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pull_.z(), 1.);
501        }
502      }
503    }
504 +  file_->cd();
505 +  ftree_->Write();
506 +  file_->Close();
507   }
508  
509   bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
# Line 544 | Line 583 | void PVStudy::printSimTrks(const Handle<
583    }
584   }
585  
586 + void PVStudy::fillHisto(res r, pull p, int ntk){
587 +  stringstream ss;
588 +  ss << ntk;
589 +  string suffix = ss.str();
590 +  if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
591 +  // Fill histograms of res and pull of ntk
592 +  h["resx_"+suffix]->Fill(r.x());
593 +  h["resy_"+suffix]->Fill(r.y());
594 +  h["resz_"+suffix]->Fill(r.z());
595 +  h["pullx_"+suffix]->Fill(p.x());
596 +  h["pully_"+suffix]->Fill(p.y());
597 +  h["pullz_"+suffix]->Fill(p.z());
598 + }
599 +
600   PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
601    map<int,double> data;
602    data.clear();
550 –  double xmin=0.,xmax=0.;
603    double fpar[2];
604    double cm2um = 10000;
605    stringstream ss;
606    ss << ntk;
607    string suffix = ss.str();
608  
609 <  for ( int i = 0; i < 6; ++i) {
610 <    switch(i){
611 <    case 0:
612 <    case 1:
613 <    case 2:
614 <      if (verbose_) cout << "Analyzing Resolution:" << endl;
615 <      xmin = -0.01; xmax = 0.01;
616 <      break;
617 <    default:
618 <      if (verbose_) cout << "Analyzing Pull:" << endl;
619 <      xmin = -10.; xmax = 10;
620 <      break;
621 <    }
622 <    TH1D *hdist = new TH1D("hdist","distribution to be fitted", 1000,xmin,xmax);
623 <
624 <    for (std::vector<PVStudy::PVAnaInfo>::const_iterator ipvinfo = pvanainfo.begin();
625 <         ipvinfo != pvanainfo.end(); ++ipvinfo){
626 <      if ( verbose_ && i==0 ) {
627 <        cout << "Num of tracks of the PV = " << ipvinfo->nTrk;
628 <        cout << "; delta x = " << ipvinfo->resx << endl;
629 <      }
630 <      if (ipvinfo->nTrk == ntk) {
631 <        switch(i) {
632 <        case 0: if (verbose_) cout << "Matched resx" << endl;
633 <          hdist->Fill(ipvinfo->resx);
634 <          h["resx_"+suffix]->Fill(ipvinfo->resx);
635 <          break;
636 <        case 1: if (verbose_) cout << "Matched resy" << endl;
585 <          hdist->Fill(ipvinfo->resy);
586 <          h["resy_"+suffix]->Fill(ipvinfo->resy);
587 <          break;
588 <        case 2: if (verbose_) cout << "Matched resz" << endl;
589 <          hdist->Fill(ipvinfo->resz);
590 <          h["resz_"+suffix]->Fill(ipvinfo->resz);
591 <          break;
592 <        case 3: if (verbose_) cout << "Matched pullx" << endl;
593 <          hdist->Fill(ipvinfo->pullx);
594 <          h["pullx_"+suffix]->Fill(ipvinfo->pullx);
595 <          break;
596 <        case 4: if (verbose_) cout << "Matched pully" << endl;
597 <          hdist->Fill(ipvinfo->pully);
598 <          h["pully_"+suffix]->Fill(ipvinfo->pully);
599 <          break;
600 <        case 5: if (verbose_) cout << "Matched pullz" << endl;
601 <          hdist->Fill(ipvinfo->pullz);
602 <          h["pullz_"+suffix]->Fill(ipvinfo->pullz);
603 <          break;
604 <        }
605 <      }
606 <      else {
607 <        //      switch(i) {
608 <        //      case 0: if (verbose_) cout << "bad resx" << endl; break;
609 <        //      case 1: if (verbose_) cout << "bad resy" << endl; break;
610 <        //      case 2: if (verbose_) cout << "bad resz" << endl; break;
611 <        //      case 3: if (verbose_) cout << "bad pullx" << endl; break;
612 <        //      case 4: if (verbose_) cout << "bad pully" << endl; break;
613 <        //      case 5: if (verbose_) cout << "bad pullz" << endl; break;
614 <        //      }
609 >  if(analyze_) {
610 >    for ( int i = 0; i < 6; ++i) {
611 >      cout << "Here" << endl;
612 >      switch (i) {
613 >      case 0:
614 >        fit(h["resx_"+suffix], fpar);
615 >        data[i] = fpar[0]*cm2um;
616 >        break;
617 >      case 1:
618 >        fit(h["resy_"+suffix], fpar);
619 >        data[i] = fpar[0]*cm2um;
620 >        break;
621 >      case 2:
622 >        fit(h["resz_"+suffix], fpar);
623 >        data[i] = fpar[0]*cm2um;
624 >        break;
625 >      case 3:
626 >        fit(h["pullx_"+suffix], fpar);
627 >        data[i] = fpar[0];
628 >        break;
629 >      case 4:
630 >        fit(h["pully_"+suffix], fpar);
631 >        data[i] = fpar[0];
632 >        break;
633 >      case 5:
634 >        fit(h["pullz_"+suffix], fpar);
635 >        data[i] = fpar[0];
636 >        break;
637        }
638 +      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
639      }
617 –    if(analyze_) {
618 –      TAxis *axis0 = hdist->GetXaxis();
619 –      int nbin = axis0->GetLast();
620 –      double nOF = hdist->GetBinContent(nbin+1);
621 –      double nUF = hdist->GetBinContent(0);
622 –      if ( verbose_ && i==0 ){
623 –        cout << "Last bin = " << nbin;
624 –        cout << "; hdist: Overflow Entries = " << nOF;
625 –        cout << "; Underflow Entries = " << nUF;
626 –        cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
627 –      }
628 –      if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
629 –        if(i<3){
630 –          hdist->Fit("gaus","Q0");
631 –          TF1 *fgaus = hdist->GetFunction("gaus");
632 –          if (verbose_) cout << "got function" << endl;
633 –          fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
634 –          fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
635 –          delete fgaus;
636 –        }
637 –        else {
638 –          hdist->Fit("gausn","Q0");
639 –          TF1 *fgaus = hdist->GetFunction("gausn");
640 –          if (verbose_) cout << "got function" << endl;
641 –          fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
642 –          fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
643 –          delete fgaus;
644 –        }
645 –        if (verbose_) cout << "fpar[2] = (" << fpar[0] << "," << fpar[1] << ")" << endl;
646 –        switch (i) {
647 –        case 0:
648 –          h["resx_"+suffix]->Fit("gaus","Q");
649 –          data[i] = fpar[1]*cm2um;
650 –          break;
651 –        case 1:
652 –          h["resy_"+suffix]->Fit("gaus","Q");
653 –          data[i] = fpar[1]*cm2um;
654 –          break;
655 –        case 2:
656 –          h["resz_"+suffix]->Fit("gaus","Q");
657 –          data[i] = fpar[1]*cm2um;
658 –          break;
659 –        case 3:
660 –          h["pullx_"+suffix]->Fit("gausn","Q");
661 –          data[i] = fpar[1];
662 –          break;
663 –        case 4:
664 –          h["pully_"+suffix]->Fit("gausn","Q");
665 –          data[i] = fpar[1];
666 –          break;
667 –        case 5:
668 –          h["pullz_"+suffix]->Fit("gausn","Q");
669 –          data[i] = fpar[1];
670 –          break;
671 –        }
672 –        if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" cm") << endl;
673 –      }
674 –      else {
675 –        if ( verbose_) cout << "Empty Histo! No Fit!!!" << endl;
676 –        data[i] = -999;
677 –      }
678 –    }
679 –    else
680 –      data[i] = -999;
681 –    delete hdist;
640    }
641 +  else
642 +    for ( int i = 0; i < 6; ++i) {
643 +      data[i] = -999;
644 +    }
645  
646 <  return PVStudy::PVAnaInfo(data[0], data[1], data[2],
647 <                            data[3], data[4], data[5],
646 >  return PVStudy::PVAnaInfo(res(data[0], data[1], data[2]),
647 >                            error(0.,0.,0.),
648 >                            pull(data[3], data[4], data[5]),
649 >                            error(0.,0.,0.),
650                              ntk);
651   }
652  
653 + void PVStudy::fit(TH1 *&hdist, double data[]){
654 +  TAxis *axis0 = hdist->GetXaxis();
655 +  int nbin = axis0->GetLast();
656 +  double nOF = hdist->GetBinContent(nbin+1);
657 +  double nUF = hdist->GetBinContent(0);
658 +  if ( verbose_ ){
659 +    cout << "Last bin = " << nbin;
660 +    cout << "; hdist: Overflow Entries = " << nOF;
661 +    cout << "; Underflow Entries = " << nUF;
662 +    cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
663 +  }
664 +  if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
665 +    hdist->Fit("gaus","Q");
666 +    TF1 *fgaus = hdist->GetFunction("gaus");
667 +    if (verbose_) cout << "got function" << endl;
668 +    data[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
669 +  }
670 + }
671 +
672   //define this as a plug-in
673   DEFINE_FWK_MODULE(PVStudy);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines