]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/TRD/AliTRDpidRefMaker.cxx
Adding 2 triggers to list of known unsupported triggers
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMaker.cxx
index b6d01e1c20b2e47442d3f22325a50f3420c1430f..c87f2f8212633e1a51a7c7ecbcd3b282eff6b18a 100644 (file)
@@ -8,12 +8,16 @@
 #include "AliLog.h"
 #include "AliESDtrack.h"
 #include "AliTrackReference.h"
+#include "AliAnalysisManager.h"
 
 #include "AliTRDReconstructor.h"
 #include "AliTRDtrackV1.h"
 #include "AliTRDseedV1.h"
 #include "AliTRDpidRefMaker.h"
+#include "AliTRDinfoGen.h"
+#include "info/AliTRDeventInfo.h"
 #include "info/AliTRDv0Info.h"
+#include "info/AliTRDpidInfo.h"
 
 
 // Defines and implements basic functionality for building reference data for TRD PID. 
 // 
 
 ClassImp(AliTRDpidRefMaker)
-ClassImp(AliTRDpidRefMaker::AliTRDpidRefData)
-ClassImp(AliTRDpidRefMaker::AliTRDpidRefDataArray)
 
 //________________________________________________________________________
-AliTRDpidRefMaker::AliTRDpidRefDataArray::AliTRDpidRefDataArray() :
-  fNtracklets(0)
+AliTRDpidRefMaker::AliTRDpidRefMaker() 
+  :AliTRDrecoTask()
+  ,fV0s(NULL)
   ,fData(NULL)
+  ,fInfo(NULL)
+  ,fPIDdataArray(NULL)
+  ,fRefPID(kV0)
+  ,fRefP(kRec)
+  ,fFreq(1.)
+  ,fP(-1.)
+  ,fPthreshold(0.)
 {
-  // Constructor of data array
-  fData = new AliTRDpidRefData[AliTRDgeometry::kNlayer];
-}
-
-//________________________________________________________________________
-AliTRDpidRefMaker::AliTRDpidRefDataArray::~AliTRDpidRefDataArray()
-{
-  // Destructor
-  delete [] fData;
-}
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::AliTRDpidRefDataArray::PushBack(Int_t ly, Int_t p, Float_t *dedx)
-{
-// Add PID data to the end of the array 
-  fData[fNtracklets].fPLbin= (ly<<4) | (p&0xf);
-  memcpy(fData[fNtracklets].fdEdx, dedx, 8*sizeof(Float_t));
-  fNtracklets++;
-}
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::AliTRDpidRefDataArray::Reset()
-{
-// Reset content
-
-  if(!fNtracklets) return;
-  while(fNtracklets--){
-    fData[fNtracklets].fPLbin = 0xff;
-    memset(fData[fNtracklets].fdEdx, 0, 8*sizeof(Float_t));
-  }
-  fNtracklets=0;
+  //
+  // Default constructor
+  //
+  SetNameTitle("PIDrefMaker", "PID Reference Maker");
 }
 
 //________________________________________________________________________
 AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title) 
   :AliTRDrecoTask(name, title)
-  ,fReconstructor(NULL)
   ,fV0s(NULL)
   ,fData(NULL)
   ,fInfo(NULL)
   ,fPIDdataArray(NULL)
-  ,fRefPID(kMC)
-  ,fRefP(kMC)
-  ,fPIDbin(0xff)
+  ,fRefPID(kV0)
+  ,fRefP(kRec)
   ,fFreq(1.)
   ,fP(-1.)
   ,fPthreshold(0.5)
@@ -90,14 +71,12 @@ AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
   // Default constructor
   //
 
-  fReconstructor = new AliTRDReconstructor();
-  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
-  memset(fdEdx, 0, 10*sizeof(Float_t));
+  memset(fdEdx, 0, AliTRDpidUtil::kNNslices*sizeof(Float_t));
   memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
 
