]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEdebugTreeTaskAOD.cxx
add EMCal trigger in eh analysis. add hists. in Raa study
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdebugTreeTaskAOD.cxx
index 384a836acea253134c200fe00560628fbc8c5386..cf5f692e0471784ce61a77108ca8d92d07641eb8 100644 (file)
@@ -26,6 +26,7 @@
 #include "AliCentrality.h"
 #include "AliAODTrack.h"
 #include "AliAODEvent.h"
+#include "AliAODPid.h"
 #include "AliHFEcuts.h"
 #include "AliHFEextraCuts.h"
 #include "AliInputEventHandler.h"
@@ -96,7 +97,7 @@ void AliHFEdebugTreeTaskAOD::UserCreateOutputObjects(){
   //printf("test\n");
   fDebugTree = new TTreeSRedirector(fFilename.Data());
 
 //printf("testa\n");
// printf("testa\n");
   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
   //printf("testb\n");
   
@@ -129,13 +130,13 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   AliPIDResponse *pid = NULL;
   AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
   if(handler){
-    //printf("testb\n");
+//    printf("testb\n");
     pid = handler->GetPIDResponse();
   } else {
     AliError("No Handler");
   }
   if(!pid){
-    printf("testc\n");
//   printf("testc\n");
     AliError("No PID response");
     return;
   }
@@ -144,22 +145,25 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     return;
   }
 
+  AliAODTrack copyTrack;
+
   // MC info
   Bool_t mcthere = kTRUE;
   AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
-  if(!aodE){ 
+  if(!aodE){
+ //        printf("testd\n");
     AliError("No AOD Event");
     return;
   }
   fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
   if(!fAODMCHeader){ 
-    mcthere = kFALSE;
-    return;
+      mcthere = kFALSE;
//   return;
   }
   fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
   if(!fAODArrayMCInfo){ 
-    mcthere = kFALSE;
-    return;
+      mcthere = kFALSE;
+  //  return;
   }
   else {
     fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
@@ -178,6 +182,9 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   }
   fExtraCuts->SetRecEventInfo(fInputEvent);
 
+
   // Get run number
   Int_t run = fInputEvent->GetRunNumber();
 
@@ -199,6 +206,16 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   AliCentrality *hicent = fInputEvent->GetCentrality();
   centrality = hicent->GetCentralityPercentile("V0M");
 
+  // Store event selection variables
+  (*fDebugTree) << "EventDebug"
+                << "Centrality="              << centrality
+                << "VertexZ="                 << vtx[2]
+                << "NumberOfContributors="    << ncontrib
+                << "run="                 << run
+                << "\n";
+
+
+
   // Look for kink mother
   Int_t numberofvertices = aodE->GetNumberOfVertices();
   Double_t listofmotherkink[numberofvertices];
@@ -219,6 +236,73 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     }
   }
   //printf("Number of kink mother in the events %d\n",numberofmotherkink);
+
+  //
+  // Loop on MC tracks only
+  //
+  AliAODMCParticle *mctrack = NULL;
+  // Monte-Carlo info
+  Double_t eR,vx,vy,vz;
+  Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
+  Int_t source,pdg,signal;
+  if(mcthere) {
+    for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
+      mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
+      if(!mctrack) continue;
+      signal = 0;
+      if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
+      // Kinematics
+      chargemc = mctrack->Charge() > 0. ? 1. : -1.;
+      momentummc = mctrack->P() * chargemc;
+      transversemomentummc = mctrack->Pt() * chargemc;
+      etamc = mctrack->Eta();
+      phimc = mctrack->Phi();
+      pdg = mctrack->GetPdgCode();
+      
+      // Get Production Vertex in radial direction
+      vx = mctrack->Xv();
+      vy = mctrack->Yv(); 
+      vz = mctrack->Zv(); 
+      eR = TMath::Sqrt(vx*vx+vy*vy);
+
+      // Get Mother PDG code of the particle
+      Int_t motherPdg = 0;
+      Int_t motherlabel = TMath::Abs(mctrack->GetMother());
+      if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
+        AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
+        if(mother) motherPdg = mother->GetPdgCode();
+      }
+      
+      // derive source
+      source = 5;
+      if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
+      else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
+      else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
+      else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
+      else if(TMath::Abs(pdg) == 11) source = 4;
+      else source = 5;
+
+      (*fDebugTree) << "MCDebug"
+                    << "centrality="          << centrality
+                    << "MBtrigger="           << isMBTrigger 
+                    << "CentralTrigger="      << isCentralTrigger
+                    << "SemicentralTrigger="  << isSemicentralTrigger
+                    << "EMCALtrigger="        << isEMCALTrigger
+                    << "run="                 << run
+                    << "p="                   << momentummc
+                    << "pt="                  << transversemomentummc
+                    << "eta="                 << etamc
+                    << "phi="                 << phimc
+                    << "pdg="                 << pdg
+                   << "px="                  << vx
+                   << "py="                  << vy
+                   << "pz="                  << vz
+                    << "ProductionVertex="    << eR
+                    << "motherPdg="           << motherPdg
+                    << "source="              << source
+                    << "\n";
+    }
+  }
   
   // Common variables
   Double_t charge, eta, phi, momentum, momentumTPC, transversemomentum;
@@ -228,21 +312,28 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   //
 
   TArrayI *arraytrack = new TArrayI(fInputEvent->GetNumberOfTracks());
