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