ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FirstData/Macros/ProjectQQ.C
Revision: 1.1
Committed: Thu Dec 17 18:09:59 2009 UTC (15 years, 4 months ago) by hwoehri
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 hwoehri 1.1 #define ProjectQQ_cxx
2     #include "ProjectQQ.h"
3     #include <TH2.h>
4     #include <TStyle.h>
5     #include <TCanvas.h>
6     #include <TLorentzVector.h>
7     #include "TClonesArray.h"
8    
9     #include "commonVar.inc"
10    
11     //dimuon related histograms
12     TH1F *hReco_mass[kNbCharge][kNbCath][kNbCuts];
13     TH1F *hReco_pT[kNbDimuonSet][kNbCharge][kNbCath][kNbCuts];
14     TH1F *hReco_rap[kNbDimuonSet][kNbCharge][kNbCath][kNbCuts];
15     TH1F *hRecoBest_mass[kNbCharge][kNbCath];
16     TH1F *hRecoBest_pT[kNbDimuonSet][kNbCharge][kNbCath];
17     TH1F *hRecoBest_rap[kNbDimuonSet][kNbCharge][kNbCath];
18     TH1F *hMult_QQ[kNbCharge][kNbCath][kNbCuts];
19     TH2F *hReco_eta1_eta2[kNbDimuonSet][kNbCharge][kNbCath+1][kNbCuts];
20     TH1F *hRecoNoBest;
21    
22     //single muon histos for two dimuon mass windows:
23     TH1F *hMuon_eta[kNbDimuonSet][kNbCharge][kNbCath];
24     TH1F *hMuon_pT[kNbDimuonSet][kNbCharge][kNbCath];
25     TH2F *hMuon_pT_eta[kNbDimuonSet][kNbCharge][kNbMuCath];
26     // quality histogams
27     TH1F *hChi2GlobalFit[kNbDimuonSet][kNbCharge];
28     TH1F *hChi2TrackerFit[kNbDimuonSet][kNbCharge][kNbMuCath];
29     TH1F *hTM2DCompatibilityTight[kNbDimuonSet][kNbCharge];
30     TH1F *hTMLastStationOptimizedLowPtLoose[kNbDimuonSet][kNbCharge];
31     TH1F *hTrackerMuonArbitrated[kNbDimuonSet][kNbCharge];//tracker only
32     TH1F *hCaloCompatibility[kNbDimuonSet][kNbCharge][kNbMuCath];
33     TH1F *hNHitsSilicon[kNbDimuonSet][kNbCharge][kNbMuCath];
34     TH1F *hDimuVtxProb[kNbDimuonSet][kNbCharge][kNbCath];
35    
36     //=======================================
37     void ProjectQQ::Loop(Bool_t removeQQ, Bool_t matchMC)
38     {
39     if (fChain == 0) return;
40    
41     Long64_t nentries = fChain->GetEntries();
42     // nentries = 10000;
43    
44     Long64_t nb = 0;
45    
46     Int_t countAllRecoQQEvents = 0, countNotFoundBestQQEvents = 0;
47     for (Long64_t jentry=0; jentry<nentries;jentry++) {
48    
49     if(jentry % 1000 == 0) printf("event %d/%d\n", (Int_t) jentry, (Int_t) nentries);
50    
51     Long64_t ientry = LoadTree(jentry);
52     if (ientry < 0) break;
53     nb = fChain->GetEntry(jentry);
54    
55     if(removeQQ && Mc_QQ_size > 0) continue; //Mc_QQ_size contains the nb. of MC J/psi's
56    
57     Int_t chargeID;
58     Int_t nComb[kNbCharge][kNbCath][kNbCuts];
59     for(int iCharge = 0; iCharge < kNbCharge; iCharge++)
60     for(int iCat = 0; iCat < kNbCath; iCat++)
61     for(int iCut = 0; iCut < kNbCuts; iCut++)
62     nComb[iCharge][iCat][iCut] = 0;
63    
64     //loop over all quarkonia
65     for(int iQQ = 0; iQQ < Reco_QQ_size; iQQ++){
66    
67     chargeID = 0;
68     if(Reco_QQ_sign[iQQ] < 0) chargeID = 1;
69     else if(Reco_QQ_sign[iQQ] > 0) chargeID = 2;
70    
71     TLorentzVector *Reco_QQ = (TLorentzVector *) Reco_QQ_4mom->At(iQQ);
72     Double_t mass_Reco_QQ = Reco_QQ->M();
73     Double_t pT_Reco_QQ = Reco_QQ->Pt();
74     Double_t rap_Reco_QQ = Reco_QQ->Rapidity();
75    
76     //get the individual muons:
77     TLorentzVector *Reco_QQ_posMu, *Reco_QQ_negMu;
78     Bool_t etaCorr = kFALSE;
79     if(Reco_QQ_type[iQQ] == 0){ //global-global
80     if(Reco_QQ_mupl[iQQ] < Reco_mu_glb_size && Reco_QQ_mumi[iQQ] < Reco_mu_glb_size){
81     Reco_QQ_posMu = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_mupl[iQQ]);
82     Reco_QQ_negMu = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_mumi[iQQ]);
83     etaCorr = kTRUE;
84     }
85     }
86     else if(Reco_QQ_type[iQQ] == 2){ //tracker-tracker
87     if(Reco_QQ_mupl[iQQ] < Reco_mu_trk_size && Reco_QQ_mumi[iQQ] < Reco_mu_trk_size){
88     Reco_QQ_posMu = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_mupl[iQQ]);
89     Reco_QQ_negMu = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_mumi[iQQ]);
90     etaCorr = kTRUE;
91     }
92     }
93    
94     Bool_t isMatched[2];//only relevant for MC
95     for (int iCut=0; iCut<kNbCuts; iCut++)
96     {
97     //start with the MC matching:
98     // if (iCut == 0 && matchMC){
99     if (matchMC){
100    
101     if(Mc_QQmupl_indx[0] > 1 || Mc_QQmumi_indx[0] > 1) continue; //sanity check
102    
103     isMatched[0] = kFALSE;
104     isMatched[1] = kFALSE;
105    
106     switch (Reco_QQ_type[iQQ]) {
107     case 0: // global-global
108     //first MC muon
109     if(Reco_QQ_muhpt[iQQ] < Reco_mu_glb_size){//sanity check
110     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
111     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
112     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
113     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
114     isMatched[0] = kTRUE;
115     }
116     //second MC muon
117     if(Reco_QQ_mulpt[iQQ] < Reco_mu_glb_size){//sanity check
118     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
119     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
120     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
121     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
122     isMatched[1] = kTRUE;
123     }
124     break;
125     case 1: // global-tracker
126     //first MC muon
127     if(Reco_QQ_muhpt[iQQ] < Reco_mu_glb_size){//sanity check
128     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
129     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
130     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
131     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
132     isMatched[0] = kTRUE;
133     }
134     //second MC muon
135     if(Reco_QQ_mulpt[iQQ] < Reco_mu_trk_size){//sanity check
136     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
137     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
138     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
139     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
140     isMatched[1] = kTRUE;
141     }
142     break;
143     case 2: // tracker-tracker
144     //first MC muon
145     if(Reco_QQ_muhpt[iQQ] < Reco_mu_trk_size){//sanity check
146     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
147     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
148     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
149     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
150     isMatched[0] = kTRUE;
151     }
152     //second MC muon
153     if(Reco_QQ_mulpt[iQQ] < Reco_mu_trk_size){//sanity check
154     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
155     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
156     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
157     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
158     isMatched[1] = kTRUE;
159     }
160     break;
161     case 3: // global-calo
162     //first MC muon
163     if(Reco_QQ_muhpt[iQQ] < Reco_mu_glb_size){//sanity check
164     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
165     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
166     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
167     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
168     isMatched[0] = kTRUE;
169     }
170     //second MC muon
171     if(Reco_QQ_mulpt[iQQ] < Reco_mu_cal_size){//sanity check
172     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
173     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
174     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
175     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
176     isMatched[1] = kTRUE;
177     }
178     break;
179     case 4: // tracker-calo
180     //first MC muon
181     if(Reco_QQ_muhpt[iQQ] < Reco_mu_trk_size){//sanity check
182     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
183     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
184     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
185     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
186     isMatched[0] = kTRUE;
187     }
188     //second MC muon
189     if(Reco_QQ_mulpt[iQQ] < Reco_mu_cal_size){//sanity check
190     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
191     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
192     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
193     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
194     isMatched[1] = kTRUE;
195     }
196     break;
197     case 5: // calo-calo
198     //first MC muon
199     if(Reco_QQ_muhpt[iQQ] < Reco_mu_cal_size){//sanity check
200     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
201     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR ||
202     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
203     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_muhpt[iQQ])) < MAX_deltaR)
204     isMatched[0] = kTRUE;
205     }
206     //second MC muon
207     if(Reco_QQ_mulpt[iQQ] < Reco_mu_cal_size){//sanity check
208     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
209     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR ||
210     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
211     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iQQ])) < MAX_deltaR)
212     isMatched[1] = kTRUE;
213     }
214     break;
215     default:
216     break;
217     }
218     //check the MC matching before continuing...
219     if(isMatched[0] == kFALSE || isMatched[1] == kFALSE)
220     continue;
221     }
222    
223     if (iCut>0) {
224     // distance cut
225     switch (Reco_QQ_type[iQQ]) {
226     case 0: // global-global
227     if ( !(fabs(Reco_mu_glb_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
228     fabs(Reco_mu_glb_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
229     fabs(Reco_mu_glb_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
230     fabs(Reco_mu_glb_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
231     continue;
232     break;
233     case 1: // global-tracker
234     if ( !(fabs(Reco_mu_glb_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
235     fabs(Reco_mu_glb_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
236     fabs(Reco_mu_trk_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
237     fabs(Reco_mu_trk_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
238     continue;
239     break;
240     case 2: // tracker-tracker
241     if ( !(fabs(Reco_mu_trk_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
242     fabs(Reco_mu_trk_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
243     fabs(Reco_mu_trk_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
244     fabs(Reco_mu_trk_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
245     continue;
246     break;
247     case 3: // global-calo
248     if ( !(fabs(Reco_mu_glb_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
249     fabs(Reco_mu_glb_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
250     fabs(Reco_mu_cal_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
251     fabs(Reco_mu_cal_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
252     continue;
253     break;
254     case 4: // tracker-calo
255     if ( !(fabs(Reco_mu_trk_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
256     fabs(Reco_mu_trk_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
257     fabs(Reco_mu_cal_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
258     fabs(Reco_mu_cal_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
259     continue;
260     break;
261     case 5: // calo-calo
262     if ( !(fabs(Reco_mu_cal_d0[Reco_QQ_muhpt[iQQ]]) < MAX_d0_trk &&
263     fabs(Reco_mu_cal_dz[Reco_QQ_muhpt[iQQ]]) < MAX_dz_trk &&
264     fabs(Reco_mu_cal_d0[Reco_QQ_mulpt[iQQ]]) < MAX_d0_trk &&
265     fabs(Reco_mu_cal_dz[Reco_QQ_mulpt[iQQ]]) < MAX_dz_trk) )
266     continue;
267     break;
268     default:
269     break;
270     }
271     }
272    
273     if (iCut>1) {
274     // for tracker muons: TM2DCompatibilityTight || TMLastStationOptimizedLowPtTight
275     // for calo muons: calo compatibility
276     switch (Reco_QQ_type[iQQ]) {
277     case 0: // global-global
278     break;
279     case 1: // global-tracker
280     if ((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5) == 0 &&
281     (Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8) == 0)
282     continue;
283     break;
284     case 2: // tracker-tracker
285     if ( ((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5) == 0 &&
286     (Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8) == 0) ||
287     ((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5) == 0 &&
288     (Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8) == 0) )
289     continue;
290     break;
291     case 3: // global-calo
292     if ( !(Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]] > MIN_caloComp) )
293     continue;
294     break;
295     case 4: // tracker-calo
296     if ( ((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5) == 0 &&
297     (Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8) == 0) ||
298     Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]] < MIN_caloComp)
299     continue;
300     break;
301     case 5: // calo-calo
302     if ( !(Reco_mu_cal_caloComp[Reco_QQ_muhpt[iQQ]] > MIN_caloComp &&
303     Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]] > MIN_caloComp) )
304     continue;
305     break;
306     default:
307     break;
308     }
309     }
310    
311     if (iCut>2) {
312     // for global: chi2/ndf (global fit)
313     // for tracker and calo: chi2/ndf (track fit) && #hits
314     switch (Reco_QQ_type[iQQ]) {
315     case 0: // global-global
316     if ( !(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_glb &&
317     Reco_mu_glb_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_glb) )
318     continue;
319     break;
320     case 1: // global-tracker
321     if ( !(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_glb &&
322     Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_trk &&
323     Reco_mu_trk_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) )
324     continue;
325     break;
326     case 2: // tracker-tracker
327     if ( !(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_trk &&
328     Reco_mu_trk_nhitstrack[Reco_QQ_muhpt[iQQ]] > MIN_nhits_trk &&
329     Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_trk &&
330     Reco_mu_trk_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) )
331     continue;
332     break;
333     case 3: // global-calo
334     if ( !(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_glb &&
335     Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
336     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) )
337     continue;
338     break;
339     case 4: // tracker-calo
340     if ( !(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_trk &&
341     Reco_mu_trk_nhitstrack[Reco_QQ_muhpt[iQQ]] > MIN_nhits_trk &&
342     Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
343     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) )
344     continue;
345     break;
346     case 5: // calo-calo
347     if ( !(Reco_mu_cal_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_cal &&
348     Reco_mu_cal_nhitstrack[Reco_QQ_muhpt[iQQ]] > MIN_nhits_trk &&
349     Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
350     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) )
351     continue;
352     break;
353     default:
354     break;
355     }
356     }
357    
358     if (iCut>3) {
359     // dimuon vertex probability
360     if (Reco_QQ_probChi2[iQQ]<= MIN_vtxprob_jpsi)
361     continue;
362     }
363    
364     nComb[chargeID][Reco_QQ_type[iQQ]][iCut]++;
365    
366     //mass histograms
367     hReco_mass[chargeID][Reco_QQ_type[iQQ]][iCut]->Fill(mass_Reco_QQ);
368     if(chargeID == 1 || chargeID == 2) //LS histo
369     hReco_mass[3][Reco_QQ_type[iQQ]][iCut]->Fill(mass_Reco_QQ);
370    
371     //fill pT and rap histos only for the signal:
372     for(int iMassW = 0; iMassW < kNbDimuonSet; iMassW++){
373     if(mass_Reco_QQ > massMIN[iMassW] && mass_Reco_QQ < massMAX[iMassW]){
374     //pT histograms
375     hReco_pT[iMassW][chargeID][Reco_QQ_type[iQQ]][iCut]->Fill(pT_Reco_QQ);
376     if(chargeID == 1 || chargeID == 2) //LS histo
377     hReco_pT[iMassW][3][Reco_QQ_type[iQQ]][iCut]->Fill(pT_Reco_QQ);
378     //rap histograms
379     hReco_rap[iMassW][chargeID][Reco_QQ_type[iQQ]][iCut]->Fill(rap_Reco_QQ);
380     if(chargeID == 1 || chargeID == 2) //LS histo
381     hReco_rap[iMassW][3][Reco_QQ_type[iQQ]][iCut]->Fill(rap_Reco_QQ);
382    
383     //correlating the single muons:
384     if(etaCorr)
385     hReco_eta1_eta2[iMassW][chargeID][Reco_QQ_type[iQQ]][iCut]->Fill(Reco_QQ_posMu->Eta(),Reco_QQ_negMu->Eta());
386    
387     // prepare single muon indices to avoid double counting
388     int glMuID[10000], trMuID[10000],caMuID[10000];
389     for (int i=0; i<10000; i++){
390     glMuID[i] = -1;
391     trMuID[i] = -1;
392     caMuID[i] = -1;
393     }
394    
395     // fill single muon quantities
396     // if (iCut==0) { // d0 and dz
397     // switch (Reco_QQ_type[iQQ]) {
398     // case 0: // global+global
399     // hChi2GlobalFit[0][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
400     // hChi2GlobalFit[0][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_mulpt[iQQ]]);
401     // break;
402     // case 1: // global+tracker
403     // hChi2GlobalFit[0][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
404     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]]);
405     // break;
406     // case 2: // tracker+tracker
407     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]]);
408     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]]);
409     // break;
410     // case 3: // global+calo
411     // hChi2GlobalFit[0][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
412     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
413     // break;
414     // case 4: // tracker+calo
415     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]]);
416     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
417     // break;
418     // case 5: // calo+calo
419     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_cal_normChi2[Reco_QQ_muhpt[iQQ]]);
420     // hChi2TrackerFit[0][chargeID]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
421     // break;
422     // default:
423     // break;
424     // }
425     // }
426    
427     if (iCut==1) { // comp, after d0 and dz
428     switch (Reco_QQ_type[iQQ]) {
429     case 0: // global+global
430     break;
431     case 1: // global+tracker
432     if (trMuID[Reco_QQ_mulpt[iQQ]]!=1) {
433     hTM2DCompatibilityTight[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5));
434     hTMLastStationOptimizedLowPtLoose[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8));
435     hTrackerMuonArbitrated[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int)pow(2,1))/(int)pow(2,1));
436     hCaloCompatibility[iMassW][chargeID][1]->Fill(Reco_mu_trk_caloComp[Reco_QQ_mulpt[iQQ]]);
437     trMuID[Reco_QQ_mulpt[iQQ]]=1;
438     }
439     break;
440     case 2: // tracker+tracker
441     if (trMuID[Reco_QQ_muhpt[iQQ]]!=1) {
442     hTM2DCompatibilityTight[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5));
443     hTMLastStationOptimizedLowPtLoose[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8));
444     hTrackerMuonArbitrated[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int)pow(2,1))/(int)pow(2,1));
445    
446     hCaloCompatibility[iMassW][chargeID][1]->Fill(Reco_mu_trk_caloComp[Reco_QQ_muhpt[iQQ]]);
447     trMuID[Reco_QQ_muhpt[iQQ]]=1;
448     }
449     if (trMuID[Reco_QQ_mulpt[iQQ]]!=1) {
450     hTM2DCompatibilityTight[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5));
451     hTMLastStationOptimizedLowPtLoose[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8));
452     hTrackerMuonArbitrated[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_mulpt[iQQ]] & (int)pow(2,1))/(int)pow(2,1));
453    
454     hCaloCompatibility[iMassW][chargeID][1]->Fill(Reco_mu_trk_caloComp[Reco_QQ_mulpt[iQQ]]);
455     trMuID[Reco_QQ_mulpt[iQQ]]=1;
456     }
457     break;
458     case 3: // global+calo
459     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
460     if (Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
461     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk)
462     hCaloCompatibility[iMassW][chargeID][2]->Fill(Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]]);
463     caMuID[Reco_QQ_mulpt[iQQ]]=1;
464     }
465     break;
466     case 4: // tracker+calo
467     if (trMuID[Reco_QQ_muhpt[iQQ]]!=1) {
468     hTM2DCompatibilityTight[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,5))/(int) pow(2,5));
469     hTMLastStationOptimizedLowPtLoose[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int) pow(2,8))/(int) pow(2,8));
470     hTrackerMuonArbitrated[iMassW][chargeID]->Fill((Reco_mu_trk_PIDmask[Reco_QQ_muhpt[iQQ]] & (int)pow(2,1))/(int)pow(2,1));
471     hCaloCompatibility[iMassW][chargeID][1]->Fill(Reco_mu_trk_caloComp[Reco_QQ_muhpt[iQQ]]);
472     trMuID[Reco_QQ_muhpt[iQQ]]=1;
473     }
474     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
475     if (Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
476     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk)
477     hCaloCompatibility[iMassW][chargeID][2]->Fill(Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]]);
478     caMuID[Reco_QQ_mulpt[iQQ]]=1;
479     }
480     break;
481     case 5: // calo+calo
482     if (caMuID[Reco_QQ_muhpt[iQQ]]!=1) {
483     if (Reco_mu_cal_normChi2[Reco_QQ_muhpt[iQQ]] < MAX_normchi2_cal &&
484     Reco_mu_cal_nhitstrack[Reco_QQ_muhpt[iQQ]] > MIN_nhits_trk) {
485     hCaloCompatibility[iMassW][chargeID][2]->Fill(Reco_mu_cal_caloComp[Reco_QQ_muhpt[iQQ]]);
486     caMuID[Reco_QQ_muhpt[iQQ]]=1;
487     }
488     }
489     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
490     if (Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]] < MAX_normchi2_cal &&
491     Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]] > MIN_nhits_trk) {
492     hCaloCompatibility[iMassW][chargeID][2]->Fill(Reco_mu_cal_caloComp[Reco_QQ_mulpt[iQQ]]);
493     caMuID[Reco_QQ_mulpt[iQQ]]=1;
494     }
495     }
496     break;
497     default:
498     break;
499     }
500     }
501    
502     if (iCut==2) { // chi2, after comp
503     switch (Reco_QQ_type[iQQ]) {
504     case 0: // global+global
505     if (glMuID[Reco_QQ_muhpt[iQQ]]!=1) {
506     hChi2GlobalFit[iMassW][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
507     hNHitsSilicon[iMassW][chargeID][0]->Fill(Reco_mu_glb_nhitstrack[Reco_QQ_muhpt[iQQ]]);
508     glMuID[Reco_QQ_muhpt[iQQ]]=1;
509     }
510     if (glMuID[Reco_QQ_mulpt[iQQ]]!=1) {
511     hChi2GlobalFit[iMassW][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_mulpt[iQQ]]);
512     hNHitsSilicon[iMassW][chargeID][0]->Fill(Reco_mu_glb_nhitstrack[Reco_QQ_mulpt[iQQ]]);
513     glMuID[Reco_QQ_mulpt[iQQ]]=1;
514     }
515     break;
516     case 1: // global+tracker
517     if (glMuID[Reco_QQ_muhpt[iQQ]]!=1) {
518     hChi2GlobalFit[iMassW][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
519     hNHitsSilicon[iMassW][chargeID][0]->Fill(Reco_mu_glb_nhitstrack[Reco_QQ_muhpt[iQQ]]);
520     glMuID[Reco_QQ_mulpt[iQQ]]=1;
521     }
522     if (trMuID[Reco_QQ_mulpt[iQQ]]!=1) {
523     hChi2TrackerFit[iMassW][chargeID][1]->Fill(Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]]);
524     hNHitsSilicon[iMassW][chargeID][1]->Fill(Reco_mu_trk_nhitstrack[Reco_QQ_mulpt[iQQ]]);
525     trMuID[Reco_QQ_mulpt[iQQ]]=1;
526     }
527     break;
528     case 2: // tracker+tracker
529     if (trMuID[Reco_QQ_muhpt[iQQ]]!=1) {
530     hChi2TrackerFit[iMassW][chargeID][1]->Fill(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]]);
531     hNHitsSilicon[iMassW][chargeID][1]->Fill(Reco_mu_trk_nhitstrack[Reco_QQ_muhpt[iQQ]]);
532     trMuID[Reco_QQ_muhpt[iQQ]]=1;
533     }
534     if (trMuID[Reco_QQ_mulpt[iQQ]]!=1) {
535     hChi2TrackerFit[iMassW][chargeID][1]->Fill(Reco_mu_trk_normChi2[Reco_QQ_mulpt[iQQ]]);
536     hNHitsSilicon[iMassW][chargeID][1]->Fill(Reco_mu_trk_nhitstrack[Reco_QQ_mulpt[iQQ]]);
537     trMuID[Reco_QQ_mulpt[iQQ]]=1;
538     }
539     break;
540     case 3: // global+calo
541     if (glMuID[Reco_QQ_muhpt[iQQ]]!=1) {
542     hChi2GlobalFit[iMassW][chargeID]->Fill(Reco_mu_glb_normChi2[Reco_QQ_muhpt[iQQ]]);
543     hNHitsSilicon[iMassW][chargeID][0]->Fill(Reco_mu_glb_nhitstrack[Reco_QQ_muhpt[iQQ]]);
544     glMuID[Reco_QQ_muhpt[iQQ]]=1;
545     }
546     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
547     hChi2TrackerFit[iMassW][chargeID][2]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
548     hNHitsSilicon[iMassW][chargeID][2]->Fill(Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]]);
549     caMuID[Reco_QQ_mulpt[iQQ]]=1;
550     }
551     break;
552     case 4: // tracker+calo
553     if (trMuID[Reco_QQ_muhpt[iQQ]]!=1) {
554     hChi2TrackerFit[iMassW][chargeID][1]->Fill(Reco_mu_trk_normChi2[Reco_QQ_muhpt[iQQ]]);
555     hNHitsSilicon[iMassW][chargeID][1]->Fill(Reco_mu_trk_nhitstrack[Reco_QQ_muhpt[iQQ]]);
556     trMuID[Reco_QQ_muhpt[iQQ]]=1;
557     }
558     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
559     hChi2TrackerFit[iMassW][chargeID][2]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
560     hNHitsSilicon[iMassW][chargeID][2]->Fill(Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]]);
561     caMuID[Reco_QQ_mulpt[iQQ]]=1;
562     }
563     break;
564     case 5: // calo+calo
565     if (caMuID[Reco_QQ_muhpt[iQQ]]!=1) {
566     hChi2TrackerFit[iMassW][chargeID][2]->Fill(Reco_mu_cal_normChi2[Reco_QQ_muhpt[iQQ]]);
567     hNHitsSilicon[iMassW][chargeID][2]->Fill(Reco_mu_cal_nhitstrack[Reco_QQ_muhpt[iQQ]]);
568     caMuID[Reco_QQ_muhpt[iQQ]]=1;
569     }
570     if (caMuID[Reco_QQ_mulpt[iQQ]]!=1) {
571     hChi2TrackerFit[iMassW][chargeID][2]->Fill(Reco_mu_cal_normChi2[Reco_QQ_mulpt[iQQ]]);
572     hNHitsSilicon[iMassW][chargeID][2]->Fill(Reco_mu_cal_nhitstrack[Reco_QQ_mulpt[iQQ]]);
573     caMuID[Reco_QQ_mulpt[iQQ]]=1;
574     }
575     break;
576     default:
577     break;
578     }
579     }
580    
581     if (iCut==3) { // dimuon vertex prob, after comp
582     hDimuVtxProb[iMassW][chargeID][Reco_QQ_type[iQQ]]->Fill(Reco_QQ_probChi2[iQQ]);
583     if(chargeID == 1 || chargeID == 2) //LS histo
584     hDimuVtxProb[iMassW][3][Reco_QQ_type[iQQ]]->Fill(Reco_QQ_probChi2[iQQ]);
585     }
586     } // fill for signal
587     }//different mass windows
588     }//cuts
589     }//#Reco_QQ combinations
590    
591     for(int iCharge = 0; iCharge < kNbCharge; iCharge++)
592     for(int iCat = 0; iCat < kNbCath; iCat++)
593     for(int iCut = 0; iCut < kNbCuts; iCut++)
594     if(nComb[iCharge][iCat][iCut] > 0)
595     hMult_QQ[chargeID][iCat][iCut]->Fill(nComb[iCharge][iCat][iCut]);
596     for(int iCat = 0; iCat < kNbCath; iCat++)
597     for(int iCut = 0; iCut < kNbCuts; iCut++)
598     hMult_QQ[3][iCat][iCut]->Fill(nComb[1][iCat][iCut]+nComb[2][iCat][iCut]);
599    
600     //now fill the histos for the "best" QQ pair only
601     Int_t iDBestQQ = -1;
602     if(Reco_QQ_size > 0){
603     countAllRecoQQEvents++;
604     iDBestQQ = theBestQQ();
605     // if(iDBestQQ >= 0){
606     // printf("best QQ is %d (out of %d) and is of type %s\n",
607     // iDBestQQ, Reco_QQ_size, oniaCatName[Reco_QQ_type[iDBestQQ]]);
608     // }
609     // else{
610     // printf("best QQ not identified... %d total QQbar\n", Reco_QQ_size);
611     // countNotFoundBestQQEvents++;
612     // }
613     if(iDBestQQ < 0){
614     countNotFoundBestQQEvents++;
615     for(int iRec = 0; iRec < Reco_QQ_size; iRec++){
616    
617     if(Reco_QQ_sign[iRec] != 0) continue;
618     //we count ALL reco QQpairs in the event:
619     hRecoNoBest->Fill(Reco_QQ_type[iRec]);
620     }
621     }
622     else{
623    
624     TLorentzVector *Reco_QQ = (TLorentzVector *) Reco_QQ_4mom->At(iDBestQQ);
625     Double_t mass_Reco_QQ = Reco_QQ->M();
626     Double_t pT_Reco_QQ = Reco_QQ->Pt();
627     Double_t rap_Reco_QQ = Reco_QQ->Rapidity();
628    
629     chargeID = 0;
630     if(Reco_QQ_sign[iDBestQQ] < 0) chargeID = 1;
631     else if(Reco_QQ_sign[iDBestQQ] > 0) chargeID = 2;
632    
633     Bool_t isMatched[2] = {kTRUE, kTRUE};//only relevant for MC; for data must be set to kTRUE
634     if (matchMC){
635    
636     isMatched[0] = kFALSE;
637     isMatched[1] = kFALSE;
638    
639     if(Mc_QQmupl_indx[0] < 2 && Mc_QQmumi_indx[0] < 2){ //sanity check
640    
641     switch (Reco_QQ_type[iDBestQQ]) {
642     case 0: // global-global
643     //first MC muon
644     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_glb_size){//sanity check
645     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
646     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
647     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
648     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
649     isMatched[0] = kTRUE;
650     }
651     //second MC muon
652     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_glb_size){//sanity check
653     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
654     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
655     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
656     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
657     isMatched[1] = kTRUE;
658     }
659     break;
660     case 1: // global-tracker
661     //first MC muon
662     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_glb_size){//sanity check
663     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
664     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
665     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
666     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
667     isMatched[0] = kTRUE;
668     }
669     //second MC muon
670     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_trk_size){//sanity check
671     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
672     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
673     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
674     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
675     isMatched[1] = kTRUE;
676     }
677     break;
678     case 2: // tracker-tracker
679     //first MC muon
680     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_trk_size){//sanity check
681     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
682     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
683     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
684     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
685     isMatched[0] = kTRUE;
686     }
687     //second MC muon
688     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_trk_size){//sanity check
689     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
690     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
691     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
692     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
693     isMatched[1] = kTRUE;
694     }
695     break;
696     case 3: // global-calo
697     //first MC muon
698     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_glb_size){//sanity check
699     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
700     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
701     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
702     (TLorentzVector*)Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
703     isMatched[0] = kTRUE;
704     }
705     //second MC muon
706     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_cal_size){//sanity check
707     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
708     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
709     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
710     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
711     isMatched[1] = kTRUE;
712     }
713     break;
714     case 4: // tracker-calo
715     //first MC muon
716     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_trk_size){//sanity check
717     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
718     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
719     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
720     (TLorentzVector*)Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
721     isMatched[0] = kTRUE;
722     }
723     //second MC muon
724     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_cal_size){//sanity check
725     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
726     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
727     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
728     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
729     isMatched[1] = kTRUE;
730     }
731     break;
732     case 5: // calo-calo
733     //first MC muon
734     if(Reco_QQ_muhpt[iDBestQQ] < Reco_mu_cal_size){//sanity check
735     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
736     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR ||
737     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
738     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_muhpt[iDBestQQ])) < MAX_deltaR)
739     isMatched[0] = kTRUE;
740     }
741     //second MC muon
742     if(Reco_QQ_mulpt[iDBestQQ] < Reco_mu_cal_size){//sanity check
743     if(deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmupl_indx[0]),
744     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR ||
745     deltaR((TLorentzVector*)Mc_mu_4mom->At(Mc_QQmumi_indx[0]),
746     (TLorentzVector*)Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ])) < MAX_deltaR)
747     isMatched[1] = kTRUE;
748     }
749     break;
750     default:
751     break;
752     }
753     }
754     }
755     //check the MC matching before continuing...
756     if(isMatched[0] == kTRUE && isMatched[1] == kTRUE){
757    
758     //mass histograms
759     hRecoBest_mass[chargeID][Reco_QQ_type[iDBestQQ]]->Fill(mass_Reco_QQ);
760     if(chargeID == 1 || chargeID == 2) //LS histo
761     hRecoBest_mass[3][Reco_QQ_type[iDBestQQ]]->Fill(mass_Reco_QQ);
762    
763     //fill pT and rap histos only for the signal:
764     for(int iMassW = 0; iMassW < kNbDimuonSet; iMassW++){
765     if(mass_Reco_QQ > massMIN[iMassW] && mass_Reco_QQ < massMAX[iMassW]){
766     //pT histograms
767     hRecoBest_pT[iMassW][chargeID][Reco_QQ_type[iDBestQQ]]->Fill(pT_Reco_QQ);
768     if(chargeID == 1 || chargeID == 2) //LS histo
769     hRecoBest_pT[iMassW][3][Reco_QQ_type[iDBestQQ]]->Fill(pT_Reco_QQ);
770     //rap histograms
771     hRecoBest_rap[iMassW][chargeID][Reco_QQ_type[iDBestQQ]]->Fill(rap_Reco_QQ);
772     if(chargeID == 1 || chargeID == 2) //LS histo
773     hRecoBest_rap[iMassW][3][Reco_QQ_type[iDBestQQ]]->Fill(rap_Reco_QQ);
774     }
775     //===============================================================
776     //single muon histograms for best dimuon within 3.0 < M < 3.2 GeV
777     //and for dimuons with M > 2 GeV
778     //===============================================================
779     TLorentzVector *Reco_mu_1, *Reco_mu_2;
780     Int_t iMuCat[2] = {-1,-1};
781     if(Reco_QQ_type[iDBestQQ] == 0){//gl-gl
782     Reco_mu_1 = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
783     Reco_mu_2 = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
784     iMuCat[0] = 0; iMuCat[1] = 0;
785     }
786     else if(Reco_QQ_type[iDBestQQ] == 1){//gl-tr
787     Reco_mu_1 = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
788     Reco_mu_2 = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
789     iMuCat[0] = 0; iMuCat[1] = 1;
790     }
791     else if(Reco_QQ_type[iDBestQQ] == 2){//tr-tr
792     Reco_mu_1 = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
793     Reco_mu_2 = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
794     iMuCat[0] = 1; iMuCat[1] = 1;
795     }
796     else if(Reco_QQ_type[iDBestQQ] == 3){//gl-calo
797     Reco_mu_1 = (TLorentzVector *) Reco_mu_glb_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
798     Reco_mu_2 = (TLorentzVector *) Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
799     iMuCat[0] = 0; iMuCat[1] = 2;
800     }
801     else if(Reco_QQ_type[iDBestQQ] == 4){//tr-calo
802     Reco_mu_1 = (TLorentzVector *) Reco_mu_trk_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
803     Reco_mu_2 = (TLorentzVector *) Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
804     iMuCat[0] = 1; iMuCat[1] = 2;
805     }
806     else if(Reco_QQ_type[iDBestQQ] == 5){//calo-calo
807     Reco_mu_1 = (TLorentzVector *) Reco_mu_cal_4mom->At(Reco_QQ_muhpt[iDBestQQ]);
808     Reco_mu_2 = (TLorentzVector *) Reco_mu_cal_4mom->At(Reco_QQ_mulpt[iDBestQQ]);
809     iMuCat[0] = 2; iMuCat[1] = 2;
810     }
811    
812     Double_t etaMu[2] = {Reco_mu_1->Eta(), Reco_mu_2->Eta()};
813     Double_t pTMu[2] = {Reco_mu_1->Pt(), Reco_mu_2->Pt()};
814     for(int iMu = 0; iMu < 2; iMu++){
815     hMuon_eta[iMassW][chargeID][iMuCat[iMu]]->Fill(etaMu[iMu]);
816     hMuon_pT[iMassW][chargeID][iMuCat[iMu]]->Fill(pTMu[iMu]);
817     hMuon_pT_eta[iMassW][chargeID][iMuCat[iMu]]->Fill(etaMu[iMu], pTMu[iMu]);
818     }
819     }//mass windows
820     }//matched: relevant for MC
821     }//best QQ
822     }//Reco_QQ_size > 0
823     }//events
824    
825     if(countAllRecoQQEvents > 0){
826     hRecoNoBest->Scale(1./countAllRecoQQEvents);
827     printf("\n\n The number of events w/o identification of the best QQbar is %d out of %d, i.e. %1.1f %%\n", countNotFoundBestQQEvents, countAllRecoQQEvents, 100.*countNotFoundBestQQEvents/countAllRecoQQEvents);
828     }
829     }
830    
831     //=============================
832     int ProjectQQ::theBestQQ() {
833    
834     int theBest = -1;
835     float thehighestPt = -1.;
836    
837     for (int iqq=0; iqq<Reco_QQ_size; iqq++) {
838    
839     if (Reco_QQ_sign[iqq] == 0 && Reco_QQ_type[iqq] == 0 ) {
840    
841     int thehptMu = Reco_QQ_muhpt[iqq]; if (thehptMu >= Reco_mu_glb_size) continue;
842     int thelptMu = Reco_QQ_mulpt[iqq]; if (thelptMu >= Reco_mu_glb_size) continue;
843     if (Reco_QQ_probChi2[iqq] > MIN_vtxprob_jpsi &&
844     //Reco_mu_glb_nhitstrack[thehptMu] > MIN_nhits_trk &&
845     Reco_mu_glb_normChi2[thehptMu] < MAX_normchi2_glb &&
846     //(((Reco_mu_glb_nhitsPixB[thehptMu] + Reco_mu_glb_nhitsPixE[thehptMu]) > MIN_nhits_pixel) || ((Reco_mu_glb_nhitsPixB[thehptMu] + Reco_mu_glb_nhitsPixE[thehptMu]) > MIN_nhits_pixel-1 && Reco_mu_glb_nhitsPix1Hit[thehptMu] == 1)) &&
847     fabs(Reco_mu_glb_d0[thehptMu]) < MAX_d0_trk &&
848     fabs(Reco_mu_glb_dz[thehptMu]) < MAX_dz_trk &&
849     //Reco_mu_glb_nhitstrack[thelptMu] > MIN_nhits_trk &&
850     Reco_mu_glb_normChi2[thelptMu] < MAX_normchi2_glb &&
851     //(((Reco_mu_glb_nhitsPixB[thelptMu] + Reco_mu_glb_nhitsPixE[thelptMu]) > MIN_nhits_pixel) || ((Reco_mu_glb_nhitsPixB[thelptMu] + Reco_mu_glb_nhitsPixE[thelptMu]) > MIN_nhits_pixel-1 && Reco_mu_glb_nhitsPix1Hit[thelptMu] == 1)) &&
852     fabs(Reco_mu_glb_d0[thelptMu]) < MAX_d0_trk &&
853     fabs(Reco_mu_glb_dz[thelptMu]) < MAX_dz_trk
854     ) {
855     return iqq;
856     }
857     }
858     }
859    
860     for (int iqq=0; iqq<Reco_QQ_size; iqq++) {
861    
862     if (Reco_QQ_sign[iqq] == 0 && Reco_QQ_type[iqq] == 1 ) {
863    
864     int thehptMu = Reco_QQ_muhpt[iqq]; if (thehptMu >= Reco_mu_glb_size) continue;
865     int thelptMu = Reco_QQ_mulpt[iqq]; if (thelptMu >= Reco_mu_trk_size) continue;
866    
867     if ( Reco_QQ_probChi2[iqq] > MIN_vtxprob_jpsi &&
868     //Reco_mu_glb_nhitstrack[thehptMu] > MIN_nhits_trk &&
869     Reco_mu_glb_normChi2[thehptMu] < MAX_normchi2_glb &&
870     //(((Reco_mu_glb_nhitsPixB[thehptMu] + Reco_mu_glb_nhitsPixE[thehptMu]) > MIN_nhits_pixel) || ((Reco_mu_glb_nhitsPixB[thehptMu] + Reco_mu_glb_nhitsPixE[thehptMu]) > MIN_nhits_pixel-1 && Reco_mu_glb_nhitsPix1Hit[thehptMu] == 1)) &&
871     fabs(Reco_mu_glb_d0[thehptMu]) < MAX_d0_trk &&
872     fabs(Reco_mu_glb_dz[thehptMu]) < MAX_dz_trk &&
873     Reco_mu_trk_nhitstrack[thelptMu] > MIN_nhits_trk &&
874     ((Reco_mu_trk_PIDmask[thelptMu] & (int)pow(2,5))/(int)pow(2,5) > 0 || (Reco_mu_trk_PIDmask[thelptMu] & (int)pow(2,8))/(int)pow(2,8) > 0) &&
875     //(Reco_mu_trk_nhitsPixB[thelptMu] + Reco_mu_trk_nhitsPixE[thelptMu]) > MIN_nhits_pixel &&
876     Reco_mu_trk_normChi2[thelptMu] < MAX_normchi2_trk &&
877     fabs(Reco_mu_trk_d0[thelptMu]) < MAX_d0_trk &&
878     fabs(Reco_mu_trk_dz[thelptMu]) < MAX_dz_trk) {
879    
880     TLorentzVector *theTrMumom = (TLorentzVector*)Reco_mu_trk_4mom->At(thelptMu);
881     if (theTrMumom->Perp() > thehighestPt) {
882     thehighestPt = theTrMumom->Perp();
883     theBest = iqq;
884     }
885     }
886     }
887     }
888    
889     if (theBest >= 0) return theBest;
890    
891     for (int iqq=0; iqq<Reco_QQ_size; iqq++) {
892    
893     if (Reco_QQ_sign[iqq] == 0 && Reco_QQ_type[iqq] == 2 ) {
894    
895     int thehptMu = Reco_QQ_muhpt[iqq]; if (thehptMu >= Reco_mu_trk_size) continue;
896     int thelptMu = Reco_QQ_mulpt[iqq]; if (thelptMu >= Reco_mu_trk_size) continue;
897    
898     if ( Reco_QQ_probChi2[iqq] > MIN_vtxprob_jpsi &&
899     Reco_mu_trk_nhitstrack[thehptMu] > MIN_nhits_trk &&
900     ((Reco_mu_trk_PIDmask[thehptMu] & (int)pow(2,5))/(int)pow(2,5) > 0 || (Reco_mu_trk_PIDmask[thehptMu] & (int)pow(2,8))/(int)pow(2,8) > 0) &&
901     //(Reco_mu_trk_nhitsPixB[thehptMu] + Reco_mu_trk_nhitsPixE[thehptMu]) > MIN_nhits_pixel &&
902     Reco_mu_trk_normChi2[thehptMu] < MAX_normchi2_trk &&
903     fabs(Reco_mu_trk_d0[thehptMu]) < MAX_d0_trk &&
904     fabs(Reco_mu_trk_dz[thehptMu]) < MAX_dz_trk &&
905     Reco_mu_trk_nhitstrack[thelptMu] > MIN_nhits_trk &&
906     ((Reco_mu_trk_PIDmask[thelptMu] & (int)pow(2,5))/(int)pow(2,5) > 0 || (Reco_mu_trk_PIDmask[thelptMu] & (int)pow(2,8))/(int)pow(2,8) > 0) &&
907     //(Reco_mu_trk_nhitsPixB[thelptMu] + Reco_mu_trk_nhitsPixE[thelptMu]) > MIN_nhits_pixel &&
908     Reco_mu_trk_normChi2[thelptMu] < MAX_normchi2_trk &&
909     fabs(Reco_mu_trk_d0[thelptMu]) < MAX_d0_trk &&
910     fabs(Reco_mu_trk_dz[thelptMu]) < MAX_dz_trk) {
911    
912     TLorentzVector *theTrMumom = (TLorentzVector*)Reco_mu_trk_4mom->At(thehptMu);
913     if (theTrMumom->Perp() > thehighestPt) {
914     thehighestPt = theTrMumom->Perp();
915     theBest = iqq;
916     }
917     }
918     }
919     }
920    
921     if (theBest >= 0) return theBest;
922    
923     for (int iqq=0; iqq<Reco_QQ_size; iqq++) {
924    
925     if (Reco_QQ_sign[iqq] == 0 && Reco_QQ_type[iqq] == 3 ) {
926    
927     int thehptMu = Reco_QQ_muhpt[iqq];
928     int thelptMu = Reco_QQ_mulpt[iqq];
929     if (thelptMu >= Reco_mu_cal_size) {
930     // cout << "Non deve succedere! cmIndex = " << thelptMu+1 << " cmSize = " << Reco_mu_cal_size << endl;
931     continue;
932     }
933    
934     if ( //Reco_mu_glb_nhitstrack[thehptMu] > MIN_nhits_trk &&
935     Reco_mu_glb_normChi2[thehptMu] < MAX_normchi2_glb &&
936     Reco_mu_cal_nhitstrack[thelptMu] > MIN_nhits_trk &&
937     Reco_mu_cal_normChi2[thelptMu] < MAX_normchi2_cal &&
938     Reco_mu_cal_caloComp[thelptMu] > MIN_caloComp) {
939    
940     TLorentzVector *theCaMumom = (TLorentzVector*)Reco_mu_cal_4mom->At(thelptMu);
941     if (theCaMumom->Perp() > thehighestPt) {
942     thehighestPt = theCaMumom->Perp();
943     theBest = iqq;
944     }
945     }
946     }
947     }
948    
949     return theBest;
950     }
951     //=========================================
952     Double_t ProjectQQ::deltaR(TLorentzVector* t, TLorentzVector* u){
953    
954     return sqrt(pow(t->Eta()-u->Eta(),2) +pow(PhiInRange(t->Phi()-u->Phi()),2));
955     }
956     //=========================================
957     double ProjectQQ::PhiInRange(double phi){
958    
959     double phiout = phi;
960    
961     if( phiout > 2*M_PI || phiout < -2*M_PI) {
962     phiout = fmod( phiout, 2*M_PI);
963     }
964     if (phiout <= -M_PI) phiout += 2*M_PI;
965     else if (phiout > M_PI) phiout -= 2*M_PI;
966    
967     return phiout;
968     }