]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEdebugTreeTask.cxx
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdebugTreeTask.cxx
index 972d2e4bb507c9580c9fe66a861df778e4ab24d4..bc15d26c0585a64e6fd2c7db0fd82d1fc29af2ab 100644 (file)
 #include "AliMCEvent.h"
 #include "AliMCEventHandler.h"
 #include "AliMCParticle.h"
+#include "AliStack.h"
 #include "AliPIDResponse.h"
+#include "AliTrackReference.h"
 #include "AliVEvent.h"
+#include "AliHFEpidTPC.h"
+#include "AliHFEpidTRD.h"
+#include "AliHFEmcQA.h"
 #include "TTreeStream.h"
 
 #include "AliHFEdebugTreeTask.h"
@@ -46,11 +51,15 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask():
   AliAnalysisTaskSE(),
   fTrackCuts(NULL),
   fSignalCuts(NULL),
+  fTRDpid(NULL),
+  fTPCpid(NULL),
+  fExtraCuts(NULL),
   fNclustersTPC(70),
   fNclustersTPCPID(0),
   fNclustersITS(2),
   fFilename("HFEtree.root"),
-  fDebugTree(NULL)
+  fDebugTree(NULL),
+  fNparents(-1)
 {
 
 }
@@ -59,17 +68,24 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
   AliAnalysisTaskSE(name),
   fTrackCuts(NULL),
   fSignalCuts(NULL),
+  fTRDpid(NULL),
+  fTPCpid(NULL),
+  fExtraCuts(NULL),
   fNclustersTPC(70),
   fNclustersTPCPID(0),
   fNclustersITS(2),
   fFilename("HFEtree.root"),
-  fDebugTree(NULL)
+  fDebugTree(NULL),
+  fNparents(-1)
 {
-
+  fTRDpid = new AliHFEpidTRD("QAtrdPID");
+  fTPCpid = new AliHFEpidTPC("QAtpcPID");
 }
 
 AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
-  if(fDebugTree) delete fDebugTree;
+    if(fDebugTree) delete fDebugTree;
+    if(fTRDpid) delete fTRDpid;
+    if(fTPCpid) delete fTPCpid;
 }
 
 void AliHFEdebugTreeTask::UserCreateOutputObjects(){
@@ -92,6 +108,9 @@ void AliHFEdebugTreeTask::UserCreateOutputObjects(){
   fTrackCuts->SetUseMixedVertex(kTRUE);
   fTrackCuts->SetVertexRange(10.);
   fTrackCuts->Initialize();
+
+  fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+
 }
 
 void AliHFEdebugTreeTask::UserExec(Option_t *){
@@ -108,15 +127,21 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     AliError("No Handler");
   }
   if(!pid){
-    AliError("No PID response\n");
+    AliError("No PID response");
     return;
   }
+  if(!fInputEvent) {
+    AliError("No Input event");
+    return;
+  }
+
+  AliESDtrack copyTrack;
   
   fTrackCuts->SetRecEvent(fInputEvent);
 
   // Cut event
   if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
-   AliDebug(1, "Event flagged as pileup\n");
+  AliDebug(1, "Event flagged as pileup\n");
     return;
   }
   if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
@@ -141,12 +166,25 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     fSignalCuts->SetMCEvent(fMCEvent);
   }
 
+  // MC get stack
+  AliStack* stack = 0x0;
+  if(mcthere){
+      stack = fMCEvent->Stack();
+      if(!stack) AliError("No Stack");
+  }
+
   // Get Primary Vertex
   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
+  Double_t vtx[3];
+  vertex->GetXYZ(vtx);
+  Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
 
   // Get centrality
   Float_t centrality = -1.;
-  TString beamtype = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetESDRun()->GetBeamType();
+  AliESDEvent *event = (dynamic_cast<AliESDEvent *>(fInputEvent));
+  if(!event) return;
+  TString beamtype = event->GetBeamType();
+  //printf("Beam type %s\n",(const char*)beamtype);
   if(!beamtype.CompareTo("Pb-Pb") || !beamtype.CompareTo("A-A")){
     // Heavy ion run
     AliDebug(1, "Heavy-Ion event\n");
@@ -154,6 +192,19 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     centrality = hicent->GetCentralityPercentile("V0M");
   }
 