-  DefineInput(1, TObjArray::Class());
-  DefineInput(2, TObjArray::Class());
-  DefineOutput(1, TTree::Class());
+  DefineInput(3, TObjArray::Class()); // v0 list
+  DefineInput(4, TObjArray::Class()); // pid info list
+  DefineOutput(2, TTree::Class());
 }
 
 
@@ -105,7 +84,6 @@ AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
 AliTRDpidRefMaker::~AliTRDpidRefMaker() 
 {
   if(fPIDdataArray) delete fPIDdataArray;
-  if(fReconstructor) delete fReconstructor;
 }
 
 //________________________________________________________________________
@@ -126,20 +104,11 @@ Float_t* AliTRDpidRefMaker::CookdEdx(AliTRDseedV1 *trklt)
 
 
 //________________________________________________________________________
-void AliTRDpidRefMaker::ConnectInputData(Option_t *opt)
-{
-  AliTRDrecoTask::ConnectInputData(opt);
-  fV0s  = dynamic_cast<TObjArray*>(GetInputData(1));
-  fInfo = dynamic_cast<TObjArray*>(GetInputData(2));
-} 
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::CreateOutputObjects()
+void AliTRDpidRefMaker::UserCreateOutputObjects()
 {
   // Create histograms
   // Called once
 
-  OpenFile(0, "RECREATE");
   fContainer = new TObjArray();
   fContainer->SetName(Form("Moni%s", GetName()));
 
@@ -150,85 +119,90 @@ void AliTRDpidRefMaker::CreateOutputObjects()
   for(Int_t is=AliPID::kSPECIES; is--;) ax->SetBinLabel(is+1, AliPID::ParticleShortName(is));
   h2->GetYaxis()->SetTitle("P bins");
   h2->GetYaxis()->SetNdivisions(511);
-
   fContainer->AddAt(h2, 0);
+  PostData(1, fContainer);
 
-  fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
-  fData->Branch("s", &fPIDbin, "s/b");
+  OpenFile(2);
+  fData = new TTree("RefPID", "Reference data for PID");
   fData->Branch("data", &fPIDdataArray);
+  PostData(2, fData);
 }
 
 //________________________________________________________________________
-void AliTRDpidRefMaker::Exec(Option_t *) 
+void AliTRDpidRefMaker::UserExec(Option_t *) 
 {
   // Main loop
   // Called for each event
+  Int_t ev((Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry());
+  if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))){
+    AliDebug(3, Form("Missing tracks container in ev %d", ev));
+    return;
+  }
+  if(!(fEvent = dynamic_cast<AliTRDeventInfo*>(GetInputData(2)))){
+    AliDebug(3, Form("Missing Event Info container in ev %d", ev));
+    return;
+  }
+  if(!(fV0s    = dynamic_cast<TObjArray*>(GetInputData(3)))){
+    AliDebug(3, Form("Missing v0 container in ev %d", ev)); 
+    return;
+  }
+  if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(4)))){
+    AliDebug(3, Form("Missing pid info container in ev %d", ev)); 
+    return;
+  }
 
+  AliDebug(1, Form("Entries: Ev[%d] Tracks[%d] V0[%d] PID[%d]", ev, fTracks->GetEntriesFast(), fV0s->GetEntriesFast(), fInfo->GetEntriesFast()));
   AliTRDtrackInfo     *track = NULL;
-  AliTRDtrackV1    *trackTRD = NULL;
+  //AliTRDtrackV1    *trackTRD = NULL;
   AliTrackReference     *ref = NULL;
   const AliTRDtrackInfo::AliESDinfo *infoESD = NULL;
   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
     if(!track->HasESDtrack()) continue;
-    trackTRD = track->GetTrack();
+    //trackTRD = track->GetTrack();
     infoESD  = track->GetESDinfo();
     Double32_t *infoPID = infoESD->GetSliceIter();
     Int_t n = infoESD->GetNSlices() - AliTRDgeometry::kNlayer;
+    if(n==0){
+      AliWarning(Form("dEdx info missing in ESD track %d", itrk));
+      continue;
+    }
     Double32_t *p = &infoPID[n];
