]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/TRD/AliTRDcheckPID.cxx
set reco param on an event by event basis
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckPID.cxx
index 7d743dbe42f866e69c29e6d9d38a4fdf948f67d7..cd37e4adeb3bb902c5cb5993f01e5c964b90f37e 100644 (file)
 #include "Cal/AliTRDCalPID.h"
 #include "Cal/AliTRDCalPIDNN.h"
 #include "AliTRDcheckPID.h"
+#include "AliTRDinfoGen.h"
 #include "info/AliTRDtrackInfo.h"
 #include "info/AliTRDpidInfo.h"
+#include "info/AliTRDv0Info.h"
 
 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
 
@@ -60,10 +62,10 @@ ClassImp(AliTRDcheckPID)
 //________________________________________________________________________
 AliTRDcheckPID::AliTRDcheckPID() 
   :AliTRDrecoTask()
-  ,fReconstructor(NULL)
   ,fUtil(NULL)
   ,fGraph(NULL)
   ,fPID(NULL)
+  ,fV0s(NULL)
   ,fMomentumAxis(NULL)
   ,fMinNTracklets(AliTRDgeometry::kNlayer)
   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
@@ -71,16 +73,17 @@ AliTRDcheckPID::AliTRDcheckPID()
   //
   // Default constructor
   //
-  SetNameTitle("checkPID", "Check TRD PID");
+  SetNameTitle("TRDcheckPID", "Check TRD PID");
+  LocalInit();
 }
 
 //________________________________________________________________________
 AliTRDcheckPID::AliTRDcheckPID(char* name ) 
   :AliTRDrecoTask(name, "Check TRD PID")
-  ,fReconstructor(NULL)
   ,fUtil(NULL)
   ,fGraph(NULL)
   ,fPID(NULL)
+  ,fV0s(NULL)
   ,fMomentumAxis(NULL)
   ,fMinNTracklets(AliTRDgeometry::kNlayer)
   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
@@ -89,8 +92,18 @@ AliTRDcheckPID::AliTRDcheckPID(char* name )
   // Default constructor
   //
 
-  fReconstructor = new AliTRDReconstructor();
-  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
+  LocalInit();
+  InitFunctorList();
+
+  DefineInput(3, TObjArray::Class());  // v0 list
+  DefineOutput(2, TObjArray::Class()); // pid info list
+}
+
+
+//________________________________________________________________________
+void AliTRDcheckPID::LocalInit() 
+{
+// Initialize working data
 
   // Initialize momentum axis with default values
   Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
@@ -101,19 +114,13 @@ AliTRDcheckPID::AliTRDcheckPID(char* name )
   memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
 
   fUtil = new AliTRDpidUtil();
-  InitFunctorList();
-
-  DefineOutput(2, TObjArray::Class());
-  
 }
 
-
 //________________________________________________________________________
 AliTRDcheckPID::~AliTRDcheckPID() 
 {
   if(fPID){fPID->Delete(); delete fPID;}
   if(fGraph){fGraph->Delete(); delete fGraph;}
-  if(fReconstructor) delete fReconstructor;
   if(fUtil) delete fUtil;
 }
 
@@ -124,11 +131,10 @@ void AliTRDcheckPID::UserCreateOutputObjects()
   // Create histograms
   // Called once
 
-  OpenFile(1, "RECREATE");
-  fContainer = Histos();
-
+  AliTRDrecoTask::UserCreateOutputObjects();
   fPID = new TObjArray();
   fPID->SetOwner(kTRUE);
+  PostData(2, fPID);
 }
 
 //________________________________________________________
@@ -138,11 +144,10 @@ void AliTRDcheckPID::UserExec(Option_t *opt)
   // Execution part
   //
 
+  fV0s = dynamic_cast<TObjArray *>(GetInputData(3));
   fPID->Delete();
 
   AliTRDrecoTask::UserExec(opt);
-
-  PostData(2, fPID);
 }
 
 