+  if(!fExtraCuts){
+    fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+  }
+  fExtraCuts->SetRecEventInfo(event);
+
+  // Store event selection variables
+  (*fDebugTree) << "EventDebug"
+                << "Centrality="              << centrality
+                << "VertexZ="                 << vtx[2]
+                << "NumberOfContributors="    << ncontrib
+                << "run="                     << run
+                << "\n";
+
   // Common variables
   Double_t charge, eta, phi, momentum, transversemomentum;
   Int_t source;
@@ -175,7 +226,7 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
       
       // Get Production Vertex in radial direction
       Double_t vx = mcpart->Particle()->Vx(), 
-               vy = mcpart->Particle()->Vy(); 
+              vy = mcpart->Particle()->Vy(); 
       Double_t productionVertex = TMath::Sqrt(vx*vx+vy*vy);
       
       // Get Mother PDG code of the particle
@@ -195,7 +246,17 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
       else if(TMath::Abs(pdg) == 11) source = 4;
       else source = 5;
       
-      (*fDebugTree) << "MCDebug" 
+      // Momemntum at the inner wall of the TPC
+      Double_t pTPC = 0., ptTPC = 0.;
+      AliTrackReference *ref = FindTrackReference(mcpart, 80, 270, AliTrackReference::kTPC);
+      if(ref){
+        pTPC = ref->P();
+        ptTPC = ref->Pt();
+      }
+
+
+
+    (*fDebugTree) << "MCDebug"
                     << "centrality="          << centrality
                     << "MBtrigger="           << isMBTrigger 
                     << "CentralTrigger="      << isCentralTrigger
@@ -215,6 +276,7 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
   }
   
   AliESDtrack *track;
+  Double_t mcp, mcpt, mcptTPC, mcpTPC;  // MC Variables added to the debug tree
   for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
     // fill the tree
     track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
@@ -222,15 +284,41 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     // Cut track (Only basic track cuts)
     if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     // Debug streaming of PID-related quantities
+    copyTrack.~AliESDtrack();
+    new(&copyTrack) AliESDtrack(*track);
+    if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kESDanalysis); // Apply Eta Correction on copy track
     Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
-    Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
-    if(TMath::Abs(nSigmaTOF) > 5) continue;
+    Double_t nSigmaTPC = pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron);
+    //if(TMath::Abs(nSigmaTOF) > 5) continue;
     // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
-    Double_t tPCdEdx = track->GetTPCsignal();
+    Double_t tPCdEdx = copyTrack.GetTPCsignal();
     // Signal, source and MCPID
     Bool_t signal = kTRUE;
     source = 5;