-
+    AliDebug(4, Form("n[%d] p[GeV/c]{%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f}", n, p[0], p[1], p[2], p[3], p[4], p[5]));
 
     ULong_t status = track->GetStatus();
-    if(!(status&AliESDtrack::kTPCout)) continue;
+    if(!(status&AliESDtrack::kTRDpid)) continue;
 
     // fill the pid information
-    SetRefPID(fRefPID, track, fPID);
-
+    SetRefPID(fRefPID, track, infoESD, fPID);
+    // get particle type
+    Int_t idx(TMath::Max(Long64_t(0), TMath::LocMax(AliPID::kSPECIES, fPID))); 
+    if(fPID[idx]<1.e-5) continue;
+    
     // prepare PID data array
     if(!fPIDdataArray){ 
-      fPIDdataArray = new AliTRDpidRefDataArray();
+      fPIDdataArray = new AliTRDpidInfo();
     } else fPIDdataArray->Reset();
+    fPIDdataArray->SetPID(idx);
 
     // fill PID information
     for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
 
       // fill P & dE/dx information
-      if(HasFriends()){ // from TRD track
-        if(!trackTRD) continue;
-        AliTRDseedV1 *trackletTRD(NULL);
-        trackTRD -> SetReconstructor(fReconstructor);
-        if(!(trackletTRD = trackTRD -> GetTracklet(ily))) continue;
-        if(!CheckQuality(trackletTRD)) continue;
-        if(!CookdEdx(trackletTRD)) continue;
-
-        // fill momentum information
-        fP = 0.;
-        switch(fRefP){
-        case kMC:
-          if(!(ref = track->GetTrackRef(trackletTRD))) continue;
-          fP = ref->P();
-          break;
-        case kRec:
-          fP = trackletTRD->GetMomentum();
-          break;
-        default: continue;
-        }
-      } else { // from ESD track
-        // fill momentum information
-        switch(fRefP){
-        case kMC:
-          if(!(ref = track->GetTrackRef(ily))) continue;
-          fP = ref->P();
-          break;
-        case kRec:
-          fP = p[ily];
-          break;
-        default: continue;
-        } 
-        Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
-        for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
+      switch(fRefP){
+      case kMC:
+        if(!(ref = track->GetTrackRef(ily))) continue;
+        fP = ref->P();
+        break;
+      case kRec:
+        fP = p[ily];
+        break;
+      default:
+        continue;
       }
-
+      Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
+      for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
+      
       // momentum threshold
       if(fP < fPthreshold) continue;
 
@@ -238,9 +212,6 @@ void AliTRDpidRefMaker::Exec(Option_t *)
 
     Fill();
   }
-
-  PostData(0, fContainer);
-  PostData(1, fData);
 }
 
 
@@ -249,16 +220,15 @@ void AliTRDpidRefMaker::Fill()
 {
 // Fill data tree
 
-  if(!fPIDdataArray->fNtracklets) return;
-  fPIDbin = TMath::LocMax(AliPID::kSPECIES, fPID); // get particle type
-// Fill data tree
+  if(!fPIDdataArray->GetNtracklets()) return;
+  // Fill data tree
   fData->Fill();
 
   
   // fill monitor
-  for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
-    Int_t pBin = fPIDdataArray->fData[ily].fPLbin & 0xf;
-    ((TH2*)fContainer->At(0))->Fill(fPIDbin, pBin);
+  for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
+    Int_t pBin = fPIDdataArray->GetData(itrklt)->Momentum();
+    ((TH2*)fContainer->At(0))->Fill(fPIDdataArray->GetPID(), pBin);
   }
 }
 
@@ -266,12 +236,11 @@ void AliTRDpidRefMaker::Fill()
 void AliTRDpidRefMaker::LinkPIDdata() 
 {
 // Link data tree to data members
-  fData->SetBranchAddress("s", &fPIDbin);
   fData->SetBranchAddress("data", &fPIDdataArray);
 }
 
 //________________________________________________________________________
