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