]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/TRD/AliTRDcheckPID.cxx
fix error during merging (https://savannah.cern.ch/bugs/index.php?70729)
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckPID.cxx
index c58cb4d8239f5a89376fcfd786f25d446da711c2..c91948d0c9d58be652de802d15e59e7ef3e6ba57 100644 (file)
@@ -51,6 +51,8 @@
 #include "Cal/AliTRDCalPIDNN.h"
 #include "AliTRDcheckPID.h"
 #include "info/AliTRDtrackInfo.h"
+#include "info/AliTRDpidInfo.h"
+#include "info/AliTRDv0Info.h"
 
 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
 
@@ -58,18 +60,51 @@ ClassImp(AliTRDcheckPID)
 
 //________________________________________________________________________
 AliTRDcheckPID::AliTRDcheckPID() 
-  :AliTRDrecoTask("checkPID", "TRD PID checker")
-  ,fReconstructor(0x0)
-  ,fUtil(0x0)
-  ,fGraph(0x0)
-  ,fMomentumAxis(0x0)
+  :AliTRDrecoTask()
+  ,fReconstructor(NULL)
+  ,fUtil(NULL)
+  ,fGraph(NULL)
+  ,fPID(NULL)
+  ,fV0s(NULL)
+  ,fMomentumAxis(NULL)
   ,fMinNTracklets(AliTRDgeometry::kNlayer)
   ,fMaxNTracklets(AliTRDgeometry::kNlayer)
  {
   //
   // Default constructor
   //
+  SetNameTitle("checkPID", "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)
+ {
+  //
+  // Default constructor
+  //
+
+  LocalInit();
+  InitFunctorList();
+
+  DefineInput(2, TObjArray::Class());  // v0 list
+  DefineOutput(2, TObjArray::Class()); // pid info list
+}
+
 
+//________________________________________________________________________
+void AliTRDcheckPID::LocalInit() 
+{
+// Initialize working data
   fReconstructor = new AliTRDReconstructor();
   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 
@@ -82,27 +117,41 @@ AliTRDcheckPID::AliTRDcheckPID()
   memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
 
   fUtil = new AliTRDpidUtil();
-  InitFunctorList();
 }
 
-
 //________________________________________________________________________
 AliTRDcheckPID::~AliTRDcheckPID() 
 {
- if(fGraph){fGraph->Delete(); delete fGraph;}
- if(fReconstructor) delete fReconstructor;
- if(fUtil) delete fUtil;
+  if(fPID){fPID->Delete(); delete fPID;}
+  if(fGraph){fGraph->Delete(); delete fGraph;}
+  if(fReconstructor) delete fReconstructor;
+  if(fUtil) delete fUtil;
 }
 
 
 //________________________________________________________________________
-void AliTRDcheckPID::CreateOutputObjects()
+void AliTRDcheckPID::UserCreateOutputObjects()
 {
   // Create histograms
   // Called once
 
-  OpenFile(0, "RECREATE");
-  fContainer = Histos();
+  AliTRDrecoTask::UserCreateOutputObjects();
+  fPID = new TObjArray();
+  fPID->SetOwner(kTRUE);
+  PostData(2, fPID);
+}
+
+//________________________________________________________
+void AliTRDcheckPID::UserExec(Option_t *opt)
+{
+  //
+  // Execution part
+  //
+
+  fV0s = dynamic_cast<TObjArray *>(GetInputData(2));
+  fPID->Delete();
+
+  AliTRDrecoTask::UserExec(opt);
 }
 
 
@@ -118,7 +167,7 @@ TObjArray * AliTRDcheckPID::Histos(){
   fContainer = new TObjArray(); fContainer->Expand(kNPlots);
 
   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
-  TH1 *h = 0x0;
+  TH1 *h = NULL;
 
   const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
   // histos with posterior probabilities for all particle species
@@ -198,6 +247,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;
 }
 
@@ -252,13 +308,13 @@ TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
   //
 
   if(!fkESD){
-    AliWarning("No ESD info available.");
+    AliDebug(2, "No ESD info available.");
     return NULL;
   }
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
 
@@ -317,13 +373,13 @@ TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
   //
 
   if(!fkESD){
-    AliWarning("No ESD info available.");
+    AliDebug(2, "No ESD info available.");
     return NULL;
   }
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -382,13 +438,13 @@ TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
   //
 
   if(!fkESD){
-    AliWarning("No ESD info available.");
+    AliDebug(2, "No ESD info available.");
     return NULL;
   }
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -427,7 +483,7 @@ TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
       AliWarning("No Histogram defined.");
       return NULL;
     }
-  
+
     hPID->Fill(FindBin(species, momentum), pidESD[is]);
   }
   return hPID;
@@ -444,7 +500,7 @@ TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -468,25 +524,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;
 
+  // Init exchange container
+  Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
+  AliTRDpidInfo *pid = new AliTRDpidInfo(s);
+
   (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
-//   (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
-  Float_t sumdEdx = 0;
-  Int_t iBin = FindBin(species, momentum);
-  AliTRDseedV1 *tracklet = 0x0;
-  for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-    sumdEdx = 0;
-    tracklet = cTrack.GetTracklet(iChamb);
+
+  Float_t sumdEdx(0.);
+  Int_t iBin = FindBin(s, momentum);
+  AliTRDseedV1 *tracklet = NULL;
+  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;
 }
 
@@ -500,7 +564,7 @@ TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -530,7 +594,7 @@ TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
   Int_t iMomBin = fMomentumAxis->FindBin(momentum);
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   Float_t *fdEdx;
-  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDseedV1 *tracklet = NULL;
   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
     tracklet = cTrack.GetTracklet(iChamb);
     if(!tracklet) continue;
@@ -554,13 +618,13 @@ TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
   if(!CheckTrackQuality(fkTrack)) return NULL;
   
-  TObjArray *arr = 0x0;
+  TObjArray *arr = NULL;
   TProfile2D *hPHX, *hPHT;
   if(!(arr = dynamic_cast<TObjArray *>(fContainer->At(kPH)))){
     AliWarning("No Histogram defined.");
@@ -583,8 +647,8 @@ TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
   }
   if(!IsInRange(momentum)) return NULL;;
 
-  AliTRDseedV1 *tracklet = 0x0;
-  AliTRDcluster *cluster = 0x0;
+  AliTRDseedV1 *tracklet = NULL;
+  AliTRDcluster *cluster = NULL;
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   Int_t iBin = FindBin(species, momentum);
   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
@@ -610,7 +674,7 @@ TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -639,7 +703,7 @@ TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
 
   Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
   Int_t iBin = FindBin(species, momentum); 
-  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDseedV1 *tracklet = NULL;
   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
     tracklet = fkTrack->GetTracklet(iChamb);
     if(!tracklet) continue;
@@ -658,7 +722,7 @@ TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -677,7 +741,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);
@@ -697,7 +761,7 @@ TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -736,7 +800,7 @@ TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
 
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(2, "No Track defined.");
     return NULL;
   }
   
@@ -765,6 +829,42 @@ TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
   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(), AliPID::ParticleShortName(fkMC->GetPID()), fkMC->GetPDG()));
+
+  TH1 *h=dynamic_cast<TH1F*>(fContainer->At(kV0));
+  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)
@@ -772,12 +872,12 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
 // Steering function to retrieve performance plots
 
   Bool_t kFIRST = kTRUE;
-  TGraphErrors *g = 0x0;
-  TAxis *ax = 0x0;
-  TObjArray *arr = 0x0;
-  TH1 *h1 = 0x0, *h=0x0;
-  TH2 *h2 = 0x0;
-  TList *content = 0x0;
+  TGraphErrors *g = NULL;
+  TAxis *ax = NULL;
+  TObjArray *arr = NULL;
+  TH1 *h1 = NULL, *h=NULL;
+  TH2 *h2 = NULL;
+  TList *content = NULL;
   switch(ifig){
   case kEfficiency:{
     gPad->Divide(2, 1, 1.e-5, 1.e-5);
@@ -1003,7 +1103,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");