source: huonglan/ANA2011/domatching.C

Last change on this file was 6, checked in by huonglan, 13 years ago

first full version of 2011 analysis
File size: 15.8 KB
Line 
1void domatching(bool topprime)
2{
3  int pdgTop = 0;
4
5  if (topprime)
6    pdgTop = 8;
7  else 
8    pdgTop = 6;
9
10  TLorentzVector vtop, vtopbar;
11  TLorentzVector W1, W2;
12  for (int i = 0; i < mc_n; i++){
13    if (mc_status->at(i)==3 && mc_pdgId->at(i)==pdgTop)
14      vtop.SetPtEtaPhiM(mc_pt->at(i),mc_eta->at(i),mc_phi->at(i),mc_m->at(i));
15    if (mc_status->at(i)==3 && mc_pdgId->at(i)==-pdgTop)
16      vtopbar.SetPtEtaPhiM(mc_pt->at(i),mc_eta->at(i),mc_phi->at(i),mc_m->at(i));
17    if (mc_status->at(i)==3 && mc_pdgId->at(i)==24)
18      W1.SetPtEtaPhiM(mc_pt->at(i),mc_eta->at(i),mc_phi->at(i),mc_m->at(i));
19    if (mc_status->at(i)==3 && mc_pdgId->at(i)==-24)
20      W2.SetPtEtaPhiM(mc_pt->at(i),mc_eta->at(i),mc_phi->at(i),mc_m->at(i));
21  }
22 
23  //TLorentzVector vttpairMC = vtop + vtopbar;
24  //m_pTttpairMC->Fill(fabs(vttpairMC.Pt())/1000.);
25 
26  //############## USE MC DATA TO FIND JETS ##############
27  TLorentzVector vec0(0,0,0,0);
28  TLorentzVector vje, vneu;
29  TLorentzVector vj1, vj2, vqj1, vqj2;
30  vector<TLorentzVector> matchedjets;
31
32  TLorentzVector vq1,vq2,ve1,vn1;
33  TLorentzVector vqi,vqj,ve2,vn2;
34  TLorentzVector vq11,vq22,vele;
35  TLorentzVector vqt1, vqt2, vqt11, vqt22;
36  double pTqt1 = 0, pTqt2 = 0;
37  int pdgqt1 = 0, pdgqt2 = 0;
38  double dRmin = 1e9;
39
40  int barele = -1;
41
42  double match[4] = {-1,-1,-1,-1};
43 
44  //find the ELECTRON with stat=1
45  for (int i = 0; i < mc_n; i++){
46    int sta = mc_status->at(i);
47    int pdg = mc_pdgId->at(i);
48    if (sta==3 && abs(pdg)==11){
49      double pTe3 = mc_pt->at(i);
50      int pdge3 = pdg;
51      ve1.SetPtEtaPhiM(pTe3,mc_eta->at(i),mc_phi->at(i),mc_m->at(i));       
52      for (int j = 0; j < mc_n; j++){
53        int sta1 = mc_status->at(j);
54        int pdg1 = mc_pdgId->at(j);
55        if (sta1==1 && pdg1==pdge3){
56          double pTe1 = mc_pt->at(j);
57          ve2.SetPtEtaPhiM(pTe1,mc_eta->at(j),mc_phi->at(j),mc_m->at(j));
58          double dRee = ve1.DeltaR(ve2);
59          if (dRee < dRmin){
60            dRmin = dRee;
61            vele  = ve2;
62            barele = mc_barcode->at(j);
63          }         
64        } 
65      }
66    }
67  }
68   
69  if (!(vele==vec0)) {
70    //re-put return kTRUE
71    //if (vele==vec0) return kTRUE;
72
73    m_loose1->Fill(2);//Normally here we have an efficiency of 100%
74
75    //find the NEUTRINO with stat=1
76    dRmin = 1e9;
77    for (int i = 0; i < mc_n; i++){
78      if (mc_status->at(i)==3 && abs(mc_pdgId->at(i))==12){
79        double pTn3 = mc_pt->at(i);
80        int pdgn3 = mc_pdgId->at(i);
81        vn1.SetPtEtaPhiM(pTn3,mc_eta->at(i),mc_phi->at(i),mc_m->at(i));       
82        for (int j = 0; j < mc_n; j++){
83          if (mc_status->at(j)==1 && mc_pdgId->at(j)==pdgn3){
84            double pTn1 = mc_pt->at(j);
85            vn2.SetPtEtaPhiM(pTn1,mc_eta->at(j),mc_phi->at(j),mc_m->at(j));
86            double dRnn = vn1.DeltaR(vn2);
87            if (dRnn < dRmin){
88              dRmin = dRnn;
89              vneu  = vn2;
90            }         
91          } 
92        }
93      }
94    }
95
96    if (!(vneu==vec0)) {
97      //re-put return kTRUE
98      //if (vneu==vec0) return kTRUE;
99     
100      m_loose1->Fill(3);//Normally here we have an efficiency of 100%
101       
102      //Find the quarks from W and from top with status=2
103      //The first particles in the event have status = 3
104      //Check if the two particle 9 and 10 are quarks or not
105      if (abs(mc_pdgId->at(9))<=5 && abs(mc_pdgId->at(10))<=5){
106        vq1.SetPtEtaPhiM(mc_pt->at(9),mc_eta->at(9),mc_phi->at(9),mc_m->at(9));
107        vq2.SetPtEtaPhiM(mc_pt->at(10),mc_eta->at(10),mc_phi->at(10),mc_m->at(10));
108        //two quarks from top
109        vqt1.SetPtEtaPhiM(mc_pt->at(8),mc_eta->at(8),mc_phi->at(8),mc_m->at(8));
110        pTqt1 = mc_pt->at(8);
111        pdgqt1 = mc_pdgId->at(8);
112        vqt2.SetPtEtaPhiM(mc_pt->at(12),mc_eta->at(12),mc_phi->at(12),mc_m->at(12));
113        pTqt2 = mc_pt->at(12);
114        pdgqt2 = mc_pdgId->at(12);
115         
116        //loop on all particles to find the right W with status =2           
117        for (int i = 0; i < mc_n; i++){
118          if (mc_status->at(i)==2 && mc_pdgId->at(i)==mc_pdgId->at(7)){
119            std::vector<int> daughter = mc_child_index->at(i);
120            //Try to find the quark comes from quark 9 in the particles with stat=2 or 10092
121            dRmin = 1e9;
122            for (unsigned int id = 0; id < daughter.size(); id++) {
123              int pdgDaughter = mc_pdgId->at(daughter.at(id));
124              TLorentzVector vdaughter;
125              if ((mc_status->at(daughter.at(id))==2 || mc_status->at(daughter.at(id))==10092) && pdgDaughter==mc_pdgId->at(9)){
126                vdaughter.SetPtEtaPhiM(mc_pt->at(daughter.at(id)),mc_eta->at(daughter.at(id)),mc_phi->at(daughter.at(id)),mc_m->at(daughter.at(id)));
127                double dR = vq1.DeltaR(vdaughter);
128                if (dR < dRmin){
129                  dRmin = dR;
130                  vq11  = vdaughter;
131                }       
132              }
133            }
134   
135            //Try to find the quark comes from quark 10 in the particles with stat=2 or 10092
136            dRmin = 1e9;
137            for (unsigned int id = 0; id < daughter.size(); id++) {
138              int pdgDaughter = mc_pdgId->at(daughter.at(id));
139              TLorentzVector vdaughter;
140              if ((mc_status->at(daughter.at(id))==2 || mc_status->at(daughter.at(id))==10092) && pdgDaughter==mc_pdgId->at(10)){
141                vdaughter.SetPtEtaPhiM(mc_pt->at(daughter.at(id)),mc_eta->at(daughter.at(id)),mc_phi->at(daughter.at(id)),mc_m->at(daughter.at(id)));
142                double dR = vq2.DeltaR(vdaughter);
143                if (dR < dRmin){
144                  dRmin = dR;
145                  vq22  = vdaughter;
146                }       
147              }
148            }     
149          }//end of if (mc_status->at(i)==2 && mc_pdgId->at(i)==mc_pdgId->at(7))
150        }//end of loop for (int i = 0; i < mc_n; i++)
151      }//end of if (abs(mc_pdgId->at(9))<=5 && abs(mc_pdgId->at(10))<=5)
152       
153      //Check if the two particle 13 and 14 are quarks or not
154      if(abs(mc_pdgId->at(13))<=5 && abs(mc_pdgId->at(14))<=5){
155        vq1.SetPtEtaPhiM(mc_pt->at(13),mc_eta->at(13),mc_phi->at(13),mc_m->at(13));
156        vq2.SetPtEtaPhiM(mc_pt->at(14),mc_eta->at(14),mc_phi->at(14),mc_m->at(14));
157        //two quarks from top
158        vqt2.SetPtEtaPhiM(mc_pt->at(8),mc_eta->at(8),mc_phi->at(8),mc_m->at(8));
159        pTqt2 = mc_pt->at(8);
160        pdgqt2 = mc_pdgId->at(8);
161        vqt1.SetPtEtaPhiM(mc_pt->at(12),mc_eta->at(12),mc_phi->at(12),mc_m->at(12));
162        pTqt1 = mc_pt->at(12);
163        pdgqt1 = mc_pdgId->at(12);
164         
165        //loop on all particles to find the right W with status =2           
166        for (int i = 0; i < mc_n; i++){
167          if (mc_status->at(i)==2 && mc_pdgId->at(i)==mc_pdgId->at(11)) {
168            std::vector<int> daughter = mc_child_index->at(i);
169            //Try to find the quark comes from quark 13 in the particles with stat=2 or 10092
170            dRmin = 1e9;
171            for (unsigned int id = 0; id < daughter.size(); id++) {
172              int pdgDaughter = mc_pdgId->at(daughter.at(id));
173              TLorentzVector vdaughter;
174              if ((mc_status->at(daughter.at(id))==2 || mc_status->at(daughter.at(id))==10092) && pdgDaughter==mc_pdgId->at(13)){
175                vdaughter.SetPtEtaPhiM(mc_pt->at(daughter.at(id)),mc_eta->at(daughter.at(id)),mc_phi->at(daughter.at(id)),mc_m->at(daughter.at(id)));
176                double dR = vq1.DeltaR(vdaughter);
177                if (dR < dRmin){
178                  dRmin = dR;
179                  vq11  = vdaughter;
180                }       
181              }
182            }       
183             
184            //Try to find the quark comes from quark 14 in the particles with stat=2 or 10092
185            dRmin = 1e9;
186            for (unsigned int id = 0; id < daughter.size(); id++) {
187              int pdgDaughter = mc_pdgId->at(daughter.at(id));
188              TLorentzVector vdaughter;
189              if ((mc_status->at(daughter.at(id))==2 || mc_status->at(daughter.at(id))==10092) && pdgDaughter==mc_pdgId->at(14)){
190                vdaughter.SetPtEtaPhiM(mc_pt->at(daughter.at(id)),mc_eta->at(daughter.at(id)),mc_phi->at(daughter.at(id)),mc_m->at(daughter.at(id)));
191                double dR = vq2.DeltaR(vdaughter);
192                if (dR < dRmin){
193                  dRmin = dR;
194                  vq22  = vdaughter;
195                }       
196              }
197            }
198          }//end of (mc_status->at(i)==2 && mc_pdgId->at(i)==mc_pdgId->at(11))
199        }//end of for (int i = 0; i < mc_n; i++)
200      }//end of if(abs(mc_pdgId->at(13))<=5 && abs(mc_pdgId->at(14))<=5)       
201
202         
203      if ((!(vq11 == vq22))&&(!(vq11 == vec0))&&(!(vq22 == vec0))&&(!(vqt1==vec0))&&(!(vqt2==vec0))) {
204        //re-put return kTRUE
205        //if ((vq11 == vq22) || (vq11 == vec0) || (vq22 == vec0) || (vqt1==vec0) || (vqt2==vec0)) return kTRUE;
206        m_loose1->Fill(4);
207         
208        //Dexu quarks top
209        dRmin = 1e9;
210        for (int j = 0; j < mc_n; j++){
211          if ((mc_status->at(j)==2 || mc_status->at(j)==10092) && mc_pdgId->at(j)==pdgqt1){
212            double pTj = mc_pt->at(j);
213            vqj.SetPtEtaPhiM(pTj,mc_eta->at(j),mc_phi->at(j),mc_m->at(j));
214            double dR = vqt1.DeltaR(vqj);
215            if (dR < dRmin){
216              dRmin = dR;
217              vqt11  = vqj;
218            }
219          }
220        }
221         
222        dRmin = 1e9;
223        for (int j = 0; j < mc_n; j++){
224          if((mc_status->at(j)==2 || mc_status->at(j)==10092) && mc_pdgId->at(j)==pdgqt2){
225            double pTj = mc_pt->at(j);
226            vqj.SetPtEtaPhiM(pTj,mc_eta->at(j),mc_phi->at(j),mc_m->at(j));
227            double dR = vqt2.DeltaR(vqj);
228            if (dR < dRmin){
229              dRmin = dR;
230              vqt22  = vqj;
231            }
232          }
233        }
234
235        if ((!(vqt11 == vqt22))&&(!(vqt11 == vec0))&&(!(vqt22 == vec0))) {
236          //re-put return kTRUE
237          //if ( (vqt11 == vqt22) || (vqt11 == vec0) || (vqt22 == vec0)) return kTRUE;
238          m_loose1->Fill(5);
239
240          dRmin = 1e9;
241          //find the real electron in el_
242          for (unsigned int i = 0; i < el_cl_pt->size(); i++) {
243            if (el_truth_barcode->at(i) == barele) 
244              vje.SetPtEtaPhiM(el_cl_pt->at(i),el_cl_eta->at(i),el_cl_phi->at(i),0);
245            //if(el_author->at(i)!=1 && el_author->at(i)!=3) continue;
246            //TLorentzVector v4e;
247            //v4e.SetPtEtaPhiM(el_cl_pt->at(i),el_cl_eta->at(i),el_cl_phi->at(i),0);
248            //double dRej  = vele.DeltaR(v4e);
249            //if (dRej < 0.4 && dRej<dRmin){
250            //dRmin = dRej;
251            //vje = v4e;     
252            //}
253          }
254
255          if (vje==vec0) {//match electron status1 with el_ but not reject the not-matched-event in order to keep the same number of event for matching partons and jets=>just fill some histograms here
256            m_pTeletruth_noreco->Fill(vele.Pt()/1000);
257            m_etaeletruth_noreco->Fill(fabs(vele.Eta()));
258            //add a 2-dim pT-eta
259            m_pTeta_eletruth_noreco->Fill(fabs(vele.Eta()),vele.Pt()/1000);
260          }
261 
262          //Secondly, after finding the 2 true quarks W and top with stat = 2 or 10092
263          //we calculate the distribution of dR between these quarks and the jets
264          vector<TLorentzVector> jets1;
265          double dRmin1 = 1e9, dRmin2 = 1e9, dRmin3 = 1e9, dRmin4 = 1e9;
266          //Find the two jets closest to the quarks
267          for (unsigned int ij = 0; ij < jet_antikt4h1topo_pt->size(); ij++) {
268            TLorentzVector v4j;
269            v4j.SetPtEtaPhiM(jet_antikt4h1topo_pt->at(ij),jet_antikt4h1topo_eta->at(ij),jet_antikt4h1topo_phi->at(ij),jet_antikt4h1topo_m->at(ij));
270             
271            //if (jet_antikt4h1topo_pt->at(ij) < 20000) continue;
272            double dRq1j = vq11.DeltaR(v4j);
273            double dRq2j = vq22.DeltaR(v4j);
274            double dRqt1j = vqt11.DeltaR(v4j);
275            double dRqt2j = vqt22.DeltaR(v4j);
276               
277            //quarkW 1
278            if (dRq1j < 0.4 && dRq1j<dRmin1){
279              dRmin1 = dRq1j;
280              vj1 = v4j;//jet matched to quarkW 1
281              match[0] = ij;
282            }   
283            //quarkW 2
284            if (dRq2j < 0.4 && dRq2j<dRmin2){
285              dRmin2 = dRq2j;
286              vj2 = v4j;//jet matched to quarkW 1
287              match[1] = ij;
288            }
289            //quark top 1
290            if (dRqt1j < 0.4 && dRqt1j<dRmin3){
291              dRmin3 = dRqt1j;
292              vqj1 = v4j;//jet matched to quark top 1
293              match[2] = ij;
294            }
295            //quark top 2
296            if (dRqt2j < 0.4 && dRqt2j<dRmin4){
297              dRmin4 = dRqt2j;
298              vqj2 = v4j;//jet matched to quark top 2
299              match[3] = ij;
300            }
301          }//end of for (unsigned int ij = 0; ij < jet_antikt4h1topo_pt->size(); ij++)
302           
303          //int a = 0, b = 0;
304          for (int i=0; i<4; i++){
305            double x = match[i];
306            if (x==-1) allPartonMatch = false;//check if we cannot find the jet matched to some quark
307            for (int j=3; j>i; j--){
308              double y = match[j];
309              if (y==x) allMatchToDiff = false;//check if there are jets merged
310            } 
311          }           
312
313          if (allPartonMatch) {//Always find jets corresponding to partons     
314            //re-put return kTRUE
315            //if (a!=0) return kTRUE;//Always find jets corresponding to partons       
316            m_loose1->Fill(6);
317            //At this step, it can be that some partons are merged in the same jet
318            //Categories
319            //Category 1 : W merged
320            if (match[0]==match[1]) 
321              m_Wmerge->Fill(EventNumber);
322            //Category 2 : 1qW and 1qtop merged
323            if (match[0]==match[2] || match[0]==match[3] || match[1]==match[2] || match[1]==match[3]) 
324              m_qWqtopmerge->Fill(EventNumber);
325            //Category 3 : 2qtop merged
326            if (match[2]==match[3])
327              m_2qtopmerge->Fill(EventNumber);
328
329            if((!(vje==vec0)) && allMatchToDiff) {  //to choose the evts where we have good matching for electron and quarks
330              //re-put return kTRUE
331              //if((vje==vec0) || b!=0) return kTRUE;  //to choose the evts where we have good matching for electron and quarks
332              m_loose1->Fill(7);
333
334              matchedjets.push_back(vqj1);
335              matchedjets.push_back(vj1);
336              matchedjets.push_back(vj2);
337              matchedjets.push_back(vqj2);
338 
339              //This histo helps us to see if choosing 4 highest pT jets is a good idea or not
340              sort(matchedjets.begin(),matchedjets.end(),sortdecreasepT);
341              TLorentzVector vmmin = matchedjets.at(3); 
342              for (unsigned int ij = 0; ij < jet_antikt4h1topo_pt->size(); ij++) {
343                TLorentzVector v4j;
344                v4j.SetPtEtaPhiM(jet_antikt4h1topo_pt->at(ij),jet_antikt4h1topo_eta->at(ij),jet_antikt4h1topo_phi->at(ij),jet_antikt4h1topo_m->at(ij));
345                if ((!(v4j==vqj1)) && (!(v4j==vj1)) && (!(v4j==vj2)) && (!(v4j==vqj2))){
346                  double dpT = jet_antikt4h1topo_pt->at(ij) - vmmin.Pt();
347                  m_pTnotmatched_matchedmin->Fill(dpT/1000.);
348                }
349              }
350
351              if (vmmin.Pt()>2e4) 
352                m_loose1->Fill(10);
353               
354
355              double pTem = vje.Pt();
356              double pXem = vje.Vect().X();
357              double pYem = vje.Vect().Y();
358              double pZem = vje.Vect().Z();
359              double Eem  = vje.E();
360
361              double EmX = MET_RefFinal_etx;
362              double EmY = MET_RefFinal_ety; 
363              double EmT = MET_RefFinal_et;
364 
365              double km  = 80.4*80.4e6 + 2.*(pXem*EmX+pYem*EmY);
366              double am  = 4.*pTem*pTem;
367              double bm  = -4.*km*pZem;
368              double cm  = 4.*Eem*Eem*EmT*EmT-km*km;
369
370              double deltam = bm*bm - 4.*am*cm;
371
372              if (deltam < 0) deltam = 0;
373 
374              double pZmism1 = (-bm + sqrt(deltam))/2./am;
375              double pZmism2 = (-bm - sqrt(deltam))/2./am;
376              TVector3 pmism31(EmX, EmY, pZmism1);
377              TVector3 pmism32(EmX, EmY, pZmism2);
378
379              TLorentzVector pmism41(pmism31,pmism31.Mag());
380              TLorentzVector pmism42(pmism32,pmism32.Mag());
381
382              TLorentzVector vWlepmatch1 = vje + pmism41;
383              TLorentzVector vWlepmatch2 = vje + pmism42;
384
385              TLorentzVector vtoplepmatch1 = vWlepmatch1 + vqj2;
386              TLorentzVector vtoplepmatch2 = vWlepmatch2 + vqj2;
387
388              //Fill the histo for matching jets
389              TLorentzVector vWhadmatch = vj1 + vj2;
390              TLorentzVector vtophadmatch = vWhadmatch + vqj1;
391              TLorentzVector vttpairmatch;
392              double mWl = 0, mtl = 0, dmt = 0;
393              m_mWhadmatch->Fill(vWhadmatch.M()/1000.);
394              m_mtophadmatch->Fill(vtophadmatch.M()/1000.);
395              //Here to fill the histos of leptonic masses we choose the solution
396              //that minimize the mass difference of hadronic and leptonic top
397              if (fabs(vtoplepmatch1.M() - vtophadmatch.M())>fabs(vtoplepmatch2.M() - vtophadmatch.M())) {
398                mWl = vWlepmatch2.M()/1000;
399                mtl = vtoplepmatch2.M()/1000;
400                dmt = (vtoplepmatch2.M() - vtophadmatch.M())/1000.;
401                //m_mWlepmatch->Fill(vWlepmatch2.M()/1000);
402                //m_mtoplepmatch->Fill(vtoplepmatch2.M()/1000);
403                //m_dmtmatch->Fill((vtoplepmatch2.M() - vtophadmatch.M())/1000.);
404                vttpairmatch = vtophadmatch + vtoplepmatch2;
405              }
406              else {
407                mWl = vWlepmatch1.M()/1000;
408                mtl = vtoplepmatch1.M()/1000;
409                dmt = (vtoplepmatch1.M() - vtophadmatch.M())/1000.;
410                vttpairmatch = vtophadmatch + vtoplepmatch1;
411              }
412              m_mWlepmatch->Fill(mWl);
413              m_mtoplepmatch->Fill(mtl);
414              m_dmtmatch->Fill(dmt);
415              m_pTttpairmatch->Fill(vttpairmatch.Pt()/1000);
416              //fill the tree for matching
417              b_mWhadmatch = vWhadmatch.M()/1000.;
418              b_mWlepmatch = mWl;
419              b_mtophadmatch = vtophadmatch.M()/1000.;
420              b_mtoplepmatch = mtl;
421            }//end of loop for 1 reconstructed electron and 4 different matched-jets
422          }//end of loop for 4 matched-jets are different to vec0
423        }//end of loop for finding 2quarks top status 2 or 10092
424      }//end of loop for finding 2quarksW status 2 or 10092
425    }//end of loop if (!(vneu==vec0))
426  }//end of loop if (!(vele==vec0))
427}
Note: See TracBrowser for help on using the repository browser.