-void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid) 
+void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, AliTRDtrackInfo *track, const AliTRDtrackInfo::AliESDinfo *infoESD, Float_t *pid) 
 {
 // Fill the reference PID values "pid" from "source" object
 // according to the option "select". Possible options are
@@ -279,20 +248,18 @@ void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Fl
 // - kMC  - MC truth [default]
 // - kRec - outside detectors
 
-
-  memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
+  if(!track){
+    AliError("No trackInfo found");
+    return;
+  }
+  memset(pid, 0, AliPID::kSPECIES*sizeof(Float_t));
   switch(select){ 
   case kV0:
     {
-      AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
-      if(!track){
-        AliError("No trackInfo found");
-        return;
-      }
-
-      //Get V0 PID decisions from the AliTRDv0Info for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) :
-      AliTRDv0Info v0;
-      for(Int_t is=AliPID::kSPECIES; is--;) fPID[is] = v0.GetV0PID(is, track);
+      //Get V0 PID decisions for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) :
+      if(!infoESD->HasV0()) return;
+      const Int_t *v0pid=infoESD->GetV0pid();
+      for(Int_t is=AliPID::kSPECIES; is--;){ pid[is] = (Float_t)v0pid[is];}
     }
     break;
   case kMC:
@@ -300,37 +267,33 @@ void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Fl
       AliError("Could not retrive reference PID from MC");
       return;
     }
-    {
-      AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
-      switch(track->GetPDG()){
-      case kElectron:
-      case kPositron:
-        fPID[AliPID::kElectron] = 1.;
-        break;
-      case kMuonPlus:
-      case kMuonMinus:
-        fPID[AliPID::kMuon] = 1.;
-        break;
-      case kPiPlus:
-      case kPiMinus:
-        fPID[AliPID::kPion] = 1.;
-        break;
-      case kKPlus:
-      case kKMinus:
-        fPID[AliPID::kKaon] = 1.;
-        break;
-      case kProton:
-      case kProtonBar:
-        fPID[AliPID::kProton] = 1.;
-        break;
-      }
+    switch(track->GetPDG()){
+    case kElectron:
+    case kPositron:
+      pid[AliPID::kElectron] = 1.;
+      break;
+    case kMuonPlus:
+    case kMuonMinus:
+      pid[AliPID::kMuon] = 1.;
+      break;
+    case kPiPlus:
+    case kPiMinus:
+      pid[AliPID::kPion] = 1.;
+      break;
+    case kKPlus:
+    case kKMinus:
+      pid[AliPID::kKaon] = 1.;
+      break;
+    case kProton:
+    case kProtonBar:
+      pid[AliPID::kProton] = 1.;
+      break;
     }
     break;
   case kRec:
     { 
-      AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
       AliTRDtrackV1 *trackTRD = track->GetTrack();
-      trackTRD -> SetReconstructor(fReconstructor);
+      trackTRD -> SetReconstructor(AliTRDinfoGen::Reconstructor());
       //fReconstructor -> SetOption("nn");
       trackTRD -> CookPID();
       for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
@@ -343,12 +306,12 @@ void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Fl
     AliWarning("PID reference source not implemented");
     return;
   }
-  AliDebug(4, Form("Ref PID [%] : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]"
-    ,AliPID::ParticleShortName(0), 1.e2*fPID[0]
-    ,AliPID::ParticleShortName(1), 1.e2*fPID[1]
-    ,AliPID::ParticleShortName(2), 1.e2*fPID[2]
-    ,AliPID::ParticleShortName(3), 1.e2*fPID[3]
-    ,AliPID::ParticleShortName(4), 1.e2*fPID[4]
+  AliDebug(4, Form("Ref PID : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]"
+    ,AliPID::ParticleShortName(0), 1.e2*pid[0]
+    ,AliPID::ParticleShortName(1), 1.e2*pid[1]
+    ,AliPID::ParticleShortName(2), 1.e2*pid[2]
+    ,AliPID::ParticleShortName(3), 1.e2*pid[3]
+    ,AliPID::ParticleShortName(4), 1.e2*pid[4]
   ));
 }