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