]>
Commit | Line | Data |
---|---|---|
fda75e2d | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | // | |
16 | // Debug tree task | |
17 | // | |
18 | // Authors: | |
19 | // Markus Fasel <M.Fasel@gsi.de> | |
20 | // | |
21 | #include <TBits.h> | |
22 | #include <TString.h> | |
23 | ||
24 | #include "AliAnalysisManager.h" | |
25 | #include "AliCentrality.h" | |
26 | #include "AliESDEvent.h" | |
27 | #include "AliESDtrack.h" | |
28 | #include "AliESDtrackCuts.h" | |
29 | #include "AliESDVertex.h" | |
30 | #include "AliHFEcuts.h" | |
31 | #include "AliHFEsignalCuts.h" | |
32 | #include "AliInputEventHandler.h" | |
33 | #include "AliLog.h" | |
34 | #include "AliMCEvent.h" | |
35 | #include "AliMCEventHandler.h" | |
36 | #include "AliMCParticle.h" | |
d44e3c84 | 37 | #include "AliStack.h" |
fda75e2d | 38 | #include "AliPIDResponse.h" |
11ff28c5 | 39 | #include "AliTrackReference.h" |
fda75e2d | 40 | #include "AliVEvent.h" |
11ff28c5 | 41 | #include "AliHFEpidTPC.h" |
a8ef1999 | 42 | #include "AliHFEpidTRD.h" |
11ff28c5 | 43 | #include "AliHFEmcQA.h" |
fda75e2d | 44 | #include "TTreeStream.h" |
45 | ||
46 | #include "AliHFEdebugTreeTask.h" | |
47 | ||
48 | ClassImp(AliHFEdebugTreeTask) | |
49 | ||
50 | AliHFEdebugTreeTask::AliHFEdebugTreeTask(): | |
51 | AliAnalysisTaskSE(), | |
52 | fTrackCuts(NULL), | |
53 | fSignalCuts(NULL), | |
a8ef1999 | 54 | fTRDpid(NULL), |
11ff28c5 | 55 | fTPCpid(NULL), |
56 | fExtraCuts(NULL), | |
fda75e2d | 57 | fNclustersTPC(70), |
58 | fNclustersTPCPID(0), | |
59 | fNclustersITS(2), | |
60 | fFilename("HFEtree.root"), | |
11ff28c5 | 61 | fDebugTree(NULL), |
62 | fNparents(-1) | |
fda75e2d | 63 | { |
64 | ||
65 | } | |
66 | ||
67 | AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name): | |
68 | AliAnalysisTaskSE(name), | |
69 | fTrackCuts(NULL), | |
70 | fSignalCuts(NULL), | |
a8ef1999 | 71 | fTRDpid(NULL), |
11ff28c5 | 72 | fTPCpid(NULL), |
73 | fExtraCuts(NULL), | |
fda75e2d | 74 | fNclustersTPC(70), |
75 | fNclustersTPCPID(0), | |
76 | fNclustersITS(2), | |
77 | fFilename("HFEtree.root"), | |
11ff28c5 | 78 | fDebugTree(NULL), |
79 | fNparents(-1) | |
fda75e2d | 80 | { |
11ff28c5 | 81 | fTRDpid = new AliHFEpidTRD("QAtrdPID"); |
82 | fTPCpid = new AliHFEpidTPC("QAtpcPID"); | |
fda75e2d | 83 | } |
84 | ||
85 | AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){ | |
a8ef1999 | 86 | if(fDebugTree) delete fDebugTree; |
87 | if(fTRDpid) delete fTRDpid; | |
11ff28c5 | 88 | if(fTPCpid) delete fTPCpid; |
fda75e2d | 89 | } |
90 | ||
91 | void AliHFEdebugTreeTask::UserCreateOutputObjects(){ | |
92 | // | |
93 | // Create debug tree, signal cuts and track cuts | |
94 | // | |
95 | fDebugTree = new TTreeSRedirector(fFilename.Data()); | |
96 | ||
97 | fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition"); | |
98 | ||
99 | fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts"); | |
100 | fTrackCuts->CreateStandardCuts(); | |
101 | // Track cuts | |
102 | fTrackCuts->SetMinNClustersTPC(fNclustersTPC); | |
103 | fTrackCuts->SetMinRatioTPCclusters(0); | |
104 | fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable); | |
105 | fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID); | |
106 | fTrackCuts->SetMinNClustersITS(fNclustersITS); | |
107 | // Event cuts | |
108 | fTrackCuts->SetUseMixedVertex(kTRUE); | |
109 | fTrackCuts->SetVertexRange(10.); | |
110 | fTrackCuts->Initialize(); | |
a8ef1999 | 111 | |
11ff28c5 | 112 | fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts"); |
a8ef1999 | 113 | |
fda75e2d | 114 | } |
115 | ||
116 | void AliHFEdebugTreeTask::UserExec(Option_t *){ | |
117 | // | |
118 | // User Exec: Fill debug Tree | |
119 | // | |
120 | ||
121 | // Get PID response | |
122 | AliPIDResponse *pid = NULL; | |
123 | AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
124 | if(handler){ | |
125 | pid = handler->GetPIDResponse(); | |
126 | } else { | |
127 | AliError("No Handler"); | |
128 | } | |
129 | if(!pid){ | |
317ae21b | 130 | AliError("No PID response"); |
131 | return; | |
132 | } | |
133 | if(!fInputEvent) { | |
134 | AliError("No Input event"); | |
fda75e2d | 135 | return; |
136 | } | |
11ff28c5 | 137 | |
138 | AliESDtrack copyTrack; | |
fda75e2d | 139 | |
140 | fTrackCuts->SetRecEvent(fInputEvent); | |
141 | ||
142 | // Cut event | |
143 | if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){ | |
0e30407a | 144 | AliDebug(1, "Event flagged as pileup\n"); |
fda75e2d | 145 | return; |
146 | } | |
147 | if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){ | |
148 | AliDebug(1, "Event rejected by the event cuts\n"); | |
149 | return; | |
150 | } | |
151 | ||
152 | // Get run number | |
153 | Int_t run = fInputEvent->GetRunNumber(); | |
154 | ||
155 | // Derive trigger | |
156 | UInt_t trigger = fInputHandler->IsEventSelected(); | |
157 | Bool_t isMBTrigger = trigger & AliVEvent::kMB; | |
158 | Bool_t isCentralTrigger = trigger & AliVEvent::kCentral; | |
159 | Bool_t isSemicentralTrigger = trigger & AliVEvent::kSemiCentral; | |
160 | Bool_t isEMCALTrigger = trigger & AliVEvent::kEMCEJE; | |
161 | ||
162 | // Check if MC information is available | |
163 | Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL; | |
164 | if(mcthere){ | |
165 | fTrackCuts->SetMCEvent(fMCEvent); | |
166 | fSignalCuts->SetMCEvent(fMCEvent); | |
167 | } | |
168 | ||
d44e3c84 | 169 | // MC get stack |
170 | AliStack* stack = 0x0; | |
171 | if(mcthere){ | |
172 | stack = fMCEvent->Stack(); | |
173 | if(!stack) AliError("No Stack"); | |
174 | } | |
175 | ||
fda75e2d | 176 | // Get Primary Vertex |
177 | const AliVVertex *vertex = fInputEvent->GetPrimaryVertex(); | |
11ff28c5 | 178 | Double_t vtx[3]; |
179 | vertex->GetXYZ(vtx); | |
180 | Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors(); | |
fda75e2d | 181 | |
182 | // Get centrality | |
183 | Float_t centrality = -1.; | |
317ae21b | 184 | AliESDEvent *event = (dynamic_cast<AliESDEvent *>(fInputEvent)); |
185 | if(!event) return; | |
186 | TString beamtype = event->GetBeamType(); | |
187 | //printf("Beam type %s\n",(const char*)beamtype); | |
fda75e2d | 188 | if(!beamtype.CompareTo("Pb-Pb") || !beamtype.CompareTo("A-A")){ |
189 | // Heavy ion run | |
190 | AliDebug(1, "Heavy-Ion event\n"); | |
191 | AliCentrality *hicent = fInputEvent->GetCentrality(); | |
192 | centrality = hicent->GetCentralityPercentile("V0M"); | |
193 | } | |
194 | ||
11ff28c5 | 195 | if(!fExtraCuts){ |
196 | fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts"); | |
197 | } | |
198 | fExtraCuts->SetRecEventInfo(event); | |
199 | ||
200 | // Store event selection variables | |
201 | (*fDebugTree) << "EventDebug" | |
202 | << "Centrality=" << centrality | |
203 | << "VertexZ=" << vtx[2] | |
204 | << "NumberOfContributors=" << ncontrib | |
959ea9d8 | 205 | << "run=" << run |
11ff28c5 | 206 | << "\n"; |
207 | ||
fda75e2d | 208 | // Common variables |
209 | Double_t charge, eta, phi, momentum, transversemomentum; | |
210 | Int_t source; | |
211 | ||
212 | // Monte-Carlo loop | |
213 | if(mcthere){ | |
214 | AliMCParticle * mcpart; | |
215 | for(Int_t itrk = 0; itrk < fMCEvent->GetNumberOfTracks(); itrk++){ | |
216 | mcpart = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrk)); | |
217 | if(!mcpart) continue; | |
218 | if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mcpart)) continue; | |
219 | // Kinematics | |
220 | charge = mcpart->Charge() > 0. ? 1. : -1.; | |
221 | momentum = mcpart->P() * charge; | |
222 | transversemomentum = mcpart->Pt() * charge; | |
223 | eta = mcpart->Eta(); | |
224 | phi = mcpart->Phi(); | |
225 | Int_t pdg = mcpart->Particle()->GetPdgCode(); | |
226 | ||
227 | // Get Production Vertex in radial direction | |
228 | Double_t vx = mcpart->Particle()->Vx(), | |
0e30407a | 229 | vy = mcpart->Particle()->Vy(); |
fda75e2d | 230 | Double_t productionVertex = TMath::Sqrt(vx*vx+vy*vy); |
231 | ||
232 | // Get Mother PDG code of the particle | |
233 | Int_t motherPdg = 0; | |
234 | Int_t motherlabel = mcpart->Particle()->GetFirstMother(); | |
235 | if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){ | |
236 | AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel)); | |
237 | if(mother) motherPdg = mother->Particle()->GetPdgCode(); | |
238 | } | |
239 | ||
240 | // derive source | |
241 | source = 5; | |
242 | if(fSignalCuts->IsCharmElectron(mcpart)) source = 0; | |
243 | else if(fSignalCuts->IsBeautyElectron(mcpart)) source = 1; | |
244 | else if(fSignalCuts->IsGammaElectron(mcpart)) source = 2; | |
245 | else if(fSignalCuts->IsNonHFElectron(mcpart)) source = 3; | |
246 | else if(TMath::Abs(pdg) == 11) source = 4; | |
247 | else source = 5; | |
248 | ||
11ff28c5 | 249 | // Momemntum at the inner wall of the TPC |
250 | Double_t pTPC = 0., ptTPC = 0.; | |
251 | AliTrackReference *ref = FindTrackReference(mcpart, 80, 270, AliTrackReference::kTPC); | |
252 | if(ref){ | |
253 | pTPC = ref->P(); | |
254 | ptTPC = ref->Pt(); | |
255 | } | |
256 | ||
257 | ||
258 | ||
0e30407a | 259 | (*fDebugTree) << "MCDebug" |
fda75e2d | 260 | << "centrality=" << centrality |
261 | << "MBtrigger=" << isMBTrigger | |
262 | << "CentralTrigger=" << isCentralTrigger | |
263 | << "SemicentralTrigger=" << isSemicentralTrigger | |
264 | << "EMCALtrigger=" << isEMCALTrigger | |
265 | << "run=" << run | |
266 | << "p=" << momentum | |
267 | << "pt=" << transversemomentum | |
268 | << "eta=" << eta | |
269 | << "phi=" << phi | |
270 | << "pdg=" << pdg | |
271 | << "ProductionVertex=" << productionVertex | |
272 | << "motherPdg=" << motherPdg | |
cedf0381 | 273 | << "source=" << source |
fda75e2d | 274 | << "\n"; |
275 | } | |
276 | } | |
277 | ||
278 | AliESDtrack *track; | |
11ff28c5 | 279 | Double_t mcp, mcpt, mcptTPC, mcpTPC; // MC Variables added to the debug tree |
fda75e2d | 280 | for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){ |
281 | // fill the tree | |
282 | track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack)); | |
283 | if(!track) continue; | |
284 | // Cut track (Only basic track cuts) | |
285 | if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue; | |
286 | // Debug streaming of PID-related quantities | |
0e30407a | 287 | copyTrack.~AliESDtrack(); |
11ff28c5 | 288 | new(©Track) AliESDtrack(*track); |
289 | if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(©Track, AliHFEpidObject::kESDanalysis); // Apply Eta Correction on copy track | |
fda75e2d | 290 | Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron); |
11ff28c5 | 291 | Double_t nSigmaTPC = pid->NumberOfSigmasTPC(©Track, AliPID::kElectron); |
292 | //if(TMath::Abs(nSigmaTOF) > 5) continue; | |
fda75e2d | 293 | // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC |
cedf0381 | 294 | Double_t tPCdEdx = copyTrack.GetTPCsignal(); |
fda75e2d | 295 | // Signal, source and MCPID |
296 | Bool_t signal = kTRUE; | |
297 | source = 5; | |
11ff28c5 | 298 | mcp = mcpt = mcpTPC = mcptTPC = 0.; |
299 | ||
300 | ||
301 | ||
302 | Double_t bgcategory = 0.; | |
303 | Int_t mArr = -1; | |
304 | Int_t mesonID = -999; | |
305 | Double_t xr[3]={-999,-999,-999}; | |
306 | Double_t eR=-999; | |
307 | Double_t eZ=-999; | |
308 | Double_t unique=-999; | |
309 | Double_t mesonunique=-999; | |
310 | Double_t mesonR=-999; | |
311 | Double_t mesonZ=-999; | |
312 | Double_t mesonMomPdg=-999; | |
313 | Double_t mesonMomPt=-999; | |
314 | Double_t mesonGMomPdg=-999; | |
315 | Double_t mesonGGMomPdg=-999; | |
316 | Double_t mesonPt = -999; | |
317 | Double_t mceta = -999; | |
318 | Double_t mcphi = -999; | |
319 | Int_t mcpdg; | |
320 | ||
321 | if(mcthere){ | |
fda75e2d | 322 | // Signal |
323 | AliMCParticle *mctrack; | |
324 | if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){ | |
325 | if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; | |
326 | } | |
327 | ||
328 | // Source | |
329 | if(fSignalCuts->IsCharmElectron(track)) source = 0; | |
330 | else if(fSignalCuts->IsBeautyElectron(track)) source = 1; | |
331 | else if(fSignalCuts->IsGammaElectron(track)) source = 2; | |
332 | else if(fSignalCuts->IsNonHFElectron(track)) source = 3; | |
333 | else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4; | |
334 | else source = 5; | |
11ff28c5 | 335 | |
3024f297 | 336 | if(!mctrack) continue; |
337 | ||
11ff28c5 | 338 | // Kinematics |
339 | mcpt = mctrack->Pt(); | |
340 | mcp = mctrack->P(); | |
341 | mceta = mctrack->Eta(); | |
342 | mcphi = mctrack->Phi(); | |
343 | mcpdg = mctrack->Particle()->GetPdgCode(); | |
344 | ||
345 | AliTrackReference *ref = FindTrackReference(mctrack, 80, 270, AliTrackReference::kTPC); | |
346 | if(ref){ | |
347 | mcpTPC = ref->P(); | |
348 | mcptTPC = ref->Pt(); | |
349 | } | |
350 | ||
0e30407a | 351 | |
11ff28c5 | 352 | TParticle *mctrack1 = mctrack->Particle(); |
353 | mesonID=GetElecSourceMC(mctrack1); | |
354 | if(mesonID==AliHFEmcQA::kGammaPi0 || mesonID==AliHFEmcQA::kPi0) mArr=0; //pion | |
355 | else if(mesonID==AliHFEmcQA::kGammaEta || mesonID==AliHFEmcQA::kEta) mArr=1; //eta | |
356 | else if(mesonID==AliHFEmcQA::kGammaOmega || mesonID==AliHFEmcQA::kOmega) mArr=2; //omega | |
357 | else if(mesonID==AliHFEmcQA::kGammaPhi || mesonID==AliHFEmcQA::kPhi) mArr=3; //phi | |
358 | else if(mesonID==AliHFEmcQA::kGammaEtaPrime || mesonID==AliHFEmcQA::kEtaPrime) mArr=4; //etaprime | |
359 | else if(mesonID==AliHFEmcQA::kGammaRho0 || mesonID==AliHFEmcQA::kRho0) mArr=5; //rho | |
360 | ||
361 | mctrack->XvYvZv(xr); | |
0e30407a | 362 | |
11ff28c5 | 363 | eR= TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); |
364 | eZ = xr[2]; | |
365 | TParticle *mctrackt = mctrack->Particle(); | |
366 | unique=mctrackt->GetUniqueID(); | |
367 | ||
368 | AliMCParticle *mctrackmother = NULL; | |
cedf0381 | 369 | AliMCParticle *mctrackmother2 = NULL; |
11ff28c5 | 370 | |
371 | if(!(mArr<0)){ | |
0e30407a | 372 | if(mesonID>=AliHFEmcQA::kGammaPi0) { // conversion electron, be careful with the enum odering |
373 | Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label | |
374 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
375 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label | |
376 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
377 | mesonPt = mctrackmother->Pt(); //meson pt | |
378 | bgcategory = 1.; | |
379 | mctrackmother->XvYvZv(xr); | |
380 | mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); | |
381 | mesonZ = xr[2]; | |
382 | ||
383 | mctrackt = mctrackmother->Particle(); | |
384 | if(mctrackt){ | |
cedf0381 | 385 | mesonunique = mctrackt->GetUniqueID(); |
0e30407a | 386 | } |
cedf0381 | 387 | |
388 | Int_t glabel2=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother | |
389 | if((mctrackmother2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel2)))){ | |
390 | mesonMomPdg=mctrackmother2->PdgCode(); | |
391 | mesonMomPt=mctrackmother2->Pt(); | |
392 | } | |
393 | ||
0e30407a | 394 | if(glabel>fMCEvent->GetNumberOfPrimaries()) { |
395 | bgcategory = 2.; | |
396 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother | |
397 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
398 | mesonMomPdg=mctrackmother->PdgCode(); | |
399 | mesonMomPt=mctrackmother->Pt(); | |
400 | if(TMath::Abs(mctrackmother->PdgCode())==310){ | |
401 | bgcategory = 3.; | |
402 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother | |
403 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
404 | mesonGMomPdg=mctrackmother->PdgCode(); | |
405 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother | |
406 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
407 | mesonGGMomPdg=mctrackmother->PdgCode(); | |
408 | } | |
409 | } | |
410 | } | |
411 | } | |
412 | } | |
413 | } | |
414 | } | |
415 | } | |
416 | else{ // nonHFE except for the conversion electron | |
417 | Int_t glabel=TMath::Abs(mctrack->GetMother()); | |
418 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
419 | mesonPt = mctrackmother->Pt(); //meson pt | |
420 | bgcategory = -1.; | |
421 | mctrackmother->XvYvZv(xr); | |
422 | mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); | |
423 | mesonZ = xr[2]; | |
424 | ||
425 | mctrackt = mctrackmother->Particle(); | |
426 | if(mctrackt){ | |
427 | mesonunique = mctrackt->GetUniqueID(); | |
428 | } | |
429 | if(glabel>fMCEvent->GetNumberOfPrimaries()) { | |
430 | bgcategory = -2.; | |
431 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother | |
432 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
433 | mesonMomPdg=mctrackmother->PdgCode(); | |
434 | mesonMomPt=mctrackmother->Pt(); | |
435 | if(TMath::Abs(mctrackmother->PdgCode())==310){ | |
436 | bgcategory = -3.; | |
437 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother | |
438 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
11ff28c5 | 439 | mesonGMomPdg=mctrackmother->PdgCode(); |
0e30407a | 440 | glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother |
441 | if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){ | |
442 | mesonGGMomPdg=mctrackmother->PdgCode(); | |
443 | } | |
444 | } | |
445 | } | |
446 | } | |
447 | } | |
448 | } | |
449 | } | |
11ff28c5 | 450 | } |
451 | ||
452 | ||
453 | ||
fda75e2d | 454 | } |
455 | // Get V0 tag (if available) | |
456 | Int_t v0pid = -1; | |
457 | if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron; | |
458 | else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion; | |
459 | else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton; | |
460 | // Kinematics | |
461 | charge = track->Charge() > 0 ? 1. : -1.; | |
462 | eta = track->Eta(); | |
463 | phi = track->Phi(); | |
464 | momentum = track->P() * charge; | |
465 | transversemomentum = track->Pt() * charge; | |
11ff28c5 | 466 | Double_t momentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : 0.; |
467 | Double_t transversemomentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->Pt() : 0.; | |
fda75e2d | 468 | // ITS number of clusters |
469 | UChar_t nclustersITS = track->GetITSclusters(NULL); | |
470 | Double_t chi2matching = track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex)); | |
471 | Double_t chi2PerClusterITS = 0.0; | |
472 | if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS); | |
473 | // TPC number of clusters (different definitions) | |
76d0b522 | 474 | UChar_t nclustersTPCfit = track->GetTPCNcls(); |
475 | UChar_t nclustersTPCall = 0; | |
476 | const TBits &clusterTPC = track->GetTPCClusterMap(); | |
477 | nclustersTPCall = clusterTPC.CountBits(); | |
fda75e2d | 478 | UChar_t nclustersTPCPID = track->GetTPCsignalN(); |
479 | UChar_t nfindableTPC = track->GetTPCNclsF(); | |
76d0b522 | 480 | Double_t clusterRatioTPCfit = 0.0; |
481 | if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC); | |
482 | Double_t clusterRatioTPCall = 0.0; | |
483 | if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC); | |
fda75e2d | 484 | UChar_t nclustersTPCshared = 0; |
485 | Float_t ncrossedRowsTPC = track->GetTPCCrossedRows(); | |
486 | const TBits &sharedTPC = track->GetTPCSharedMap(); | |
487 | for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++; | |
488 | // TRD clusters and tracklets | |
489 | UChar_t nclustersTRD = track->GetTRDncls(); | |
a8ef1999 | 490 | UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID(); |
fda75e2d | 491 | // ITS and TRD acceptance maps |
91e50e2b | 492 | UChar_t hasClusterITS[6], statusITS[6], hasTrackletTRD[6]; |
fda75e2d | 493 | UChar_t itsPixel = track->GetITSClusterMap(); |
cedf0381 | 494 | for(Int_t icl = 0; icl < 6; icl++){ |
495 | hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0; | |
496 | if(CheckITSstatus(track, icl)) statusITS[icl] = 1; | |
497 | else statusITS[icl] = 0; | |
498 | } | |
a8ef1999 | 499 | Double_t trddEdxSum[6]; |
500 | for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;} | |
fda75e2d | 501 | for(Int_t itl = 0; itl < 6; itl++){ |
0e30407a | 502 | Int_t nSliceNonZero = 0; |
11ff28c5 | 503 | trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge |
fda75e2d | 504 | for(Int_t islice = 0; islice < 8; islice++){ |
505 | if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++; | |
506 | } | |
507 | hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0; | |
508 | } | |
509 | // TRD PID | |
510 | Double_t pidprobs[5]; | |
511 | track->GetTRDpid(pidprobs); | |
512 | Double_t likeEleTRD = pidprobs[0]; | |
513 | Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]); | |
a8ef1999 | 514 | Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6); |
515 | Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6); | |
516 | ||
fda75e2d | 517 | // DCA |
518 | Float_t b[2] = {0.,0.}; | |
519 | Float_t bCov[3] = {0.,0.,0.}; | |
520 | track->GetImpactParameters(b,bCov); | |
521 | Double_t dca = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); // impact parameter space | |
522 | Double_t dcaSR=0, dcaSZ=0; | |
523 | if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy | |
524 | if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z | |
525 | Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma | |
11ff28c5 | 526 | |
527 | // HFE DCA | |
528 | Double_t hfeb[2] = {-99.,-99.}; | |
529 | Double_t hfebCov[3] = {-999.,-999.,-999.}; | |
530 | fExtraCuts->GetHFEImpactParameters(track, hfeb, hfebCov); | |
531 | ||
0e30407a | 532 | Double_t tofdx= -999.0; |
533 | Double_t tofdz= -999.0; | |
534 | tofdx=track->GetTOFsignalDx(); | |
535 | tofdz=track->GetTOFsignalDz(); | |
536 | ||
d44e3c84 | 537 | // TOF track status |
538 | UInt_t status = 0; | |
539 | status = track->GetStatus(); | |
540 | Bool_t hasTOFout = status&AliESDtrack::kTOFout; | |
541 | Bool_t hasTOFtime = status&AliESDtrack::kTIME; | |
542 | Bool_t hasTOFpid = status&AliESDtrack::kTOFpid; | |
543 | Bool_t hasgoodTOF = kFALSE; | |
544 | if (hasTOFout && hasTOFtime && hasTOFpid) hasgoodTOF = kTRUE; | |
545 | ||
546 | // TRD track status | |
547 | Bool_t hasTRDin = status&AliESDtrack::kTRDin; | |
548 | ||
549 | ||
550 | // TOF mismatch (particle spectra group) | |
551 | Int_t mismatchlevel=0; // default value; in data always 0 | |
552 | if(mcthere){ | |
553 | Int_t tofLabel[3]; | |
554 | track->GetTOFLabel(tofLabel); | |
555 | if(TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) mismatchlevel=1; | |
556 | TParticle *matchedTrack = stack->Particle(TMath::Abs(tofLabel[0])); | |
557 | if(TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel())) | |
558 | { | |
559 | if(mismatchlevel==1) mismatchlevel=3; | |
560 | else mismatchlevel=2; | |
561 | } | |
562 | } | |
563 | ||
564 | ||
fda75e2d | 565 | // Fill Tree |
566 | (*fDebugTree) << "PIDdebug" | |
567 | << "centrality=" << centrality | |
568 | << "MBtrigger=" << isMBTrigger | |
569 | << "CentralTrigger=" << isCentralTrigger | |
570 | << "SemicentralTrigger=" << isSemicentralTrigger | |
571 | << "EMCALtrigger=" << isEMCALTrigger | |
572 | << "signal=" << signal | |
573 | << "source=" << source | |
574 | << "v0pid=" << v0pid | |
575 | << "run=" << run | |
576 | << "p=" << momentum | |
11ff28c5 | 577 | << "ptpc=" << momentumTPC |
fda75e2d | 578 | << "pt=" << transversemomentum |
11ff28c5 | 579 | << "pttpc=" << transversemomentumTPC |
580 | << "mcp=" << mcp | |
581 | << "mcpt=" << mcpt | |
582 | << "mcpTPC=" << mcpTPC | |
583 | << "mcptTPC=" << mcptTPC | |
584 | << "mceta=" << mceta | |
585 | << "mcphi=" << mcphi | |
586 | << "mcpdg=" << mcpdg | |
fda75e2d | 587 | << "eta=" << eta |
588 | << "phi=" << phi | |
a8ef1999 | 589 | << "ntracklets=" << ntrackletsTRDPID |
76d0b522 | 590 | << "nclustersTPC=" << nclustersTPCfit |
591 | << "nclustersTPCall=" << nclustersTPCall | |
fda75e2d | 592 | << "nclustersTPCPID=" << nclustersTPCPID |
593 | << "nclustersTPCshared=" << nclustersTPCshared | |
594 | << "ncrossedRowsTPC=" << ncrossedRowsTPC | |
76d0b522 | 595 | << "clusterRatioTPC=" << clusterRatioTPCfit |
596 | << "clusterRatioTPCall=" << clusterRatioTPCall | |
fda75e2d | 597 | << "nclustersITS=" << nclustersITS |
598 | << "nclusters=" << nclustersTRD | |
599 | << "chi2matching=" << chi2matching | |
cedf0381 | 600 | << "chi2PerClusterITS=" << chi2PerClusterITS |
fda75e2d | 601 | << "its0=" << hasClusterITS[0] |
602 | << "its1=" << hasClusterITS[1] | |
603 | << "its2=" << hasClusterITS[2] | |
604 | << "its3=" << hasClusterITS[3] | |
605 | << "its4=" << hasClusterITS[4] | |
606 | << "its5=" << hasClusterITS[5] | |
cedf0381 | 607 | << "statusITS0=" << statusITS[0] |
608 | << "statusITS1=" << statusITS[1] | |
609 | << "statusITS2=" << statusITS[2] | |
610 | << "statusITS3=" << statusITS[3] | |
611 | << "statusITS4=" << statusITS[4] | |
612 | << "statusITS5=" << statusITS[5] | |
fda75e2d | 613 | << "trd0=" << hasTrackletTRD[0] |
614 | << "trd1=" << hasTrackletTRD[1] | |
615 | << "trd2=" << hasTrackletTRD[2] | |
616 | << "trd3=" << hasTrackletTRD[3] | |
617 | << "trd4=" << hasTrackletTRD[4] | |
618 | << "trd5=" << hasTrackletTRD[5] | |
cedf0381 | 619 | << "TRDdEdxl0=" << trddEdxSum[0] |
620 | << "TRDdEdxl1=" << trddEdxSum[1] | |
621 | << "TRDdEdxl2=" << trddEdxSum[2] | |
622 | << "TRDdEdxl3=" << trddEdxSum[3] | |
623 | << "TRDdEdxl4=" << trddEdxSum[4] | |
624 | << "TRDdEdxl5=" << trddEdxSum[5] | |
fda75e2d | 625 | << "TOFsigmaEl=" << nSigmaTOF |
626 | << "TPCsigmaEl=" << nSigmaTPC | |
627 | << "TPCdEdx=" << tPCdEdx | |
628 | << "TRDlikeEl=" << likeEleTRD | |
629 | << "TRDlikeEln=" << likeEleTRDn | |
cedf0381 | 630 | << "trdtruncmean1=" << trdtruncmean1 |
a8ef1999 | 631 | << "trdtruncmean2=" << trdtruncmean2 |
fda75e2d | 632 | << "dcaR=" << b[0] |
633 | << "dcaZ=" << b[1] | |
634 | << "dca=" << dca | |
635 | << "dcaSR=" << dcaSR | |
636 | << "dcaSZ=" << dcaSZ | |
637 | << "dcaS=" << dcaS | |
11ff28c5 | 638 | << "hfedcaR=" << hfeb[0] |
639 | << "hfedcaZ=" << hfeb[1] | |
640 | << "hfedcacovR=" << hfebCov[0] | |
641 | << "hfedcacovZ=" << hfebCov[2] | |
fda75e2d | 642 | << "vx=" << vtx[0] |
643 | << "vy=" << vtx[1] | |
0e30407a | 644 | << "vz=" << vtx[2] |
cedf0381 | 645 | << "tofdx=" << tofdx |
d44e3c84 | 646 | << "tofdz=" << tofdz |
647 | << "statusTOFtracking=" << hasgoodTOF | |
648 | << "TOFmismatchlevel=" << mismatchlevel | |
649 | << "statusTRDtracking=" << hasTRDin | |
cedf0381 | 650 | << "ncontrib=" << ncontrib |
651 | << "mesonID=" << mesonID | |
652 | << "eR=" << eR | |
653 | << "mesonR=" << mesonR | |
654 | << "eZ=" << eZ | |
655 | << "mesonZ=" << mesonZ | |
656 | << "unique=" << unique | |
657 | << "mesonunique=" << mesonunique | |
658 | << "bgcategory=" << bgcategory | |
659 | << "mesonpt=" << mesonPt | |
11ff28c5 | 660 | << "mesonMomPdg=" << mesonMomPdg |
661 | << "mesonGMomPdg=" << mesonGMomPdg | |
662 | << "mesonGGMomPdg=" << mesonGGMomPdg | |
663 | << "mesonMomPt=" << mesonMomPt | |
fda75e2d | 664 | << "\n"; |
665 | } | |
666 | } | |
667 | ||
668 | ||
669 | void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; } | |
11ff28c5 | 670 | //___________________________________________________________ |
671 | AliTrackReference *AliHFEdebugTreeTask::FindTrackReference(AliMCParticle *track, Float_t minRadius, Float_t maxRadius, Int_t detectorID) | |
672 | { | |
673 | // | |
674 | // Find the track reference | |
675 | // | |
676 | AliTrackReference *ref = NULL, *reftmp; | |
677 | Float_t radius; | |
678 | for(Int_t iref = 0; iref < track->GetNumberOfTrackReferences(); iref++){ | |
679 | reftmp = track->GetTrackReference(iref); | |
680 | if(reftmp->DetectorId() != detectorID) continue; | |
681 | radius = reftmp->R(); | |
682 | if(radius >= minRadius && radius < maxRadius){ | |
683 | ref = reftmp; | |
684 | break; | |
685 | } | |
686 | if(radius > maxRadius) break; | |
687 | } | |
688 | return ref; | |
689 | } | |
690 | ||
691 | //__________________________________________ | |
692 | Int_t AliHFEdebugTreeTask::GetElecSourceMC(TParticle * const mcpart) | |
693 | { | |
694 | // decay particle's origin | |
695 | ||
696 | if(!mcpart){ | |
697 | AliDebug(1, "no mcparticle, return\n"); | |
698 | return -1; | |
699 | } | |
700 | ||
701 | if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return AliHFEmcQA::kMisID; | |
702 | ||
703 | Int_t origin = -1; | |
704 | Bool_t isFinalOpenCharm = kFALSE; | |
705 | ||
706 | Int_t iLabel = mcpart->GetFirstMother(); | |
707 | if (iLabel<0){ | |
708 | AliDebug(1, "Stack label is negative, return\n"); | |
709 | return -1; | |
710 | } | |
711 | ||
712 | AliMCParticle *mctrack = NULL; | |
713 | Int_t tmpMomLabel=0; | |
714 | if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; | |
715 | TParticle *partMother = mctrack->Particle(); | |
716 | TParticle *partMotherCopy = mctrack->Particle(); | |
717 | Int_t maPdgcode = partMother->GetPdgCode(); | |
718 | ||
0e30407a | 719 | // if the mother is charmed hadron |
720 | if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kCharm || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kCharm ) { | |
11ff28c5 | 721 | |
0e30407a | 722 | for (Int_t i=0; i<fNparents; i++){ |
11ff28c5 | 723 | if (abs(maPdgcode)==fParentSelect[0][i]){ |
724 | isFinalOpenCharm = kTRUE; | |
725 | } | |
0e30407a | 726 | } |
727 | if (!isFinalOpenCharm) return -1; | |
11ff28c5 | 728 | |
0e30407a | 729 | // iterate until you find B hadron as a mother or become top ancester |
730 | for (Int_t i=1; i<fgkMaxIter; i++){ | |
11ff28c5 | 731 | |
732 | Int_t jLabel = partMother->GetFirstMother(); | |
733 | if (jLabel == -1){ | |
734 | origin = AliHFEmcQA::kDirectCharm; | |
735 | return origin; | |
736 | } | |
737 | if (jLabel < 0){ // safety protection | |
738 | AliDebug(1, "Stack label is negative, return\n"); | |
739 | return -1; | |
740 | } | |
741 | ||
742 | // if there is an ancester | |
743 | if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; | |
744 | TParticle* grandMa = mctrack->Particle(); | |
745 | Int_t grandMaPDG = grandMa->GetPdgCode(); | |
746 | ||
747 | for (Int_t j=0; j<fNparents; j++){ | |
0e30407a | 748 | if (abs(grandMaPDG)==fParentSelect[1][j]){ |
749 | origin = AliHFEmcQA::kBeautyCharm; | |
750 | return origin; | |
751 | } | |
11ff28c5 | 752 | } |
753 | ||
754 | partMother = grandMa; | |
0e30407a | 755 | } // end of iteration |
756 | } // end of if | |
757 | else if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kBeauty || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kBeauty ) { | |
758 | for (Int_t i=0; i<fNparents; i++){ | |
11ff28c5 | 759 | if (abs(maPdgcode)==fParentSelect[1][i]){ |
760 | origin = AliHFEmcQA::kDirectBeauty; | |
761 | return origin; | |
762 | } | |
0e30407a | 763 | } |
764 | } // end of if | |
765 | else if ( abs(maPdgcode) == 22 ) { //conversion | |
766 | ||
767 | tmpMomLabel = partMotherCopy->GetFirstMother(); | |
768 | if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1; | |
769 | partMother = mctrack->Particle(); | |
770 | maPdgcode = partMother->GetPdgCode(); | |
771 | if ( abs(maPdgcode) == 111 ) { | |
772 | origin = AliHFEmcQA::kGammaPi0; | |
773 | return origin; | |
774 | } | |
775 | else if ( abs(maPdgcode) == 221 ) { | |
776 | origin = AliHFEmcQA::kGammaEta; | |
777 | return origin; | |
778 | } | |
779 | else if ( abs(maPdgcode) == 223 ) { | |
780 | origin = AliHFEmcQA::kGammaOmega; | |
781 | return origin; | |
782 | } | |
783 | else if ( abs(maPdgcode) == 333 ) { | |
784 | origin = AliHFEmcQA::kGammaPhi; | |
785 | return origin; | |
786 | } | |
787 | else if ( abs(maPdgcode) == 331 ) { | |
788 | origin = AliHFEmcQA::kGammaEtaPrime; | |
789 | return origin; | |
790 | } | |
791 | else if ( abs(maPdgcode) == 113 ) { | |
792 | origin = AliHFEmcQA::kGammaRho0; | |
793 | return origin; | |
794 | } | |
795 | else origin = AliHFEmcQA::kElse; | |
796 | //origin = kGamma; // finer category above | |
797 | return origin; | |
798 | ||
799 | } // end of if | |
800 | else if ( abs(maPdgcode) == 111 ) { | |
801 | origin = AliHFEmcQA::kPi0; | |
802 | return origin; | |
803 | } // end of if | |
804 | else if ( abs(maPdgcode) == 221 ) { | |
805 | origin = AliHFEmcQA::kEta; | |
806 | return origin; | |
807 | } // end of if | |
808 | else if ( abs(maPdgcode) == 223 ) { | |
809 | origin = AliHFEmcQA::kOmega; | |
810 | return origin; | |
811 | } // end of if | |
812 | else if ( abs(maPdgcode) == 333 ) { | |
813 | origin = AliHFEmcQA::kPhi; | |
814 | return origin; | |
815 | } // end of if | |
816 | else if ( abs(maPdgcode) == 331 ) { | |
817 | origin = AliHFEmcQA::kEtaPrime; | |
818 | return origin; | |
819 | } // end of if | |
820 | else if ( abs(maPdgcode) == 113 ) { | |
821 | origin = AliHFEmcQA::kRho0; | |
822 | return origin; | |
823 | } // end of if | |
824 | else{ | |
11ff28c5 | 825 | origin = AliHFEmcQA::kElse; |
0e30407a | 826 | } |
827 | return origin; | |
11ff28c5 | 828 | } |
cedf0381 | 829 | |
830 | //______________________________________________________ | |
831 | Bool_t AliHFEdebugTreeTask::CheckITSstatus( const AliESDtrack * const esdtrack, Int_t layer) const { | |
832 | // | |
833 | // Check whether ITS area is dead | |
834 | // | |
835 | Int_t itsStatus = 0; | |
836 | Int_t det; | |
837 | Float_t xloc, zloc; | |
838 | esdtrack->GetITSModuleIndexInfo(layer, det, itsStatus, xloc, zloc); | |
839 | Bool_t status; | |
840 | switch(itsStatus){ | |
841 | case 2: status = kFALSE; break; | |
842 | case 3: status = kFALSE; break; | |
843 | case 7: status = kFALSE; break; | |
844 | default: status = kTRUE; | |
845 | } | |
846 | return status; | |
847 | } | |
848 |