]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/TRD/AliTRDcheckPID.cxx
Adding 2 triggers to list of known unsupported triggers
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckPID.cxx
index 93200e9619766455a191d2de0f2ae5d6c00fb041..3cbc966b951d67f386233adc14261672107c55d4 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,17 +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)
@@ -93,7 +95,7 @@ AliTRDcheckPID::AliTRDcheckPID(char* name )
   LocalInit();
   InitFunctorList();
 
-  DefineInput(2, TObjArray::Class());  // v0 list
+  DefineInput(3, TObjArray::Class());  // v0 list
   DefineOutput(2, TObjArray::Class()); // pid info list
 }
 
@@ -102,8 +104,6 @@ AliTRDcheckPID::AliTRDcheckPID(char* name )
 void AliTRDcheckPID::LocalInit() 
 {
 // Initialize working data
-  fReconstructor = new AliTRDReconstructor();
-  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 
   // Initialize momentum axis with default values
   Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
@@ -121,7 +121,6 @@ AliTRDcheckPID::~AliTRDcheckPID()
 {
   if(fPID){fPID->Delete(); delete fPID;}
   if(fGraph){fGraph->Delete(); delete fGraph;}
-  if(fReconstructor) delete fReconstructor;
   if(fUtil) delete fUtil;
 }
 
@@ -132,11 +131,10 @@ void AliTRDcheckPID::UserCreateOutputObjects()
   // Create histograms
   // Called once
 
-  if(!HasFunctorList()) InitFunctorList();
-  fContainer = Histos();
-
+  AliTRDrecoTask::UserCreateOutputObjects();
   fPID = new TObjArray();
   fPID->SetOwner(kTRUE);
+  PostData(2, fPID);
 }
 
 //________________________________________________________
@@ -146,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);
 }
 
 
@@ -163,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;
@@ -206,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);
 
@@ -246,6 +243,19 @@ 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);
+
+  // dQ/dl for 1D-Likelihood
+  if(!(h = (TH1F *)gROOT->FindObject("dQdl"))){
+    h = new TH2F("dQdl", "dQ/dl per layer;p*species;dQ/dl [a.u.]", xBins, -0.5, xBins - 0.5, 800, 0., 40000.);
+  } else h->Reset();
+  fContainer->AddAt(h, kdQdl);
+
   return fContainer;
 }
 
@@ -270,9 +280,9 @@ Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
 {
 // Documentation to come
 
 track -> SetReconstructor(fReconstructor);
-  (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
-  track -> CookPID();
/* 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)){
     return kElectron;
@@ -317,7 +327,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.;
@@ -332,7 +342,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);
@@ -382,7 +392,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.;
@@ -396,7 +406,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;
 
@@ -455,7 +465,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);
   }
@@ -476,8 +486,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;
@@ -485,6 +493,49 @@ TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
 
 
 
+//_______________________________________________________
+TH1 *AliTRDcheckPID::PlotdQdl(const AliTRDtrackV1 *track){
+  //
+  // Plot the total charge for the 1D Likelihood method
+  //
+  if(track) fkTrack = track;
+  if(!fkTrack){
+    AliDebug(2, "No Track defined");
+    return NULL;
+  }
+  TH2 *hdQdl = dynamic_cast<TH2F *>(fContainer->At(kdQdl));
+  if(!hdQdl){
+    AliWarning("No Histogram defined");
+    return NULL;
+  }
+
+  if(!CheckTrackQuality(fkTrack)) return NULL;
+
+  Int_t pdg = 0;
+  Float_t momentum = 0.;
+  AliTRDtrackV1 cTrack(*fkTrack);
+  if(fkMC){
+    if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
+    pdg = fkMC->GetPDG();
+  } else {
+    //AliWarning("No MC info available!");
+    momentum = cTrack.GetMomentum(0);
+    pdg = CalcPDG(&cTrack);
+  }
+  if(!IsInRange(momentum)) return NULL;
+
+  // Init exchange container
+  Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
+  Int_t ibin = FindBin(s, momentum);
+
+  AliTRDseedV1 *tracklet = NULL;
+  for(Int_t iseed = 0; iseed < 6; iseed++){
+    if(!((tracklet = fkTrack->GetTracklet(iseed)) && tracklet->IsOK())) continue;
+    hdQdl->Fill(ibin, tracklet->GetdQdl());
+  }
+  return hdQdl;
+}
+
 //_______________________________________________________
 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
 {
@@ -507,7 +558,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){
@@ -518,25 +569,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;
 }
 
@@ -563,7 +622,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){
@@ -576,7 +635,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;
@@ -627,7 +686,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);
   }
@@ -641,10 +700,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;
@@ -681,7 +740,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);
   }
@@ -719,7 +778,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){
@@ -727,7 +786,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);
@@ -768,7 +827,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);
   }
@@ -808,13 +867,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)
@@ -1020,7 +1116,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;
     }
 
@@ -1136,7 +1232,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);