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