+  Int_t counterdc=0;
   
   AliAODTrack *track = 0x0;
-  AliAODMCParticle *mctrack = NULL;
   for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
     // fill the tree
     track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
     if(!track) continue;
     // Cut track (Only basic track cuts)
-    //printf("testv\n");
+    // printf("testv\n");
     if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     //
     //printf("testu\n");
     Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
     Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
     Double_t tPCdEdx = track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.;
+    // Eta correction
+    copyTrack.~AliAODTrack();
+    new(&copyTrack) AliAODTrack(*track);
+    if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
+    if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
+    Double_t nSigmaTPCcorr = pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron);
+    Double_t tPCdEdxcorr = copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.;
    
     // Kinematics
     charge = track->Charge() > 0 ? 1. : -1.;
@@ -258,6 +349,8 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) itsrefit = 1;
     Int_t tpcrefit=0;
     if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) tpcrefit = 1;
+    Int_t tofpid=0;
+    if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) tofpid = 1;
 
     // ITS number of clusters
     UChar_t nclustersITS = track->GetITSNcls();
@@ -279,8 +372,29 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     // TRD clusters and tracklets
     UChar_t nclustersTRD = track->GetTRDncls();
     UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
-    Int_t   nslicesTRD = track->GetNumberOfTRDslices();
+    Int_t   nslicesTRD = track->GetDetPid()->GetTRDnSlices();
     Int_t   chi2TRD = track->GetTRDchi2();
+    AliAODPid* aodpid= track->GetDetPid();
+    Double_t* arraytrdsignals;
+    arraytrdsignals=aodpid->GetTRDslices();
+    Int_t nslicetemp=0;
+    Int_t trdlayer[6];
+    for(Int_t iplane = 0; iplane < 6; iplane++){
+       trdlayer[iplane]=0;
+    }
+
+    for(Int_t iplane = 0; iplane < 6; iplane++){
+       nslicetemp=0;
+       for(Int_t b=(iplane*8);b<((iplane*8)+8);b++)
+       {
+           if(ntrackletsTRDPID>0)
+           {
+               if(arraytrdsignals[b]>0.001) nslicetemp++;
+           }
+       }
+       if(nslicetemp > 0) trdlayer[iplane]=1;
+    }
+
     // ITS and TRD acceptance maps
     UChar_t itsPixel = track->GetITSClusterMap();
     Bool_t statusL0 = kFALSE;
@@ -311,14 +425,15 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
 
     // Double counted
     Int_t doublec = 0;
-    for(Int_t l=0; l < itrack; l++){
+    for(Int_t l=0; l < counterdc; l++){
       Int_t iTrack2 = arraytrack->At(l);
       if(iTrack2==id) doublec=1;
     }
     //printf("Doublec %d\n",doublec);
 
     // Add the id at this place
-    arraytrack->AddAt(id,itrack);
+    arraytrack->AddAt(id,counterdc);
+    counterdc++;
 
     // Filter
     Int_t filter[20];
@@ -354,10 +469,6 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
 
     //printf("track\n");
 
-    // Monte-Carlo info
-    Double_t eR,vx,vy,vz;
-    Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
-    Int_t source,pdg,signal;
     signal = 0;
     if(mcthere){
       Int_t label = TMath::Abs(track->GetLabel());
@@ -367,8 +478,8 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
       if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
       // Kinematics
       chargemc = mctrack->Charge() > 0. ? 1. : -1.;
-      momentummc = mctrack->P() * charge;
-      transversemomentummc = mctrack->Pt() * charge;
+      momentummc = mctrack->P() * chargemc;
+      transversemomentummc = mctrack->Pt() * chargemc;
       etamc = mctrack->Eta();
       phimc = mctrack->Phi();
       pdg = mctrack->GetPdgCode();
@@ -422,6 +533,7 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
                   << "phi="                 << phi
                  << "itsrefit="            << itsrefit
                  << "tpcrefit="            << tpcrefit
+                 << "tofpid="              << tofpid
                  << "nclustersTPC="        << nclustersTPCfit
                  << "nclustersTPCall="     << nclustersTPCall
                   << "nclustersTPCPID="     << nclustersTPCPID
@@ -433,12 +545,20 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
                  << "nclustersTRD="        << nclustersTRD
                  << "ntrackletsTRD="       << ntrackletsTRDPID
                  << "nslicesTRD="          << nslicesTRD
+                  << "trd0="                << trdlayer[0]
+                  << "trd1="                << trdlayer[1]
+                  << "trd2="                << trdlayer[2]
+                  << "trd3="                << trdlayer[3]
+                  << "trd4="                << trdlayer[4]
+                  << "trd5="                << trdlayer[5]
                   << "chi2TRD="             << chi2TRD
                  << "statusITS0="          << statusL0
                   << "statusITS1="          << statusL1
                  << "TOFsigmaEl="          << nSigmaTOF
                   << "TPCsigmaEl="          << nSigmaTPC
+                 << "TPCsigmaElcorr="      << nSigmaTPCcorr
                   << "TPCdEdx="             << tPCdEdx
+                  << "TPCdEdxcorr="         << tPCdEdxcorr
                  << "dcaR="                << dcaxy
                   << "dcaZ="                << dcaz
                  << "kinkdaughter="        << kink