@@ -155,7 +160,7 @@ TObjArray * AliTRDcheckPID::Histos(){
   if(fContainer) return fContainer;
 
   Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins(); 
-  fContainer = new TObjArray(); fContainer->Expand(kNPlots);
+  fContainer = new TObjArray(); fContainer->Expand(kNPlots); fContainer->SetOwner(kTRUE);
 
   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
   TH1 *h = NULL;
@@ -198,15 +203,15 @@ TObjArray * AliTRDcheckPID::Histos(){
   fPH->SetOwner(); fPH->SetName("PH");
   fContainer->AddAt(fPH, kPH);
   if(!(h = (TProfile2D*)gROOT->FindObject("PHT"))){
-    h = new TProfile2D("PHT", "", 
+    h = new TProfile2D("PHT", "<PH>(tb);p*species;tb [100*ns];entries", 
       xBins, -0.5, xBins - 0.5,
-      AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
+      AliTRDseedV1::kNtb, -0.5, AliTRDseedV1::kNtb - 0.5);
   } else h->Reset();
   fPH->AddAt(h, 0);
   if(!(h = (TProfile2D*)gROOT->FindObject("PHX"))){
-    h = new TProfile2D("PHX", "", 
+    h = new TProfile2D("PHX", "<PH>(x);p*species;x_{drift} [cm];entries", 
       xBins, -0.5, xBins - 0.5,
-      AliTRDtrackerV1::GetNTimeBins(), 0., .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght());
+      40, 0., 4.5);
   } else h->Reset();
   fPH->AddAt(h, 1);
 
@@ -238,6 +243,13 @@ TObjArray * AliTRDcheckPID::Histos(){
   } else h->Reset();
   fContainer->AddAt(h, kNTracklets);
 
+  // V0 performance
+  if(!(h = (TH1F*)gROOT->FindObject("nV0"))){
+    h = new TH1F("nV0", "V0s/track;v0/track;entries", 
+      6, -0.5, 5.5);
+  } else h->Reset();
+  fContainer->AddAt(h, kV0);
+
   return fContainer;
 }
 
@@ -262,8 +274,8 @@ Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
 {
 // Documentation to come
 
-  track -> SetReconstructor(fReconstructor);
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
+  track -> SetReconstructor(AliTRDinfoGen::Reconstructor());
+  (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
   track -> CookPID();
 
   if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion)  + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
@@ -309,7 +321,7 @@ TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
   if(!CheckTrackQuality(fkTrack)) return NULL;
 
   AliTRDtrackV1 cTrack(*fkTrack);
-  cTrack.SetReconstructor(fReconstructor);
+  cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 
   Int_t pdg = 0;
   Float_t momentum = 0.;
@@ -324,7 +336,7 @@ TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
   }
   if(!IsInRange(momentum)) return NULL;
 
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
+  (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
   cTrack.CookPID();
   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
@@ -374,7 +386,7 @@ TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
   if(!CheckTrackQuality(fkTrack)) return NULL;
   
   AliTRDtrackV1 cTrack(*fkTrack);
-  cTrack.SetReconstructor(fReconstructor);
+  cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
 
   Int_t pdg = 0;
   Float_t momentum = 0.;
@@ -388,7 +400,7 @@ TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
   }
   if(!IsInRange(momentum)) return NULL;
 
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
+  (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
   cTrack.CookPID();
   if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
 
@@ -447,7 +459,7 @@ TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
   } else {
     //AliWarning("No MC info available!");
     AliTRDtrackV1 cTrack(*fkTrack);
-    cTrack.SetReconstructor(fReconstructor);
+    cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
     momentum = cTrack.GetMomentum(0);
     pdg = CalcPDG(&cTrack);
   }
@@ -468,8 +480,6 @@ TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
       return NULL;
     }
 
-    AliTRDpidInfo *pid = new AliTRDpidInfo();
-    fPID->Add(pid);
     hPID->Fill(FindBin(species, momentum), pidESD[is]);
   }
   return hPID;
