TrackProducer/test/analyze_HZZ.C

0001 {
0002   // book histograms
0003   TFile histofile("HZZ.root","RECREATE");
0004   TH1F* h_chisqtr = new TH1F("chisqtr","Track chisq",100,0.,40.);
0005   TH1F* h_pttr = new TH1F("pttr","Track pT (GeV)",100,0.0,200.0);
0006   TH1F* h_pxtr = new TH1F("pxtr","Track px (GeV)",100,0.0,200.0);
0007   TH1F* h_pytr = new TH1F("pytr","Track py (GeV)",100,0.0,200.0);
0008   TH1F* h_pztr = new TH1F("pztr","Track pz (GeV)",100,0.0,200.0);
0009   TH1F* h_etatr = new TH1F("etatr","Track#eta",60,-3.0,3.0);
0010   TH1F* h_chargetr = new TH1F("chargetr","Track charge",6,-3.0,3.0);
0011   TH1F* h_mmumu = new TH1F("mmumu","mu+mu- mass (GeV)",100,0.0,200.0);
0012   TH1F* h_noz = new TH1F("noz","num of reco Z",10,0.0,10.0);
0013   TH1F* h_mzz = new TH1F("mzz","Z0Z0 mass (GeV)",100,100.0,300.0);
0014   TH1F* h_m4mu = new TH1F("m4mu","2mu+2mu- mass (GeV)",100,100.0,300.0);
0015
0016   TTree *tree = (TTree*)file.Get("Events");
0017
0018   std::vector<reco::Track> trackCollection;
0019
0020   TBranch *branch = tree->GetBranch("recoTracks_generalTracks__RECO.obj");
0021   branch->SetAddress(&trackCollection);
0022
0023   for ( unsigned int index = 0; index < tree->GetEntries(); ++index ) {
0024     std::cout << "Event: " << index << std::endl;
0025     branch->GetEntry(index);
0026     double px4=0.0, py4=0.0, pz4=0.0, e4=0.0;
0027     int q4=0, n4=0;
0028     for ( unsigned int bindex = 0; bindex < trackCollection.size(); ++bindex ) {
0029       reco::Track* track = (reco::Track*)trackCollection[bindex];
0030       h_chisqtr->Fill(track->normalizedChi2());
0031       double pT = sqrt(track->px()*track->px()+track->py()*track->py());
0032       h_pttr->Fill(pT);
0033       h_pxtr->Fill(track->px());
0034       h_pytr->Fill(track->py());
0035       h_pztr->Fill(track->pz());
0036       h_etatr->Fill(track->eta());
0037       h_chargetr->Fill(track->charge());
0038       n4++;
0039       if (track->charge()>0.0){q4++;}else{q4--;}
0040       px4+=track->px();
0041       py4+=track->py();
0042       pz4+=track->pz();
0043       e4+=sqrt(track->px()*track->px()+track->py()*track->py()
0044               +track->pz()*track->pz()+0.011163691);
0045     }
0046     double ptot = sqrt( px4*px4 + py4*py4 + pz4*pz4 );
0047     double mz = sqrt((e4+ptot)*(e4-ptot));
0048     if ((4==n4)&&(0==q4))h_m4mu->Fill(mz);
0049     std::vector<double> Zpx;
0050     std::vector<double> Zpy;
0051     std::vector<double> Zpz;
0052     std::vector<double> Ze;
0053     std::vector<int> Zpart1;
0054     std::vector<int> Zpart2;
0055     if (trackCollection.size() >1){
0056       for ( unsigned int bindex = 0; bindex < trackCollection.size()-1; ++bindex ) {
0057         reco::Track* track1 = (reco::Track*)trackCollection[bindex];
0058         for ( unsigned int cindex = bindex+1; cindex < trackCollection.size(); ++cindex ) {
0059           reco::Track* track2 = (reco::Track*)trackCollection[cindex];
0060           if (track1->charge()*track2->charge() < 0.0){
0061             double e1 = sqrt((track1->px()*track1->px())
0062                             +(track1->py()*track1->py())
0063                             +(track1->pz()*track1->pz())+0.011163691);
0064             double e2 = sqrt((track2->px()*track2->px())
0065                             +(track2->py()*track2->py())
0066                             +(track2->pz()*track2->pz())+0.011163691);
0067             double etot = e1+e2;
0068             double pxtot = track1->px()+track2->px();
0069             double pytot = track1->py()+track2->py();
0070             double pztot = track1->pz()+track2->pz();
0071             double ptot = sqrt( pxtot*pxtot + pytot*pytot + pztot*pztot );
0072             double mz = sqrt((etot+ptot)*(etot-ptot));
0073             h_mmumu->Fill(mz);
0074             if ((mz>80.0)&&(mz<100.0)){
0075               Zpx.push_back(pxtot); Zpy.push_back(pytot); Zpz.push_back(pztot);
0076               Ze.push_back(etot); Zpart1.push_back(bindex); Zpart2.push_back(cindex);
0077             }
0078           }//tracks opposite charge
0079         }//end track2
0080       }//end track1
0081     }//end got >1 trk
0082     h_noz->Fill( Zpx.size() );
0083     if (Zpx.size() >1){
0084       for ( unsigned int bindex = 0; bindex < Zpx.size()-1; ++bindex ) {
0085         for ( unsigned int cindex = bindex+1; cindex < Zpx.size(); ++cindex ) {
0086           if ((Zpart1[bindex]!=Zpart1[cindex])&&(Zpart1[bindex]!=Zpart2[cindex])
0087             &&(Zpart2[bindex]!=Zpart1[cindex])&&(Zpart2[bindex]!=Zpart2[cindex])){
0088               double etot = Ze[bindex]+Ze[cindex];
0089               double pxtot = Zpx[bindex]+Zpx[cindex];
0090               double pytot = Zpy[bindex]+Zpy[cindex];
0091               double pztot = Zpz[bindex]+Zpz[cindex];
0092               double ptot = sqrt( pxtot*pxtot + pytot*pytot + pztot*pztot );
0093               double mh = sqrt((etot+ptot)*(etot-ptot));
0094               h_mzz->Fill(mh);
0095           }
0096         }
0097       }
0098     }
0099   }//end evt loop
0100
0101
0102   // save histograms
0103   histofile.Write();
0104   histofile.Close();
0105 }