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