@@ -499,7 +509,7 @@ TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
   }
 
   AliTRDtrackV1 cTrack(*fkTrack);
-  cTrack.SetReconstructor(fReconstructor);
+  cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
   Int_t pdg = 0;
   Float_t momentum = 0.;
   if(fkMC){
@@ -510,25 +520,33 @@ TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
     momentum = cTrack.GetMomentum(0);
     pdg = CalcPDG(&cTrack);
   }
-  Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   if(!IsInRange(momentum)) return NULL;
 
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
-//   (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
-  Float_t sumdEdx = 0;
-  Int_t iBin = FindBin(species, momentum);
+  // Init exchange container
+  Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
+  AliTRDpidInfo *pid = new AliTRDpidInfo(s);
+
+  (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
+
+  Float_t sumdEdx(0.);
+  Int_t iBin = FindBin(s, momentum);
   AliTRDseedV1 *tracklet = NULL;
-  for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-    sumdEdx = 0;
-    tracklet = cTrack.GetTracklet(iChamb);
+  for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
+    tracklet = cTrack.GetTracklet(ily);
     if(!tracklet) continue;
     tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
+
+    // fill exchange container
+    pid->PushBack(tracklet->GetPlane(), 
+                  AliTRDpidUtil::GetMomentumBin(tracklet->GetMomentum()), tracklet->GetdEdx());
+
+    sumdEdx = 0.;
     for(Int_t i = AliTRDpidUtil::kNNslices; i--;) sumdEdx += tracklet->GetdEdx()[i];
-//     tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
-//     for(Int_t i = 3; i--;) sumdEdx += tracklet->GetdEdx()[i];
     sumdEdx /= AliTRDCalPIDNN::kMLPscale;
     hdEdx -> Fill(iBin, sumdEdx);
   }
+  fPID->Add(pid);
+
   return hdEdx;
 }
 
@@ -555,7 +573,7 @@ TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
   }
 
   AliTRDtrackV1 cTrack(*fkTrack);
-  cTrack.SetReconstructor(fReconstructor);
+  cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
   Int_t pdg = 0;
   Float_t momentum = 0.;
   if(fkMC){
@@ -568,7 +586,7 @@ TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
   }
   if(!IsInRange(momentum)) return NULL;
 
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
+  (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
   Int_t iMomBin = fMomentumAxis->FindBin(momentum);
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   Float_t *fdEdx;
@@ -619,7 +637,7 @@ TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
   } else {
     //AliWarning("No MC info available!");
     AliTRDtrackV1 cTrack(*fkTrack);
-    cTrack.SetReconstructor(fReconstructor);
+    cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
     momentum = cTrack.GetMomentum(0);
     pdg = CalcPDG(&cTrack);
   }
@@ -633,10 +651,10 @@ TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
     tracklet = fkTrack->GetTracklet(iChamb);
     if(!tracklet) continue;
     Float_t x0 = tracklet->GetX0(); 
-    for(Int_t ic = 0; ic < AliTRDtrackerV1::GetNTimeBins(); ic++){
+    for(Int_t ic = 0; ic < AliTRDseedV1::kNclusters; ic++){
       if(!(cluster = tracklet->GetClusters(ic))) continue;
       hPHT -> Fill(iBin, cluster->GetLocalTimeBin(), TMath::Abs(cluster->GetQ()));
-      hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
+      if(ic<AliTRDseedV1::kNtb) hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
     }
   }
   return hPHT;
@@ -673,7 +691,7 @@ TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
   } else {
     //AliWarning("No MC info available!");
     AliTRDtrackV1 cTrack(*fkTrack);
-    cTrack.SetReconstructor(fReconstructor);
+    cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
     momentum = cTrack.GetMomentum(0);
     pdg = CalcPDG(&cTrack);
   }