-    if(mcthere){      
+    mcp = mcpt = mcpTPC = mcptTPC = 0.;
+
+
+
+    Double_t bgcategory = 0.;
+    Int_t mArr = -1;
+    Int_t mesonID = -999;
+    Double_t xr[3]={-999,-999,-999};
+    Double_t eR=-999;
+    Double_t eZ=-999;
+    Double_t unique=-999;
+    Double_t mesonunique=-999;
+    Double_t mesonR=-999;
+    Double_t mesonZ=-999;
+    Double_t mesonMomPdg=-999;
+    Double_t mesonMomPt=-999;
+    Double_t mesonGMomPdg=-999;
+    Double_t mesonGGMomPdg=-999;
+    Double_t mesonPt = -999;
+    Double_t mceta = -999;
+    Double_t mcphi = -999;
+    Int_t mcpdg;
+
+    if(mcthere){
       // Signal
       AliMCParticle *mctrack;
       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){
@@ -244,6 +332,125 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
       else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
       else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4;
       else source = 5;
+
+      if(!mctrack) continue;
+
+      // Kinematics
+      mcpt  = mctrack->Pt();
+      mcp   = mctrack->P();
+      mceta = mctrack->Eta();
+      mcphi = mctrack->Phi();
+      mcpdg = mctrack->Particle()->GetPdgCode();
+
+      AliTrackReference *ref = FindTrackReference(mctrack, 80, 270, AliTrackReference::kTPC);
+      if(ref){
+        mcpTPC = ref->P();
+        mcptTPC = ref->Pt();
+      }
+
+    
+      TParticle *mctrack1 = mctrack->Particle();
+      mesonID=GetElecSourceMC(mctrack1);
+      if(mesonID==AliHFEmcQA::kGammaPi0 || mesonID==AliHFEmcQA::kPi0) mArr=0;                //pion
+      else if(mesonID==AliHFEmcQA::kGammaEta || mesonID==AliHFEmcQA::kEta) mArr=1;           //eta
+      else if(mesonID==AliHFEmcQA::kGammaOmega || mesonID==AliHFEmcQA::kOmega) mArr=2;       //omega
+      else if(mesonID==AliHFEmcQA::kGammaPhi || mesonID==AliHFEmcQA::kPhi) mArr=3;           //phi
+      else if(mesonID==AliHFEmcQA::kGammaEtaPrime || mesonID==AliHFEmcQA::kEtaPrime) mArr=4; //etaprime
+      else if(mesonID==AliHFEmcQA::kGammaRho0 || mesonID==AliHFEmcQA::kRho0) mArr=5;         //rho
+    
+      mctrack->XvYvZv(xr);
+    
+      eR= TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+      eZ = xr[2];
+      TParticle *mctrackt = mctrack->Particle();
+      unique=mctrackt->GetUniqueID();
+    
+      AliMCParticle *mctrackmother = NULL;
+      AliMCParticle *mctrackmother2 = NULL;
+
+      if(!(mArr<0)){
+    if(mesonID>=AliHFEmcQA::kGammaPi0) {  // conversion electron, be careful with the enum odering
+        Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
+        if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+      glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
+      if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+          mesonPt = mctrackmother->Pt(); //meson pt
+          bgcategory = 1.;
+          mctrackmother->XvYvZv(xr);
+          mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+          mesonZ = xr[2];
+
+          mctrackt = mctrackmother->Particle();
+          if(mctrackt){
+            mesonunique = mctrackt->GetUniqueID();
+          }
+
+        Int_t glabel2=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+        if((mctrackmother2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel2)))){
+            mesonMomPdg=mctrackmother2->PdgCode();
+            mesonMomPt=mctrackmother2->Pt();
+        }
+
+          if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+        bgcategory = 2.;
+        glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+        if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+            mesonMomPdg=mctrackmother->PdgCode();
+            mesonMomPt=mctrackmother->Pt();
+            if(TMath::Abs(mctrackmother->PdgCode())==310){
+          bgcategory = 3.;
+          glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+          if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+              mesonGMomPdg=mctrackmother->PdgCode();
+              glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+              if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+            mesonGGMomPdg=mctrackmother->PdgCode();
+              }
+          }
+            }
+        }
+          }
+      }
+        }
+    }
+    else{ // nonHFE except for the conversion electron
+        Int_t glabel=TMath::Abs(mctrack->GetMother());
+        if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+      mesonPt = mctrackmother->Pt(); //meson pt
+      bgcategory = -1.;
+      mctrackmother->XvYvZv(xr);
+      mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+      mesonZ = xr[2];
+
+      mctrackt = mctrackmother->Particle();
+      if(mctrackt){
+          mesonunique = mctrackt->GetUniqueID();
+      }
+      if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+          bgcategory = -2.;
+          glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+          if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+        mesonMomPdg=mctrackmother->PdgCode();
+        mesonMomPt=mctrackmother->Pt();
+        if(TMath::Abs(mctrackmother->PdgCode())==310){
+            bgcategory = -3.;
+            glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+            if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                                  mesonGMomPdg=mctrackmother->PdgCode();
+          glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+          if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                                    mesonGGMomPdg=mctrackmother->PdgCode();
+          }
+            }
+        }
+          }
+      }
+        }
+    }
+      }
+
+
+
     }
     // Get V0 tag (if available)
     Int_t v0pid = -1;
