Update hfe package
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdebugTreeTask.cxx
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"
38 #include "AliTrackReference.h"
39 #include "AliVEvent.h"
40 #include "AliHFEpidTPC.h"
41 #include "AliHFEpidTRD.h"
42 #include "AliHFEmcQA.h"
43 #include "TTreeStream.h"
44
45 #include "AliHFEdebugTreeTask.h"
46
47 ClassImp(AliHFEdebugTreeTask)
48
49 AliHFEdebugTreeTask::AliHFEdebugTreeTask():
50   AliAnalysisTaskSE(),
51   fTrackCuts(NULL),
52   fSignalCuts(NULL),
53   fTRDpid(NULL),
54   fTPCpid(NULL),
55   fExtraCuts(NULL),
56   fNclustersTPC(70),
57   fNclustersTPCPID(0),
58   fNclustersITS(2),
59   fFilename("HFEtree.root"),
60   fDebugTree(NULL),
61   fNparents(-1)
62 {
63
64 }
65
66 AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
67   AliAnalysisTaskSE(name),
68   fTrackCuts(NULL),
69   fSignalCuts(NULL),
70   fTRDpid(NULL),
71   fTPCpid(NULL),
72   fExtraCuts(NULL),
73   fNclustersTPC(70),
74   fNclustersTPCPID(0),
75   fNclustersITS(2),
76   fFilename("HFEtree.root"),
77   fDebugTree(NULL),
78   fNparents(-1)
79 {
80   fTRDpid = new AliHFEpidTRD("QAtrdPID");
81   fTPCpid = new AliHFEpidTPC("QAtpcPID");
82 }
83
84 AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
85     if(fDebugTree) delete fDebugTree;
86     if(fTRDpid) delete fTRDpid;
87     if(fTPCpid) delete fTPCpid;
88 }
89
90 void 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();
110
111   fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
112
113 }
114
115 void 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){
129     AliError("No PID response");
130     return;
131   }
132   if(!fInputEvent) {
133     AliError("No Input event");
134     return;
135   }
136
137   AliESDtrack copyTrack;
138   
139   fTrackCuts->SetRecEvent(fInputEvent);
140
141   // Cut event
142   if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
143   AliDebug(1, "Event flagged as pileup\n");
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();
170   Double_t vtx[3];
171   vertex->GetXYZ(vtx);
172   Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
173
174   // Get centrality
175   Float_t centrality = -1.;
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);
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
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
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(), 
220               vy = mcpart->Particle()->Vy(); 
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       
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
250     (*fDebugTree) << "MCDebug"
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
264                     << "source="              << source
265                     << "\n";
266     }
267   }
268   
269   AliESDtrack *track;
270   Double_t mcp, mcpt, mcptTPC, mcpTPC;  // MC Variables added to the debug tree
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
278     copyTrack.~AliESDtrack();
279     new(&copyTrack) AliESDtrack(*track);
280     if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kESDanalysis); // Apply Eta Correction on copy track
281     Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
282     Double_t nSigmaTPC = pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron);
283     //if(TMath::Abs(nSigmaTOF) > 5) continue;
284     // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
285     Double_t tPCdEdx = copyTrack.GetTPCsignal();
286     // Signal, source and MCPID
287     Bool_t signal = kTRUE;
288     source = 5;
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){
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;
326
327       if(!mctrack) continue;
328
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
342     
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);
353     
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;
360       AliMCParticle *mctrackmother2 = NULL;
361
362       if(!(mArr<0)){
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){
376             mesonunique = mctrackt->GetUniqueID();
377           }
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
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)))){
430                                   mesonGMomPdg=mctrackmother->PdgCode();
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     }
441       }
442
443
444
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;
457     Double_t momentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : 0.;
458     Double_t transversemomentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->Pt() : 0.;
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)
465     UChar_t nclustersTPC = track->GetTPCncls();
466     UChar_t nclustersTPCPID = track->GetTPCsignalN();
467     UChar_t nfindableTPC =  track->GetTPCNclsF();
468     Double_t clusterRatioTPC = 0.0;
469     if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPC = static_cast<Double_t>(nclustersTPC)/static_cast<Double_t>(nfindableTPC);
470     UChar_t nclustersTPCshared = 0;
471     Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
472     const TBits &sharedTPC = track->GetTPCSharedMap();
473     for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
474     // TRD clusters and tracklets
475     UChar_t nclustersTRD = track->GetTRDncls();
476     UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
477     // ITS and TRD acceptance maps
478     UChar_t hasClusterITS[6], statusITS[5], hasTrackletTRD[6];
479     UChar_t itsPixel = track->GetITSClusterMap();
480     for(Int_t icl = 0; icl < 6; icl++){ 
481       hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
482       if(CheckITSstatus(track, icl)) statusITS[icl] = 1;
483       else statusITS[icl] = 0;
484     }
485     Double_t trddEdxSum[6];
486     for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
487     for(Int_t itl = 0; itl < 6; itl++){
488       Int_t nSliceNonZero = 0;
489       trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
490       for(Int_t islice = 0; islice < 8; islice++){
491         if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
492       }
493       hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0;
494     }
495     // TRD PID
496     Double_t pidprobs[5];
497     track->GetTRDpid(pidprobs);
498     Double_t likeEleTRD = pidprobs[0];
499     Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
500     Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
501     Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);
502
503     // DCA
504     Float_t b[2] = {0.,0.};
505     Float_t bCov[3] = {0.,0.,0.};
506     track->GetImpactParameters(b,bCov);
507     Double_t dca = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); // impact parameter space
508     Double_t dcaSR=0, dcaSZ=0;
509     if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
510     if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
511     Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
512
513     // HFE DCA
514     Double_t hfeb[2] = {-99.,-99.};
515     Double_t hfebCov[3] = {-999.,-999.,-999.};
516     fExtraCuts->GetHFEImpactParameters(track, hfeb, hfebCov);
517
518     Double_t tofdx= -999.0;
519     Double_t tofdz= -999.0;
520     tofdx=track->GetTOFsignalDx();
521     tofdz=track->GetTOFsignalDz();
522
523     // Fill Tree
524     (*fDebugTree) << "PIDdebug"
525                   << "centrality="          << centrality
526                   << "MBtrigger="           << isMBTrigger 
527                   << "CentralTrigger="      << isCentralTrigger
528                   << "SemicentralTrigger="  << isSemicentralTrigger
529                   << "EMCALtrigger="        << isEMCALTrigger
530                   << "signal="              << signal
531                   << "source="              << source
532                   << "v0pid="               << v0pid
533                   << "run="                 << run
534                   << "p="                   << momentum
535                   << "ptpc="                << momentumTPC
536                   << "pt="                  << transversemomentum
537                   << "pttpc="               << transversemomentumTPC
538                   << "mcp="                 << mcp
539                   << "mcpt="                << mcpt
540                   << "mcpTPC="              << mcpTPC
541                   << "mcptTPC="             << mcptTPC
542                   << "mceta="               << mceta
543                   << "mcphi="               << mcphi
544                   << "mcpdg="               << mcpdg
545                   << "eta="                 << eta
546                   << "phi="                 << phi
547                   << "ntracklets="          << ntrackletsTRDPID
548                   << "nclustersTPC="        << nclustersTPC
549                   << "nclustersTPCPID="     << nclustersTPCPID
550                   << "nclustersTPCshared="  << nclustersTPCshared
551                   << "ncrossedRowsTPC="     << ncrossedRowsTPC
552                   << "clusterRatioTPC="     << clusterRatioTPC
553                   << "nclustersITS="        << nclustersITS
554                   << "nclusters="           << nclustersTRD
555                   << "chi2matching="        << chi2matching
556                   << "chi2PerClusterITS="   << chi2PerClusterITS
557                   << "its0="                << hasClusterITS[0]
558                   << "its1="                << hasClusterITS[1]
559                   << "its2="                << hasClusterITS[2]
560                   << "its3="                << hasClusterITS[3]
561                   << "its4="                << hasClusterITS[4]
562                   << "its5="                << hasClusterITS[5]
563                   << "statusITS0="          << statusITS[0]
564                   << "statusITS1="          << statusITS[1]
565                   << "statusITS2="          << statusITS[2]
566                   << "statusITS3="          << statusITS[3]
567                   << "statusITS4="          << statusITS[4]
568                   << "statusITS5="          << statusITS[5]
569                   << "trd0="                << hasTrackletTRD[0]
570                   << "trd1="                << hasTrackletTRD[1]
571                   << "trd2="                << hasTrackletTRD[2]
572                   << "trd3="                << hasTrackletTRD[3]
573                   << "trd4="                << hasTrackletTRD[4]
574                   << "trd5="                << hasTrackletTRD[5]
575                   << "TRDdEdxl0="           << trddEdxSum[0]
576                   << "TRDdEdxl1="           << trddEdxSum[1]
577                   << "TRDdEdxl2="           << trddEdxSum[2]
578                   << "TRDdEdxl3="           << trddEdxSum[3]
579                   << "TRDdEdxl4="           << trddEdxSum[4]
580                   << "TRDdEdxl5="           << trddEdxSum[5]
581                   << "TOFsigmaEl="          << nSigmaTOF
582                   << "TPCsigmaEl="          << nSigmaTPC
583                   << "TPCdEdx="             << tPCdEdx
584                   << "TRDlikeEl="           << likeEleTRD
585                   << "TRDlikeEln="          << likeEleTRDn
586                   << "trdtruncmean1="       << trdtruncmean1
587                   << "trdtruncmean2="       << trdtruncmean2
588                   << "dcaR="                << b[0]
589                   << "dcaZ="                << b[1]
590                   << "dca="                 << dca
591                   << "dcaSR="               << dcaSR
592                   << "dcaSZ="               << dcaSZ
593                   << "dcaS="                << dcaS
594                   << "hfedcaR="             << hfeb[0]
595                   << "hfedcaZ="             << hfeb[1]
596                   << "hfedcacovR="          << hfebCov[0]
597                   << "hfedcacovZ="          << hfebCov[2]
598                   << "vx="                  << vtx[0]
599                   << "vy="                  << vtx[1]
600                   << "vz="                  << vtx[2]
601                   << "tofdx="                << tofdx
602                   << "tofdz="                << tofdz
603                   << "ncontrib="            << ncontrib
604                   << "mesonID="             << mesonID
605                   << "eR="                  << eR
606                   << "mesonR="              << mesonR
607                   << "eZ="                  << eZ
608                   << "mesonZ="              << mesonZ
609                   << "unique="              << unique
610                   << "mesonunique="         << mesonunique
611                   << "bgcategory="          << bgcategory
612                   << "mesonpt="             << mesonPt
613                   << "mesonMomPdg="         << mesonMomPdg
614                   << "mesonGMomPdg="         << mesonGMomPdg
615                   << "mesonGGMomPdg="         << mesonGGMomPdg
616                   << "mesonMomPt="         << mesonMomPt
617                   << "\n";
618   }
619 }
620
621
622 void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }
623 //___________________________________________________________
624 AliTrackReference *AliHFEdebugTreeTask::FindTrackReference(AliMCParticle *track, Float_t minRadius, Float_t maxRadius, Int_t detectorID)
625 {
626   //
627   // Find the track reference
628   //
629   AliTrackReference *ref = NULL, *reftmp;
630   Float_t radius;
631   for(Int_t iref = 0; iref < track->GetNumberOfTrackReferences(); iref++){
632     reftmp = track->GetTrackReference(iref);
633     if(reftmp->DetectorId() != detectorID) continue;
634     radius = reftmp->R();
635     if(radius >= minRadius && radius < maxRadius){
636       ref = reftmp;
637       break;
638     } 
639     if(radius > maxRadius) break;
640   }
641   return ref;
642 }
643
644 //__________________________________________
645 Int_t AliHFEdebugTreeTask::GetElecSourceMC(TParticle * const mcpart)
646 {
647   // decay particle's origin 
648
649   if(!mcpart){
650     AliDebug(1, "no mcparticle, return\n");
651     return -1;
652   }
653
654   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return AliHFEmcQA::kMisID;
655
656   Int_t origin = -1;
657   Bool_t isFinalOpenCharm = kFALSE;
658
659   Int_t iLabel = mcpart->GetFirstMother();
660   if (iLabel<0){
661     AliDebug(1, "Stack label is negative, return\n");
662     return -1;
663   }
664
665   AliMCParticle *mctrack = NULL;
666   Int_t tmpMomLabel=0;
667   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
668   TParticle *partMother = mctrack->Particle();
669   TParticle *partMotherCopy = mctrack->Particle();
670   Int_t maPdgcode = partMother->GetPdgCode();
671
672   // if the mother is charmed hadron  
673   if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kCharm || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kCharm ) {
674
675     for (Int_t i=0; i<fNparents; i++){
676         if (abs(maPdgcode)==fParentSelect[0][i]){
677           isFinalOpenCharm = kTRUE;
678         }
679     }
680     if (!isFinalOpenCharm) return -1;
681
682     // iterate until you find B hadron as a mother or become top ancester 
683     for (Int_t i=1; i<fgkMaxIter; i++){
684
685         Int_t jLabel = partMother->GetFirstMother();
686         if (jLabel == -1){
687           origin = AliHFEmcQA::kDirectCharm;
688           return origin;
689         }
690         if (jLabel < 0){ // safety protection
691           AliDebug(1, "Stack label is negative, return\n");
692           return -1;
693         }
694
695         // if there is an ancester
696         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
697         TParticle* grandMa = mctrack->Particle();
698         Int_t grandMaPDG = grandMa->GetPdgCode();
699
700         for (Int_t j=0; j<fNparents; j++){
701           if (abs(grandMaPDG)==fParentSelect[1][j]){
702             origin = AliHFEmcQA::kBeautyCharm;
703             return origin;
704           }
705         }
706
707         partMother = grandMa;
708     } // end of iteration 
709   } // end of if
710   else if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kBeauty || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kBeauty ) {
711     for (Int_t i=0; i<fNparents; i++){
712         if (abs(maPdgcode)==fParentSelect[1][i]){
713           origin = AliHFEmcQA::kDirectBeauty;
714           return origin;
715         }
716     }
717   } // end of if
718   else if ( abs(maPdgcode) == 22 ) { //conversion
719
720     tmpMomLabel = partMotherCopy->GetFirstMother();
721     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
722     partMother = mctrack->Particle();
723     maPdgcode = partMother->GetPdgCode();
724     if ( abs(maPdgcode) == 111 ) {
725       origin = AliHFEmcQA::kGammaPi0;
726       return origin;
727     } 
728     else if ( abs(maPdgcode) == 221 ) {
729       origin = AliHFEmcQA::kGammaEta;
730       return origin;
731     } 
732     else if ( abs(maPdgcode) == 223 ) {
733       origin = AliHFEmcQA::kGammaOmega;
734       return origin;
735     } 
736     else if ( abs(maPdgcode) == 333 ) {
737       origin = AliHFEmcQA::kGammaPhi;
738       return origin;
739     }
740     else if ( abs(maPdgcode) == 331 ) {
741       origin = AliHFEmcQA::kGammaEtaPrime;
742       return origin; 
743     }
744     else if ( abs(maPdgcode) == 113 ) {
745       origin = AliHFEmcQA::kGammaRho0;
746       return origin;
747     }
748     else origin = AliHFEmcQA::kElse;
749     //origin = kGamma; // finer category above
750     return origin;
751
752   } // end of if
753   else if ( abs(maPdgcode) == 111 ) {
754     origin = AliHFEmcQA::kPi0;
755     return origin;
756   } // end of if
757   else if ( abs(maPdgcode) == 221 ) {
758     origin = AliHFEmcQA::kEta;
759     return origin;
760   } // end of if
761   else if ( abs(maPdgcode) == 223 ) {
762     origin = AliHFEmcQA::kOmega;
763     return origin;
764   } // end of if
765   else if ( abs(maPdgcode) == 333 ) {
766     origin = AliHFEmcQA::kPhi;
767     return origin;
768   } // end of if
769   else if ( abs(maPdgcode) == 331 ) {
770     origin = AliHFEmcQA::kEtaPrime;
771     return origin;
772   } // end of if
773   else if ( abs(maPdgcode) == 113 ) {
774     origin = AliHFEmcQA::kRho0;
775     return origin;
776   } // end of if
777   else{ 
778     origin = AliHFEmcQA::kElse;
779   }
780   return origin;
781 }
782
783 //______________________________________________________
784 Bool_t AliHFEdebugTreeTask::CheckITSstatus( const AliESDtrack * const esdtrack, Int_t layer) const {
785   //
786   // Check whether ITS area is dead
787   //
788   Int_t itsStatus = 0;
789   Int_t det;
790   Float_t xloc, zloc;
791   esdtrack->GetITSModuleIndexInfo(layer, det, itsStatus, xloc, zloc);
792   Bool_t status;
793   switch(itsStatus){
794     case 2: status = kFALSE; break;
795     case 3: status = kFALSE; break;
796     case 7: status = kFALSE; break;
797     default: status = kTRUE;
798   }
799   return status;
800 }
801