]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaBtag.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Remove/comment out the code related to signal plus bacround mixing, since it is avail...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaBtag.cxx
CommitLineData
32301b07 1/**************************************************************************\r
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4* Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15/* $Id: $ */\r
16\r
17//_________________________________________________________________________\r
18//\r
ac53c372 19// Class for the electron identification and B-tagging.\r
32301b07 20// Clusters from EMCAL matched to tracks\r
21// and kept in the AOD. Few histograms produced.\r
22//\r
ac53c372 23// -- Author: T.R.P.Aronsson (Yale) J.L. Klay (Cal Poly), M. Heinz (Yale)\r
32301b07 24//////////////////////////////////////////////////////////////////////////////\r
25 \r
26// --- ROOT system --- \r
27#include <TH2F.h>\r
28#include <TH3F.h>\r
29#include <TParticle.h>\r
32301b07 30#include <TClonesArray.h>\r
31//#include <TObjString.h>\r
32//#include <Riostream.h>\r
33\r
34// --- Analysis system --- \r
35#include "AliAnaBtag.h" \r
36#include "AliCaloTrackReader.h"\r
37#include "AliMCAnalysisUtils.h"\r
38#include "AliAODCaloCluster.h"\r
39#include "AliFiducialCut.h"\r
40#include "AliAODTrack.h"\r
41#include "AliAODPid.h"\r
42#include "AliCaloPID.h"\r
43#include "AliAODMCParticle.h"\r
44#include "AliStack.h"\r
45#include "AliExternalTrackParam.h"\r
46#include "AliESDv0.h"\r
47#include "AliESDtrack.h"\r
48#include "AliAODJet.h"\r
49#include "AliAODEvent.h"\r
50#include "AliGenPythiaEventHeader.h"\r
51//#include "iostream.h"\r
52ClassImp(AliAnaBtag)\r
53 \r
54//____________________________________________________________________________\r
55AliAnaBtag::AliAnaBtag() \r
56: AliAnaPartCorrBaseClass(),fCalorimeter(""),\r
57 fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.),\r
58 fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),\r
59 fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),\r
ac53c372 60 fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),\r
32301b07 61 fhEmcalElectrons(0),fhTRDElectrons(0),fhTPCElectrons(0),fhDVM1(0),fhDVM2(0),fhJets(0),fhJetsAllEtaPhi(0),fhJetsLeadingBElectronEtaPhi(0),fhDVM1EtaPhi(0),fhBJetElectronDetaDphi(0),fhClusterEnergy(0),fhTestalle(0),fhResidual(0),fhElectrons(0),fhTracks(0)\r
62{\r
63 //default ctor\r
64 \r
65 //Initialize parameters\r
66 InitParameters();\r
67\r
32301b07 68}\r
69\r
32301b07 70//____________________________________________________________________________\r
71AliAnaBtag::~AliAnaBtag() \r
72{\r
73 //dtor\r
74\r
75}\r
76\r
77\r
78//________________________________________________________________________\r
79TList * AliAnaBtag::GetCreateOutputObjects()\r
80{ \r
81 // Create histograms to be saved in output file and \r
82 // store them in outputContainer\r
83 TList * outputContainer = new TList() ; \r
84 outputContainer->SetName("ElectronHistos") ; \r
85\r
86\r
ac53c372 87 fhEmcalElectrons = new TH1F("fhEmcalElectrons","",400,0,400);\r
32301b07 88 outputContainer->Add(fhEmcalElectrons);\r
89 \r
90 fhTRDElectrons = new TH1F("fhTRDElectrons","",400,0,400);\r
91 outputContainer->Add(fhTRDElectrons);\r
92\r
93 fhTPCElectrons = new TH1F("fhTPCElectrons","",400,0,400);\r
94 outputContainer->Add(fhTPCElectrons);\r
95\r
96 fhDVM1 = new TH1F("fhDVM1","",400,0,400);\r
97 outputContainer->Add(fhDVM1);\r
98\r
99 fhDVM2 = new TH1F("fhDVM2","",400,0,400);\r
100 outputContainer->Add(fhDVM2);\r
101\r
102\r
103 fhJets = new TH2F("fhJets","",400,0,400,20,0,20);\r
104 outputContainer->Add(fhJets);\r
105\r
32301b07 106 fhJetsAllEtaPhi = new TH2F("fhJetsAllEtaPhi","",100,-2,2,100,-2,8);\r
107 outputContainer->Add(fhJetsAllEtaPhi);\r
108\r
109 fhJetsLeadingBElectronEtaPhi = new TH2F("fhJetsLeadingBElectronEtaPhi","",100,-5,5,200,-10,10);\r
110 outputContainer->Add(fhJetsLeadingBElectronEtaPhi);\r
111\r
112 fhDVM1EtaPhi = new TH2F("fhDVM1EtaPhi","",100,-2,2,100,-2,8);\r
113 outputContainer->Add(fhDVM1EtaPhi);\r
114\r
115 fhBJetElectronDetaDphi = new TH2F("fhBJetElectronDetaDphi","",100,-5,5,200,-10,10);\r
116 outputContainer->Add(fhBJetElectronDetaDphi);\r
117\r
118 fhClusterEnergy = new TH1F("fhClusterEnergy","",100,0,10);\r
119 outputContainer->Add(fhClusterEnergy);\r
120\r
121 fhTestalle = new TH1F("fhTestalle","",400,0,400);\r
122 outputContainer->Add(fhTestalle);\r
123\r
124 fhResidual = new TH1F("fhResidual","",500,0,5);\r
125 outputContainer->Add(fhResidual);\r
126\r
127 fhElectrons = new TH2F("fhElectrons","",200,0,100,20,0,20);\r
128 outputContainer->Add(fhElectrons);\r
129\r
130 fhTracks = new TH2F("fhTracks","",200,0,100,20,0,20);\r
131 outputContainer->Add(fhTracks);\r
132\r
133\r
134 \r
135 return outputContainer ;\r
136 \r
137}\r
138\r
139//____________________________________________________________________________\r
140void AliAnaBtag::Init()\r
141{\r
142\r
143 //do some initialization\r
144 if(fCalorimeter == "PHOS") {\r
145 printf("AliAnaBtag::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");\r
146 fCalorimeter = "EMCAL";\r
147 }\r
148 if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){\r
149 printf("AliAnaBtag::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");\r
150 abort();\r
151 }\r
152\r
153}\r
154\r
155\r
156//____________________________________________________________________________\r
157void AliAnaBtag::InitParameters()\r
158{\r
159 \r
160 //Initialize the parameters of the analysis.\r
161 SetOutputAODClassName("AliAODPWG4Particle");\r
162 SetOutputAODName("PWG4Particle");\r
163\r
164 AddToHistogramsName("Btag_");\r
165\r
166 fCalorimeter = "EMCAL" ;\r
167 fpOverEmin = 0.5;\r
168 fpOverEmax = 1.2;\r
169 fResidualCut = 0.02;\r
170 fMinClusEne = 4.0;\r
171 //DVM B-tagging\r
172 fDrCut = 1.0; \r
173 fPairDcaCut = 0.02;\r
174 fDecayLenCut = 1.0;\r
175 fImpactCut = 0.5;\r
176 fAssocPtCut = 1.0;\r
177 fMassCut = 1.5;\r
178 fSdcaCut = 0.1;\r
179 fITSCut = 4;\r
180 //IPSig B-tagging\r
181 fNTagTrkCut = 2;\r
182 fIPSigCut = 3.0;\r
183 //Jet fiducial cuts\r
184 fJetEtaCut = 0.3;\r
185 fJetPhiMin = 1.8;\r
186 fJetPhiMax = 2.9;\r
187}\r
188\r
189//__________________________________________________________________\r
190void AliAnaBtag::MakeAnalysisFillAOD() \r
191{\r
192 //This reads in tracks, extrapolates to EMCAL, does p/E selectrons, identifies electron candidates\r
193 //After candidates are obtained, btagging and saving into AOD.\r
194\r
195 Double_t bfield = 0.;\r
196 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
197\r
198 //Select the calorimeter of the electron\r
199 if(fCalorimeter != "EMCAL") {\r
200 printf("This class not yet implemented for PHOS\n");\r
201 abort();\r
202 }\r
203 \r
204 TObjArray *cl = GetAODEMCAL();\r
205 \r
206 if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;\r
207 Int_t ntracks = GetAODCTS()->GetEntriesFast();\r
208 if(GetDebug() > 0)\r
209 printf("AliAnaBtag::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);\r
210\r
211 Int_t iCluster = -999;\r
212 Int_t ntot = cl->GetEntriesFast();\r
213\r
214 //CLUSTER STUFF\r
215 if(1){\r
216 for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
217 AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));\r
218 if(!clus) continue;\r
32301b07 219 fhClusterEnergy->Fill(clus->E());\r
220 }\r
221 }\r
222\r
223\r
224 for (Int_t itrk = 0; itrk < ntracks; itrk++) {////////////// track loop\r
225 iCluster = -999; //start with no match\r
226 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;\r
227\r
228 Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
229 Bool_t dcaOkay = GetDCA(track,imp,cov); //homegrown dca calculation until AOD is fixed\r
230 if(!dcaOkay) printf("AliAnaBtag::Problem computing DCA to primary vertex for track %d. Skipping it...\n",itrk);\r
231\r
232 fhTracks->Fill(track->Pt(),1);\r
233\r
234 AliAODPid* pid = (AliAODPid*) track->GetDetPid();\r
235 if(pid == 0) {\r
236 if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);\r
237 continue;\r
238 } \r
239 else {\r
240 Double_t emcpos[3];\r
241 pid->GetEMCALPosition(emcpos);\r
242 Double_t emcmom[3];\r
243 pid->GetEMCALMomentum(emcmom);\r
244 \r
245 TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);\r
246 TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);\r
247 Double_t tphi = pos.Phi();\r
248 Double_t teta = pos.Eta();\r
249 Double_t tmom = mom.Mag();\r
250\r
251 Bool_t in = kFALSE;\r
252 if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&\r
253 mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE; //Checks momentum?\r
254 //Also check the track\r
255 if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. &&\r
256 track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE;\r
257\r
258 if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD() - Track(Extrap) pt %2.2f(%2.2f), phi %2.2f(%2.2f), eta %2.2f(%2.2f) in fiducial cut %d\n",track->Pt(), mom.Pt(), track->Phi(), mom.Phi(), track->Eta(),mom.Eta(), in);\r
259\r
260 if(mom.Pt() > GetMinPt() && in) {\r
261 fhTracks->Fill(track->Pt(),3);\r
262 Double_t dEdx = pid->GetTPCsignal();\r
263 \r
264 //Int_t ntot = cl->GetEntriesFast();\r
265 Double_t res = 999.;\r
266 Double_t pOverE = -999.;\r
267 \r
268 Int_t pidProb = track->GetMostProbablePID();\r
269 Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;\r
270 Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;\r
271 Bool_t trkChgHad = kFALSE; if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) trkChgHad = kTRUE;\r
272\r
273\r
274 //Track Matching!\r
275 Bool_t emcEle = kFALSE; \r
276 double minRes=100.;\r
277 Double_t minR = 99;\r
278 Double_t minPe =-1;\r
279 Double_t minEp =-1;\r
280 Double_t minMult = -1;\r
281 Double_t minPt = -1;\r
282\r
283 for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
284 AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));\r
285 if(!clus) continue;\r
286\r
287 //As of 11-Oct-2009\r
288 //only select "good" clusters \r
289// if (clus->GetNCells() < 2 ) continue;\r
290// if (clus->GetNCells() > 30 ) continue;\r
291// if (clus->E() < fMinClusEne ) continue;\r
292// if (clus->GetDispersion() > 1 ) continue;\r
293// if (clus->GetM20() > 0.4 ) continue;\r
294// if (clus->GetM02() > 0.4 ) continue;\r
295// if (clus->GetM20() < 0.03 ) continue;\r
296// if (clus->GetM02() < 0.03 ) continue;\r
297// new optimized from ben. 2010May\r
298 if (clus->GetNCells() < 2 ) continue;\r
299 if (clus->GetNCells() > 35 ) continue;\r
300 if (clus->E() < 0 ) continue;\r
301 if (clus->GetDispersion() > 1.08 ) continue;\r
302 if (clus->GetM20() > 0.42 ) continue;\r
303 if (clus->GetM02() > 0.4 ) continue;\r
304 if (clus->GetM20() < 0 ) continue;\r
305 if (clus->GetM02() < 0.06 ) continue;\r
306 \r
307 \r
c8fe2783 308 Float_t x[3];\r
32301b07 309 clus->GetPosition(x);\r
310 TVector3 cluspos(x[0],x[1],x[2]);\r
311 Double_t deta = teta - cluspos.Eta();\r
312 Double_t dphi = tphi - cluspos.Phi();\r
313 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
314 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
315\r
316 res = sqrt(dphi*dphi + deta*deta);\r
317 \r
318 if(res<minRes){\r
319 minRes=res;\r
320 \r
321 }\r
322 \r
323 if(res < 0.0275) { // if(res < fResidualCut) { //Optimized from Ben\r
324 iCluster = iclus;\r
325 \r
326 Double_t energy = clus->E(); \r
327 if(energy > 0) pOverE = tmom/energy;\r
328 \r
329 if (res< minR) {\r
330 minR = res;\r
331 minPe = pOverE;\r
332 minEp = energy/tmom;\r
333 minMult = clus->GetNCells() ;\r
334 minPt = track->Pt();\r
335 }\r
336\r
32301b07 337\r
338 } else {\r
339 //unmatched\r
340 }//res cut\r
341 }//calo cluster loop\r
342\r
343 fhResidual->Fill(minRes);\r
344\r
345 if(minPe > 0.9 && minPe < 1.08) emcEle = kTRUE;// if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;\r
346 \r
347\r
348 \r
349 if(emcEle)\r
350 fhEmcalElectrons->Fill(track->Pt());\r
351 if(trkEle)\r
352 fhTRDElectrons->Fill(track->Pt());\r
353 if(tpcEle)\r
354 fhTPCElectrons->Fill(track->Pt());\r
355 \r
356 \r
357 //Take all emcal electrons, but the others only if pT < 10 GeV\r
358 if(emcEle && (tpcEle || trkEle) ) {\r
359 fhTestalle->Fill(track->Pt());\r
360 //B-tagging\r
361 if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");\r
362 Int_t dvmbtag = GetDVMBtag(track);\r
363 \r
364 \r
365 if(dvmbtag>0){\r
366 fhDVM1->Fill(track->Pt());\r
367 fhDVM1EtaPhi->Fill(track->Eta(),track->Phi()); \r
368 \r
369 }\r
370 if(dvmbtag>1)\r
371 fhDVM2->Fill(track->Pt());\r
372 \r
373\r
374 //Create particle to save in AOD\r
375 Double_t eMass = 0.511/1000; //mass in GeV\r
376 Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);\r
377 AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);\r
378 tr.SetLabel(track->GetLabel());\r
379 tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters\r
380 tr.SetTrackLabel(track->GetID(),-1); //sets the indices of the original tracks tr.SetTrackLabel(track->GetID(),-1) instead of itrk;\r
381\r
382 if(emcEle) {//PID determined by EMCAL\r
383 tr.SetDetector(fCalorimeter);\r
384 } else {\r
385 tr.SetDetector("CTS"); //PID determined by CTS\r
386 }\r
387\r
1e68a3f4 388 //if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);\r
32301b07 389 //Make this preserve sign of particle\r
390 if(track->Charge() < 0) tr.SetPdg(11); //electron is 11\r
391 else tr.SetPdg(-11); //positron is -11\r
392\r
393 tr.SetBtag(dvmbtag);\r
394 \r
395\r
396 //Check origin of the candiates\r
397 if(IsDataMC()){\r
398 tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));\r
399 if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());\r
400 }//Work with stack also \r
401 \r
402 AddAODParticle(tr);\r
403 \r
404 if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg()); \r
405 }//electron\r
406 }//pt, fiducial selection\r
407 }//pid check\r
408 }//track loop \r
409 \r
410\r
411\r
412\r
413 if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD() End fill AODs \n"); \r
414 \r
415}\r
416\r
417//__________________________________________________________________\r
418void AliAnaBtag::MakeAnalysisFillHistograms() \r
419{\r
420 //Do analysis and fill histograms\r
421\r
422 AliStack * stack = 0x0;\r
423// TParticle * primary = 0x0;\r
424\r
425\r
426\r
427\r
428 if(IsDataMC()) {\r
429 if(GetReader()->ReadStack()){\r
430 stack = GetMCStack() ; \r
431 if(!stack)\r
432 printf("AliAnaBtag::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
433 \r
434 }\r
435 }// is data and MC\r
436\r
437 ////////////////////////////////////\r
438 //Loop over jets and check for b-tag\r
439 ////////////////////////////////////\r
440 double maxjetEta=-4.;\r
441 double maxjetPhi=-4.;\r
442\r
443 Int_t njets = (GetReader()->GetOutputEvent())->GetNJets();\r
444 if(njets > 0) {\r
445 if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets. Performing b-jet tag analysis\n",njets);\r
446\r
447 for(Int_t ijet = 0; ijet < njets ; ijet++) {\r
448 AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;\r
449\r
450 if(ijet==0){\r
451 maxjetEta=jet->Eta();\r
452 maxjetPhi=jet->Phi();\r
453 }\r
454\r
455 fhJets->Fill(jet->Pt(),1);\r
456 fhJetsAllEtaPhi->Fill(jet->Eta(),jet->Phi());\r
457\r
458 if(jet->Pt() < 0.) continue; //This has to be adjusted depending on pp or AA!\r
459 fhJets->Fill(jet->Pt(),3); //All jets after pt cut\r
460\r
461 //Geometric EMCAL cut\r
462 if(TMath::Abs(jet->Eta()) > fJetEtaCut) continue;\r
463 if(jet->Phi() < fJetPhiMin || jet->Phi() > fJetPhiMax) continue;\r
464 fhJets->Fill(jet->Pt(),4); //All jets after geometric cut\r
465\r
466 Bool_t leadJet = kFALSE;\r
467 if (ijet==0){ \r
468 fhJets->Fill(jet->Pt(),5); //Leading jets\r
469 leadJet= kTRUE;\r
470 }\r
471\r
472\r
473 Bool_t dvmJet = kFALSE; \r
474 TRefArray* rt = jet->GetRefTracks();\r
475 Int_t ntrk = rt->GetEntries();\r
476\r
477 for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
478 AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
479 Bool_t isDVM = CheckIfBjet(jetTrack);\r
480 if(isDVM) dvmJet = kTRUE;\r
481 }\r
482\r
483 if(dvmJet)\r
484 fhJets->Fill(jet->Pt(),6);\r
485\r
486\r
487\r
488 if(IsDataMC()) {\r
489 //determine tagging efficiency & mis-tagging rate\r
490 //using b-quarks from stack\r
491 Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());\r
492 Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());\r
493 if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");\r
494 if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");\r
495 if (dvmJet && GetDebug() > 0) printf("== found DVM jet==\n");\r
496 }\r
497\r
498 } //jet loop\r
499 } //jets exist\r
500 \r
501\r
502\r
503\r
504 //Electron loop, read back electrons, fill histos\r
505 Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
506 if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);\r
507 \r
508 for(Int_t iaod = 0; iaod < naod ; iaod++){\r
509 AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
510 Int_t pdg = ele->GetPdg();\r
511\r
512\r
513 if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue; //not necessary..\r
514\r
515 //MC tag of this electron\r
516 // Int_t mctag = ele->GetTag();\r
517\r
518\r
519 fhElectrons->Fill(ele->Pt(),1); //All electrons\r
520 Bool_t photonic = kFALSE;\r
521 photonic = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
522 if(!photonic) fhElectrons->Fill(ele->Pt(),3); //nonphotonic electrons\r
523 if(photonic) fhElectrons->Fill(ele->Pt(),4); //photonic electrons\r
524\r
525 //Fill electron histograms \r
526 Float_t phiele = ele->Phi();\r
527 Float_t etaele = ele->Eta();\r
528\r
529\r
530 if(ele->GetBtag()>0){ // removed bit tag shit\r
531 fhElectrons->Fill(ele->Pt(),5);\r
532 if(!photonic) fhElectrons->Fill(ele->Pt(),6);\r
533 if(photonic) fhElectrons->Fill(ele->Pt(),7);\r
534 fhJetsLeadingBElectronEtaPhi->Fill(maxjetEta,maxjetPhi); \r
535 double deta=etaele-maxjetEta;\r
536 double dphi=phiele-maxjetPhi;\r
537 \r
538 //double r = sqrt((dphi*dphi)+(deta*deta));\r
539 fhBJetElectronDetaDphi->Fill(deta,dphi);\r
540 \r
541 }\r
542 \r
543 }//electron aod loop\r
544\r
545}\r
546\r
547//__________________________________________________________________\r
548Int_t AliAnaBtag::GetDVMBtag(AliAODTrack * tr )\r
549{\r
550 //This method uses the Displaced Vertex between electron-hadron\r
551 //pairs and the primary vertex to determine whether an electron is\r
552 //likely from a B hadron.\r
553\r
554 Int_t ncls1 = 0;\r
555 for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;\r
556\r
557\r
558 if (ncls1 < fITSCut) return 0;\r
559\r
560 Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
561 Bool_t dcaOkay = GetDCA(tr,imp,cov); //homegrown dca calculation until AOD is fixed \r
562 if(!dcaOkay) {\r
563 printf("AliAnaBtag::Problem computing DCA to primary vertex for track %d",tr->GetID());\r
564 return 0;\r
565 }\r
566\r
567 if (TMath::Abs(imp[0]) > fImpactCut ) return 0;\r
568 if (TMath::Abs(imp[1]) > fImpactCut ) return 0;\r
569\r
570 Int_t nvtx1 = 0;\r
571 Int_t nvtx2 = 0;\r
572 Int_t nvtx3 = 0;\r
573\r
574 for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {\r
575 //loop over assoc\r
576 AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
577 Int_t id1 = tr->GetID();\r
578 Int_t id2 = track2->GetID();\r
579 if(id1 == id2) continue;\r
580\r
581 Int_t ncls2 = 0;\r
582 for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;\r
583 if (ncls2 < fITSCut) continue;\r
584\r
585 if(tr->Pt()<6.&&track2->Pt() < 0.4) continue;\r
586 if(tr->Pt()>6.&&tr->Pt()<10.&&track2->Pt() < 0.2) continue;\r
587 if(tr->Pt()>10.&&track2->Pt() < 0.6) continue;\r
588\r
589\r
590 Double_t dphi = tr->Phi() - track2->Phi();\r
591 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
592 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
593 Double_t deta = tr->Eta() - track2->Eta();\r
594 Double_t dr = sqrt(deta*deta + dphi*dphi);\r
595\r
596 if(dr > fDrCut) continue;\r
597 \r
598 if(tr->Pt()<6.){\r
599 Double_t sDca = ComputeSignDca(tr, track2, 1.4,0.025);\r
600 if(sDca > 0.06) nvtx2++;\r
601 } \r
602 if(tr->Pt()>6.&&tr->Pt()<10.){\r
603 Double_t sDca = ComputeSignDca(tr, track2, 1.7,0.012);\r
604 if(sDca > 0.03) nvtx2++;\r
605 } \r
606 if(tr->Pt()>10.){\r
607 Double_t sDca = ComputeSignDca(tr, track2, 1.5,0.14);\r
608 if(sDca > 0.04) nvtx2++;\r
609 } \r
610\r
611\r
612 } //loop over hadrons\r
613\r
614 if(GetDebug() > 0) {\r
615 if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);\r
616 if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);\r
617 if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);\r
618 }\r
619\r
620\r
621\r
622\r
623 return nvtx2;\r
624\r
625}\r
626\r
627//__________________________________________________________________\r
628Double_t AliAnaBtag::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut, double pdcacut)\r
629{\r
630 //Compute the signed dca between two tracks\r
631 //and return the result\r
632\r
633 Double_t signDca=-999.;\r
634 if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);\r
635\r
636 //=====Now calculate DCA between both tracks======= \r
637 Double_t massE = 0.000511;\r
638 Double_t massK = 0.493677;\r
639\r
640 Double_t vertex[3] = {-999.,-999.,-999}; //vertex\r
641 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
642 GetReader()->GetVertex(vertex); //If only one file, get the vertex from there\r
643 //FIXME: Add a check for whether file 2 is PYTHIA or HIJING\r
644 //If PYTHIA, then set the vertex from file 2, if not, use the\r
645 //vertex from file 1\r
1e68a3f4 646 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);\r
32301b07 647 }\r
648 \r
649 TVector3 primV(vertex[0],vertex[1],vertex[2]) ;\r
650\r
651 if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;\r
652\r
653 AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);\r
654 AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);\r
655\r
656 Double_t bfield[3];\r
657 param1->GetBxByBz(bfield);\r
658 bfield[0]=0;\r
659 bfield[1]=0;\r
660 bfield[2]=5.;\r
661 Double_t bz = 5.; // = param1->GetBz();\r
662 //cout<<"In ComputeSignDCA, bfield[3] and bz: "<<bfield[0]<<" "<<bfield[1]<<" "<<bfield[2]<<" "<<bz<<endl;\r
663 Double_t xplane1 = 0.; Double_t xplane2 = 0.;\r
664 Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2);\r
665\r
666 param1->PropagateToBxByBz(xplane1,bfield);\r
667 param2->PropagateToBxByBz(xplane2,bfield);\r
668\r
669 Int_t id1 = 0, id2 = 0;\r
670 AliESDv0 bvertex(*param1,id1,*param2,id2);\r
671 Double_t vx,vy,vz;\r
672 bvertex.GetXYZ(vx,vy,vz);\r
673\r
674 Double_t emom[3];\r
675 Double_t hmom[3];\r
676 param1->PxPyPz(emom);\r
677 param2->PxPyPz(hmom);\r
678 TVector3 emomAtB(emom[0],emom[1],emom[2]);\r
679 TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);\r
680 TVector3 secvtxpt(vx,vy,vz);\r
681 TVector3 decayvector(0,0,0);\r
682 decayvector = secvtxpt - primV; //decay vector from PrimVtx\r
683 Double_t decaylength = decayvector.Mag();\r
684\r
685 //printf("\t JLK pairDCA = %2.2f\n",pairdca);\r
686\r
687 if(GetDebug() > 0) {\r
688 printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );\r
689 printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );\r
690 }\r
691\r
692\r
693\r
694 if (emomAtB.Mag()>0 && pairdca < pdcacut && decaylength < fDecayLenCut ) {\r
695 TVector3 sumMom = emomAtB+hmomAtB;\r
696 Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);\r
697 Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);\r
698 Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);\r
699 Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
700 Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));\r
701 Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();\r
702\r
703// if (masscut<1.1) fhDVMBtagQA2->Fill(sDca, mass);\r
704\r
705 if (mass > masscut && massPhot > 0.1) signDca = sDca;\r
706 \r
707 if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);\r
708 if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);\r
709 }\r
710\r
711 //clean up\r
712 delete param1;\r
713 delete param2;\r
714\r
715 return signDca;\r
716}\r
717\r
718\r
32301b07 719//__________________________________________________________________\r
720Bool_t AliAnaBtag::PhotonicV0(Int_t id) \r
721{\r
722 //This method checks to see whether a track that has been flagged as\r
723 //an electron was determined to match to a V0 candidate with\r
724 //invariant mass consistent with photon conversion\r
725\r
726 Bool_t itIS = kFALSE;\r
727\r
728 Double_t massEta = 0.547;\r
729 Double_t massRho0 = 0.770;\r
730 Double_t massOmega = 0.782;\r
731 Double_t massPhi = 1.020;\r
732 \r
733 //---Get V0s---\r
734 AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();\r
735 int nv0s = aod->GetNumberOfV0s();\r
736 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
737 AliAODv0 *v0 = aod->GetV0(iV0);\r
738 if (!v0) continue;\r
739 double radius = v0->RadiusV0();\r
740 double mass = v0->InvMass2Prongs(0,1,11,11);\r
741 if(GetDebug() > 0) {\r
742 printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );\r
743 printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);\r
744 printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );\r
745 }\r
746 if (mass < 0.100 ||\r
747 (mass > massEta-0.05 || mass < massEta+0.05) ||\r
748 (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
749 (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
750 (mass > massPhi-0.05 || mass < massPhi+0.05)) {\r
751 if ( id == v0->GetNegID() || id == v0->GetPosID()) {\r
752 itIS=kTRUE;\r
753 if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );\r
754 }\r
755 } }\r
756 return itIS;\r
757\r
758}\r
759\r
760//__________________________________________________________________\r
761Bool_t AliAnaBtag::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3]) \r
762{\r
763 //Use the Event vertex and AOD track information to get\r
764 //a real impact parameter for the track\r
765\r
766\r
767 Double_t maxD = 100000.; //max transverse IP\r
768 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
769 AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();\r
770 AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();\r
771 AliESDtrack esdTrack(track);\r
772 Double_t bfield[3];\r
773 esdTrack.GetBxByBz(bfield);\r
774 bfield[0]=0;\r
775 bfield[1]=0;\r
776 bfield[2]=5.;\r
777 \r
778 Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov);\r
779 return gotit;\r
780 }\r
781\r
782 return kFALSE;\r
783\r
784}\r
785//__________________________________________________________________\r
786Bool_t AliAnaBtag::CheckIfBjet(const AliAODTrack* track)\r
787{\r
788 Bool_t bjet = kFALSE;\r
789 Int_t trackId = track->GetID(); //get the index in the reader\r
790 Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
791 if(GetDebug() > 3) printf("AliAnaBtag::CheckIfBjet() - aod branch entries %d\n", naod);\r
792\r
793 for(Int_t iaod = 0; iaod < naod ; iaod++){\r
794 AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
795 Int_t electronlabel = ele->GetTrackLabel(0);\r
796 if(electronlabel != trackId) continue; //skip to the next one if they don't match\r
797 if(ele->GetBtag()>0)\r
798 bjet = kTRUE;\r
799 } \r
800\r
801 return bjet;\r
802} \r
803\r
804\r
805\r
806//__________________________________________________________________\r
807AliAODMCParticle* AliAnaBtag::GetMCParticle(Int_t ipart) \r
808{\r
809 //Get the MC particle at position ipart\r
810\r
811 AliAODMCParticle* aodprimary = 0x0;\r
812 TClonesArray * mcparticles0 = 0x0;\r
813 TClonesArray * mcparticles1 = 0x0;\r
814\r
815 if(GetReader()->ReadAODMCParticles()){\r
816 //Get the list of MC particles \r
817 mcparticles0 = GetReader()->GetAODMCParticles(0);\r
818 if(!mcparticles0 && GetDebug() > 0) {\r
819 printf("AliAnaBtag::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");\r
820 }\r
1e68a3f4 821// if(GetReader()->GetSecondInputAODTree()){\r
822// mcparticles1 = GetReader()->GetAODMCParticles(1);\r
823// if(!mcparticles1 && GetDebug() > 0) {\r
824// printf("AliAnaBtag::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");\r
825// }\r
826// }\r
32301b07 827\r
828 Int_t npart0 = mcparticles0->GetEntriesFast();\r
829 Int_t npart1 = 0;\r
830 if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
831 if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
832 else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);\r
833 if(!aodprimary) {\r
834 printf("AliAnaBtag::GetMCParticle() *** no primary ***: label %d \n", ipart);\r
835 return 0x0;\r
836 }\r
837\r
838 } else {\r
839 printf("AliAnaBtag::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
840 }\r
841 return aodprimary;\r
842\r
843}\r
844\r
845//__________________________________________________________________\r
846Bool_t AliAnaBtag::IsMcBJet(Double_t jeta, Double_t jphi)\r
847{\r
848 //Check the jet eta,phi against that of the b-quark\r
849 //to decide whether it is an MC B-jet\r
850 Bool_t bjet=kFALSE;\r
851\r
32301b07 852 AliStack* stack = 0x0;\r
853 \r
854 for(Int_t ipart = 0; ipart < 100; ipart++) {\r
855\r
856 Double_t pphi = -999.;\r
857 Double_t peta = -999.;\r
858 Int_t pdg = 0;\r
859 if(GetReader()->ReadStack()) {\r
860 stack = GetMCStack();\r
861 if(!stack) {\r
862 printf("AliAnaBtag::IsMCBJet() *** no stack ***: \n");\r
863 return kFALSE;\r
864 }\r
865 TParticle* primary = stack->Particle(ipart);\r
866 if (!primary) continue;\r
867 pdg = primary->GetPdgCode();\r
868 pphi = primary->Phi();\r
869 peta = primary->Eta();\r
870 } else if(GetReader()->ReadAODMCParticles()) {\r
871 AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
872 if(!aodprimary) continue;\r
873 pdg = aodprimary->GetPdgCode();\r
874 pphi = aodprimary->Phi();\r
875 peta = aodprimary->Eta();\r
876 }\r
877 if ( TMath::Abs(pdg) != 5) continue;\r
878 \r
879 // printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());\r
880 Double_t dphi = jphi - pphi;\r
881 Double_t deta = jeta - peta;\r
882 Double_t dr = sqrt(deta*deta + dphi*dphi);\r
883 \r
884 if (dr < 0.2) {\r
885 bjet=kTRUE;\r
886 //printf("MTH: **** found matching MC-Bjet: PDG=%d, pt=%f,dr=%f \n", pdgcode, primary->Pt(),dr );\r
887 break;\r
888 }\r
889 }\r
890 return bjet;\r
891\r
892}\r
893\r
894//__________________________________________________________________\r
895Bool_t AliAnaBtag::IsMcDJet(Double_t jeta, Double_t jphi)\r
896{\r
897 //Check if this jet is a charm jet\r
898 Bool_t cjet=kFALSE;\r
899\r
900 AliStack* stack = 0x0;\r
901\r
902 for(Int_t ipart = 0; ipart < 100; ipart++) {\r
903 \r
904 Double_t pphi = -999.;\r
905 Double_t peta = -999.;\r
906 Int_t pdg = 0;\r
907 if(GetReader()->ReadStack()) {\r
908 stack = GetMCStack();\r
909 if(!stack) {\r
910 printf("AliAnaBtag::IsMCDJet() *** no stack ***: \n");\r
911 return kFALSE;\r
912 }\r
913 TParticle* primary = stack->Particle(ipart);\r
914 if (!primary) continue;\r
915 pdg = primary->GetPdgCode();\r
916 pphi = primary->Phi();\r
917 peta = primary->Eta();\r
918 } else if(GetReader()->ReadAODMCParticles()) {\r
919 AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
920 if(!aodprimary) continue;\r
921 pdg = aodprimary->GetPdgCode();\r
922 pphi = aodprimary->Phi();\r
923 peta = aodprimary->Eta();\r
924 }\r
925\r
926 if ( TMath::Abs(pdg) != 4) continue;\r
927\r
928 Double_t dphi = jphi - pphi;\r
929 Double_t deta = jeta - peta;\r
930 Double_t dr = sqrt(deta*deta + dphi*dphi);\r
931 \r
932 if (dr < 0.2) {\r
933 cjet=kTRUE;\r
934 break;\r
935 }\r
936 }\r
937\r
938 return cjet;\r
939\r
940}\r
941\r
942//__________________________________________________________________\r
943void AliAnaBtag::Print(const Option_t * opt) const\r
944{\r
945 //Print some relevant parameters set for the analysis\r
946 \r
947 if(! opt)\r
948 return;\r
949 \r
950 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
951 AliAnaPartCorrBaseClass::Print(" ");\r
952\r
953 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;\r
954 printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax);\r
955 printf("residual cut = %f\n",fResidualCut);\r
956 printf("---DVM Btagging\n");\r
957 printf("max IP-cut (e,h) = %f\n",fImpactCut);\r
958 printf("min ITS-hits = %d\n",fITSCut);\r
959 printf("max dR (e,h) = %f\n",fDrCut);\r
960 printf("max pairDCA = %f\n",fPairDcaCut);\r
961 printf("max decaylength = %f\n",fDecayLenCut);\r
962 printf("min Associated Pt = %f\n",fAssocPtCut);\r
963 printf("---IPSig Btagging\n");\r
964 printf("min tag track = %d\n",fNTagTrkCut);\r
965 printf("min IP significance = %f\n",fIPSigCut);\r
966 printf(" \n") ;\r
967 \r
968} \r
969\r
970\r
971//__________________________________________________________________\r
972void AliAnaBtag::Terminate(TList* outputList)\r
973{\r
974 \r
975 //Do some plots to end\r
976 //Recover histograms from output histograms list, needed for\r
977 //distributed analysis. \r
978 //ReadHistograms(outputList);\r
979\r
980 printf(" AliAnaBtag::Terminate() *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;\r
981\r
982}\r
983\r