@@ -256,30 +463,44 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     phi = track->Phi();
     momentum = track->P() * charge;
     transversemomentum = track->Pt() * charge;
+    Double_t momentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : 0.;
+    Double_t transversemomentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->Pt() : 0.;
     // ITS number of clusters
     UChar_t nclustersITS = track->GetITSclusters(NULL);
     Double_t chi2matching =  track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex));
     Double_t chi2PerClusterITS = 0.0;
     if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS);
     // TPC number of clusters (different definitions)
-    UChar_t nclustersTPC = track->GetTPCncls();
+    UChar_t nclustersTPCfit = track->GetTPCNcls();
+    UChar_t nclustersTPCall = 0;
+    const TBits &clusterTPC = track->GetTPCClusterMap();
+    nclustersTPCall = clusterTPC.CountBits();
     UChar_t nclustersTPCPID = track->GetTPCsignalN();
     UChar_t nfindableTPC =  track->GetTPCNclsF();
-    Double_t clusterRatioTPC = 0.0;
-    if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPC = static_cast<Double_t>(nclustersTPC)/static_cast<Double_t>(nfindableTPC);
+    Double_t clusterRatioTPCfit = 0.0;
+    if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC);
+    Double_t clusterRatioTPCall = 0.0;
+    if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC);
     UChar_t nclustersTPCshared = 0;
     Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
     const TBits &sharedTPC = track->GetTPCSharedMap();
     for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
     // TRD clusters and tracklets
     UChar_t nclustersTRD = track->GetTRDncls();
-    UChar_t ntracklets = track->GetTRDntrackletsPID();
+    UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
     // ITS and TRD acceptance maps
-    UChar_t hasClusterITS[6], hasTrackletTRD[6];
+    UChar_t hasClusterITS[6], statusITS[6], hasTrackletTRD[6];
     UChar_t itsPixel = track->GetITSClusterMap();
-    for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
+    for(Int_t icl = 0; icl < 6; icl++){ 
+      hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
+      if(CheckITSstatus(track, icl)) statusITS[icl] = 1;
+      else statusITS[icl] = 0;
+    }
+    Double_t trddEdxSum[6];
+    for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
     for(Int_t itl = 0; itl < 6; itl++){
       Int_t nSliceNonZero = 0;
+      trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
       for(Int_t islice = 0; islice < 8; islice++){
         if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
       }
@@ -290,6 +511,9 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     track->GetTRDpid(pidprobs);
     Double_t likeEleTRD = pidprobs[0];
     Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
+    Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
+    Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);
+
     // DCA
     Float_t b[2] = {0.,0.};
     Float_t bCov[3] = {0.,0.,0.};
@@ -299,10 +523,45 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
     if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
     Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