@@ -711,7 +729,7 @@ TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
   }
 
   AliTRDtrackV1 cTrack(*fkTrack);
-  cTrack.SetReconstructor(fReconstructor);
+  cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
   Int_t pdg = 0;
   Float_t momentum = 0.;
   if(fkMC){
@@ -719,7 +737,7 @@ TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
     pdg = fkMC->GetPDG();
   } else {
     //AliWarning("No MC info available!");
-    momentum = cTrack.GetMomentum(0);
+    momentum = cTrack.GetMomentum();
     pdg = CalcPDG(&cTrack);
   }
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
@@ -760,7 +778,7 @@ TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
   } else {
     //AliWarning("No MC info available!");
     AliTRDtrackV1 cTrack(*fkTrack);
-    cTrack.SetReconstructor(fReconstructor);
+    cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
     momentum = cTrack.GetMomentum(0);
     pdg = CalcPDG(&cTrack);
   }
@@ -800,13 +818,50 @@ TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
   } else {
     //AliWarning("No MC info available!");
     AliTRDtrackV1 cTrack(*fkTrack);
-    cTrack.SetReconstructor(fReconstructor);
+    cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
     momentum = cTrack.GetMomentum(0);
   }
   if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
   return hMomBin;
 }
 
+//_______________________________________________________
+TH1 *AliTRDcheckPID::PlotV0(const AliTRDtrackV1 *track)
+{
+  //
+  // Plot the V0 performance against MC
+  //
+
+  if(track) fkTrack = track;
+  if(!fkTrack){
+    AliDebug(2, "No Track defined.");
+    return NULL;
+  }  
+  if(!fkESD->HasV0()) return NULL;
+  if(!HasMCdata()){ 
+    AliDebug(1, "No MC defined.");
+    return NULL;
+  }
+  if(!fContainer){
+    AliWarning("No output container defined.");
+    return NULL;
+  }
+  AliDebug(2, Form("TRACK[%d] species[%s][%d]\n", fkESD->GetId(), fkMC->GetPID()>=0?AliPID::ParticleShortName(fkMC->GetPID()):"none", fkMC->GetPDG()));
+
+  TH1 *h(NULL);
+  if(!(h = dynamic_cast<TH1F*>(fContainer->At(kV0)))) return NULL;
+  Int_t sgn(0), n(0); AliTRDv0Info *v0(NULL);
+  for(Int_t iv0(fV0s->GetEntriesFast()); iv0--;){
+    if(!(v0=(AliTRDv0Info*)fV0s->At(iv0))) continue;
+    if(!(sgn = v0->HasTrack(fkESD->GetId()))) continue;
+    //for(Int_t is=AliPID::kSPECIES; is--;) v0->GetPID(is, track);
+    //v0->Print();
+    n++;
+    //break;
+  }
+  h->Fill(n);
+  return h;
+}
 
 //________________________________________________________
 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
@@ -1012,7 +1067,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
         h1->GetXaxis()->SetTitle("x_{drift} [cm]");
         h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
       }
-      h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
+      h1->DrawClone(kFIRST ? "c" : "samec");
       kFIRST = kFALSE;
     }
 
@@ -1045,7 +1100,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
         h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
         h1->GetYaxis()->SetTitle("Prob. [%]");
         h = (TH1F*)h1->DrawClone("c");
-        h->SetMaximum(55.);
+        h->SetMaximum(20.);
         h->GetXaxis()->SetRangeUser(0., 35.);
         kFIRST = kFALSE;
       } else h = (TH1F*)h1->DrawClone("samec");
@@ -1128,7 +1183,7 @@ Bool_t AliTRDcheckPID::PostProcess()
 }
 
 //________________________________________________________________________
-void AliTRDcheckPID::EvaluateEfficiency(TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
+void AliTRDcheckPID::EvaluateEfficiency(const TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
 // Process PID information for pion efficiency
 
   fUtil->SetElectronEfficiency(electronEfficiency);