-    // Vertex
-    Double_t vtx[3];
-    vertex->GetXYZ(vtx);
-    Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
+
+    // HFE DCA
+    Double_t hfeb[2] = {-99.,-99.};
+    Double_t hfebCov[3] = {-999.,-999.,-999.};
+    fExtraCuts->GetHFEImpactParameters(track, hfeb, hfebCov);
+
+    Double_t tofdx= -999.0;
+    Double_t tofdz= -999.0;
+    tofdx=track->GetTOFsignalDx();
+    tofdz=track->GetTOFsignalDz();
+
+    // TOF track status
+    UInt_t status = 0;
+    status = track->GetStatus();
+    Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
+    Bool_t hasTOFtime = status&AliESDtrack::kTIME;
+    Bool_t hasTOFpid  = status&AliESDtrack::kTOFpid;
+    Bool_t hasgoodTOF     = kFALSE;
+    if (hasTOFout && hasTOFtime && hasTOFpid) hasgoodTOF = kTRUE;
+
+    // TRD track status
+    Bool_t hasTRDin = status&AliESDtrack::kTRDin;
+
+
+    // TOF mismatch (particle spectra group)
+    Int_t mismatchlevel=0; // default value; in data always 0
+    if(mcthere){
+       Int_t tofLabel[3];
+       track->GetTOFLabel(tofLabel);
+        if(TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) mismatchlevel=1;
+       TParticle *matchedTrack = stack->Particle(TMath::Abs(tofLabel[0]));
+       if(TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
+       {
+           if(mismatchlevel==1) mismatchlevel=3;
+            else mismatchlevel=2;
+       }
+    }
+
+
     // Fill Tree
     (*fDebugTree) << "PIDdebug"
                   << "centrality="          << centrality
@@ -315,49 +574,275 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
                   << "v0pid="               << v0pid
                   << "run="                 << run
                   << "p="                   << momentum
+                  << "ptpc="                << momentumTPC
                   << "pt="                  << transversemomentum
+                  << "pttpc="               << transversemomentumTPC
+                  << "mcp="                 << mcp
+                  << "mcpt="                << mcpt
+                  << "mcpTPC="              << mcpTPC
+                  << "mcptTPC="             << mcptTPC
+                  << "mceta="               << mceta
+                  << "mcphi="               << mcphi
+                  << "mcpdg="               << mcpdg
                   << "eta="                 << eta
                   << "phi="                 << phi
-                  << "ntracklets="          << ntracklets
-                  << "nclustersTPC="        << nclustersTPC
+                  << "ntracklets="          << ntrackletsTRDPID
+                  << "nclustersTPC="        << nclustersTPCfit
+                 << "nclustersTPCall="     << nclustersTPCall
                   << "nclustersTPCPID="     << nclustersTPCPID
                   << "nclustersTPCshared="  << nclustersTPCshared
                   << "ncrossedRowsTPC="     << ncrossedRowsTPC
-                  << "clusterRatioTPC="     << clusterRatioTPC
+                  << "clusterRatioTPC="     << clusterRatioTPCfit
+                  << "clusterRatioTPCall="  << clusterRatioTPCall
                   << "nclustersITS="        << nclustersITS
                   << "nclusters="           << nclustersTRD
                   << "chi2matching="        << chi2matching
-                 << "chi2PerClusterITS="   << chi2PerClusterITS
+                  << "chi2PerClusterITS="   << chi2PerClusterITS
                   << "its0="                << hasClusterITS[0]
                   << "its1="                << hasClusterITS[1]
                   << "its2="                << hasClusterITS[2]
                   << "its3="                << hasClusterITS[3]
                   << "its4="                << hasClusterITS[4]
                   << "its5="                << hasClusterITS[5]
+                  << "statusITS0="          << statusITS[0]
+                  << "statusITS1="          << statusITS[1]
+                  << "statusITS2="          << statusITS[2]
+                  << "statusITS3="          << statusITS[3]
+                  << "statusITS4="          << statusITS[4]
+                  << "statusITS5="          << statusITS[5]
                   << "trd0="                << hasTrackletTRD[0]
                   << "trd1="                << hasTrackletTRD[1]
                   << "trd2="                << hasTrackletTRD[2]
                   << "trd3="                << hasTrackletTRD[3]
                   << "trd4="                << hasTrackletTRD[4]
                   << "trd5="                << hasTrackletTRD[5]
+                  << "TRDdEdxl0="           << trddEdxSum[0]
+                  << "TRDdEdxl1="           << trddEdxSum[1]
+                  << "TRDdEdxl2="           << trddEdxSum[2]
+                  << "TRDdEdxl3="           << trddEdxSum[3]
+                  << "TRDdEdxl4="           << trddEdxSum[4]
+                  << "TRDdEdxl5="           << trddEdxSum[5]
                   << "TOFsigmaEl="          << nSigmaTOF
                   << "TPCsigmaEl="          << nSigmaTPC
                   << "TPCdEdx="             << tPCdEdx
                   << "TRDlikeEl="           << likeEleTRD
                   << "TRDlikeEln="          << likeEleTRDn
+                  << "trdtruncmean1="       << trdtruncmean1
+                  << "trdtruncmean2="       << trdtruncmean2
                   << "dcaR="                << b[0]
                   << "dcaZ="                << b[1]
                   << "dca="                 << dca
                   << "dcaSR="               << dcaSR
                   << "dcaSZ="               << dcaSZ
                   << "dcaS="                << dcaS
+                  << "hfedcaR="             << hfeb[0]
+                  << "hfedcaZ="             << hfeb[1]
+                  << "hfedcacovR="          << hfebCov[0]
+                  << "hfedcacovZ="          << hfebCov[2]
                   << "vx="                  << vtx[0]
                   << "vy="                  << vtx[1]
-                  << "vz="                  << vtx[2] 
+                  << "vz="                  << vtx[2]
+                  << "tofdx="                << tofdx
+                 << "tofdz="                << tofdz
+                  << "statusTOFtracking="   << hasgoodTOF
+                  << "TOFmismatchlevel="    << mismatchlevel
+                  << "statusTRDtracking="   << hasTRDin
                   << "ncontrib="            << ncontrib
+                  << "mesonID="             << mesonID
+                  << "eR="                  << eR
+                  << "mesonR="              << mesonR
+                  << "eZ="                  << eZ
+                  << "mesonZ="              << mesonZ
+                  << "unique="              << unique
+                  << "mesonunique="         << mesonunique
+                  << "bgcategory="          << bgcategory
+                  << "mesonpt="             << mesonPt
+                  << "mesonMomPdg="         << mesonMomPdg
+                  << "mesonGMomPdg="         << mesonGMomPdg
+                  << "mesonGGMomPdg="         << mesonGGMomPdg
+                  << "mesonMomPt="         << mesonMomPt
                   << "\n";
   }
 }
 
 
 void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }
+//___________________________________________________________
+AliTrackReference *AliHFEdebugTreeTask::FindTrackReference(AliMCParticle *track, Float_t minRadius, Float_t maxRadius, Int_t detectorID)
+{
+  //
+  // Find the track reference
+  //
+  AliTrackReference *ref = NULL, *reftmp;
+  Float_t radius;
+  for(Int_t iref = 0; iref < track->GetNumberOfTrackReferences(); iref++){
+    reftmp = track->GetTrackReference(iref);
+    if(reftmp->DetectorId() != detectorID) continue;
+    radius = reftmp->R();
+    if(radius >= minRadius && radius < maxRadius){
+      ref = reftmp;
+      break;
+    } 
+    if(radius > maxRadius) break;
+  }
+  return ref;
+}
+
+//__________________________________________
+Int_t AliHFEdebugTreeTask::GetElecSourceMC(TParticle * const mcpart)
+{
+  // decay particle's origin 
+
+  if(!mcpart){
+    AliDebug(1, "no mcparticle, return\n");
+    return -1;
+  }
+
+  if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return AliHFEmcQA::kMisID;
+
+  Int_t origin = -1;
+  Bool_t isFinalOpenCharm = kFALSE;
+
+  Int_t iLabel = mcpart->GetFirstMother();
+  if (iLabel<0){
+    AliDebug(1, "Stack label is negative, return\n");
+    return -1;
+  }
+
+  AliMCParticle *mctrack = NULL;
+  Int_t tmpMomLabel=0;
+  if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
+  TParticle *partMother = mctrack->Particle();
+  TParticle *partMotherCopy = mctrack->Particle();
+  Int_t maPdgcode = partMother->GetPdgCode();
+
+  // if the mother is charmed hadron  
+  if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kCharm || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kCharm ) {
+
+    for (Int_t i=0; i<fNparents; i++){
+        if (abs(maPdgcode)==fParentSelect[0][i]){
+          isFinalOpenCharm = kTRUE;
+        }
+    }
+    if (!isFinalOpenCharm) return -1;
+
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<fgkMaxIter; i++){
+
+        Int_t jLabel = partMother->GetFirstMother();
+        if (jLabel == -1){
+          origin = AliHFEmcQA::kDirectCharm;
+          return origin;
+        }
+        if (jLabel < 0){ // safety protection
+          AliDebug(1, "Stack label is negative, return\n");
+          return -1;
+        }
+
+        // if there is an ancester
+        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
+        TParticle* grandMa = mctrack->Particle();
+        Int_t grandMaPDG = grandMa->GetPdgCode();
+
+        for (Int_t j=0; j<fNparents; j++){
+          if (abs(grandMaPDG)==fParentSelect[1][j]){
+            origin = AliHFEmcQA::kBeautyCharm;
+            return origin;
+          }
+        }
+
+        partMother = grandMa;
+    } // end of iteration 
+  } // end of if
+  else if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kBeauty || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kBeauty ) {
+    for (Int_t i=0; i<fNparents; i++){
+        if (abs(maPdgcode)==fParentSelect[1][i]){
+          origin = AliHFEmcQA::kDirectBeauty;
+          return origin;
+        }
+    }
+  } // end of if
+  else if ( abs(maPdgcode) == 22 ) { //conversion
+
+    tmpMomLabel = partMotherCopy->GetFirstMother();
+    if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
+    partMother = mctrack->Particle();
+    maPdgcode = partMother->GetPdgCode();
+    if ( abs(maPdgcode) == 111 ) {
+      origin = AliHFEmcQA::kGammaPi0;
+      return origin;
+    } 
+    else if ( abs(maPdgcode) == 221 ) {
+      origin = AliHFEmcQA::kGammaEta;
+      return origin;
+    } 
+    else if ( abs(maPdgcode) == 223 ) {
+      origin = AliHFEmcQA::kGammaOmega;
+      return origin;
+    } 
+    else if ( abs(maPdgcode) == 333 ) {
+      origin = AliHFEmcQA::kGammaPhi;
+      return origin;
+    }
+    else if ( abs(maPdgcode) == 331 ) {
+      origin = AliHFEmcQA::kGammaEtaPrime;
+      return origin; 
+    }
+    else if ( abs(maPdgcode) == 113 ) {
+      origin = AliHFEmcQA::kGammaRho0;
+      return origin;
+    }
+    else origin = AliHFEmcQA::kElse;
+    //origin = kGamma; // finer category above
+    return origin;
+
+  } // end of if
+  else if ( abs(maPdgcode) == 111 ) {
+    origin = AliHFEmcQA::kPi0;
+    return origin;
+  } // end of if
+  else if ( abs(maPdgcode) == 221 ) {
+    origin = AliHFEmcQA::kEta;
+    return origin;
+  } // end of if
+  else if ( abs(maPdgcode) == 223 ) {
+    origin = AliHFEmcQA::kOmega;
+    return origin;
+  } // end of if
+  else if ( abs(maPdgcode) == 333 ) {
+    origin = AliHFEmcQA::kPhi;
+    return origin;
+  } // end of if
+  else if ( abs(maPdgcode) == 331 ) {
+    origin = AliHFEmcQA::kEtaPrime;
+    return origin;
+  } // end of if
+  else if ( abs(maPdgcode) == 113 ) {
+    origin = AliHFEmcQA::kRho0;
+    return origin;
+  } // end of if
+  else{ 
+    origin = AliHFEmcQA::kElse;
+  }
+  return origin;
+}
+
+//______________________________________________________
+Bool_t AliHFEdebugTreeTask::CheckITSstatus( const AliESDtrack * const esdtrack, Int_t layer) const {
+  //
+  // Check whether ITS area is dead
+  //
+  Int_t itsStatus = 0;
+  Int_t det;
+  Float_t xloc, zloc;
+  esdtrack->GetITSModuleIndexInfo(layer, det, itsStatus, xloc, zloc);
+  Bool_t status;
+  switch(itsStatus){
+    case 2: status = kFALSE; break;
+    case 3: status = kFALSE; break;
+    case 7: status = kFALSE; break;
+    default: status = kTRUE;
+  }
+  return status;
+}
+