pidRefMaker for multidim LQ PID ready for mass production

- automatic tuning of bucket size/statistics
- monitoring plots
- control plots

diff --git a/TRD/qaRec/AliTRDpidRefMaker.cxx b/TRD/qaRec/AliTRDpidRefMaker.cxx
index 6517f4e391585205155e2923672fc90cc3d0ca31..c7657ab0379c5da6cec4fd73e591bd9732173a3c 100644 (file)
--- a/TRD/qaRec/AliTRDpidRefMaker.cxx
+++ b/TRD/qaRec/AliTRDpidRefMaker.cxx
@@ -3,7 +3,7 @@
 #include "TTree.h"
 #include "TEventList.h"
 #include "TObjArray.h"
-#include "TH1.h"
+#include "TH2.h"
 
 #include "AliLog.h"
 #include "AliESDtrack.h"
@@ -90,7 +90,16 @@ void AliTRDpidRefMaker::CreateOutputObjects()
   OpenFile(0, "RECREATE");
   fContainer = new TObjArray();
   fContainer->SetName(Form("Moni%s", GetName()));
-  fContainer->AddAt(new TH1F("hPDG","hPDG",AliPID::kSPECIES,-0.5,5.5),0);
+
+  TH2 *h2 = new TH2I("hPDG","Particle abundance", AliPID::kSPECIES, -0.5, 4.5, AliTRDCalPID::kNMom, -0.5, AliTRDCalPID::kNMom-0.5);
+  TAxis *ax = h2->GetXaxis();
+  ax->SetNdivisions(505);
+  ax->SetTitle("Particle species");
+  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);
 }
 
 //________________________________________________________________________

diff --git a/TRD/qaRec/AliTRDpidRefMakerLQ.cxx b/TRD/qaRec/AliTRDpidRefMakerLQ.cxx
index 2cbb7a75406f1e958c55b6668bd6a33252180291..d6f062d92f7e93b895a372f8d8762a907f730934 100644 (file)
--- a/TRD/qaRec/AliTRDpidRefMakerLQ.cxx
+++ b/TRD/qaRec/AliTRDpidRefMakerLQ.cxx
@@ -57,12 +57,11 @@ AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
   :AliTRDpidRefMaker("PidRefMakerLQ", "PID(LQ) Reference Maker")
   ,fPbin(-1)
   ,fSbin(-1)
+  ,fResults(0x0)
 {
   //
   // AliTRDpidRefMakerLQ default constructor
   //
-
-  memset(fH2dEdx, 0x0, AliPID::kSPECIES*sizeof(TH2*));
 }
 
 //__________________________________________________________________
@@ -71,10 +70,10 @@ AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
   //
   // AliTRDCalPIDQRef destructor
   //
-
-  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
-    if(fH2dEdx[ispec]) delete fH2dEdx[ispec];
-  }    
+  if(fResults) {
+    fResults->Delete();
+    delete fResults;
+  }
 }
 
 //________________________________________________________________________
@@ -85,6 +84,21 @@ void AliTRDpidRefMakerLQ::CreateOutputObjects()
 
   AliTRDpidRefMaker::CreateOutputObjects();
 
+  // save dE/dx references
+  TH2 *h2 = 0x0;
+  for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
+    TObjArray *arr = new TObjArray(AliPID::kSPECIES);
+    arr->SetName(Form("Pbin%02d", ip));
+    for(Int_t is=AliPID::kSPECIES; is--;) {
+      h2 = new TH2I(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
+      h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
+      h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
+      h2->GetZaxis()->SetTitle("#");
+      arr->AddAt(h2, is);
+    }
+    fContainer->AddAt(arr, 1+ip);
+  }
+
   fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
   fData->Branch("s", &fSbin, "l/b");
   fData->Branch("p", &fPbin, "p/b");
@@ -106,6 +120,40 @@ Float_t* AliTRDpidRefMakerLQ::GetdEdx(AliTRDseedV1 *trklt)
 }
 
 
+//__________________________________________________________________
+Bool_t AliTRDpidRefMakerLQ::GenerateOCDBEntry(Option_t *)
+{
+  return kTRUE;
+}
+
+//__________________________________________________________________
+Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
+{
+  if(ifig<0 || ifig>AliTRDCalPID::kNMom-1){ 
+    AliError("Ref fig requested outside definition.");
+    return kFALSE;
+  }
+  if(!fResults){ 
+    AliError("No results processed.");
+    return kFALSE;
+  }
+  if(!gPad){
+    AliWarning("Please provide a canvas to draw results.");
+    return kFALSE;
+  }
+
+  TObjArray *arr = (TObjArray*)fResults->At(ifig);
+  gPad->Divide(3, 2, 1.e-5, 1.e-5); 
+  TList *l=gPad->GetListOfPrimitives(); 
+  for(Int_t is=0; is<AliPID::kSPECIES; is++){
+    ((TVirtualPad*)l->At(is))->cd();
+    ((TH2*)arr->At(is))->Draw("colz");
+  }
+
+  return kTRUE;
+}
+
+
 //________________________________________________________________________
 void AliTRDpidRefMakerLQ::Fill()
 {
@@ -115,6 +163,11 @@ void AliTRDpidRefMakerLQ::Fill()
   fPbin = AliTRDpidUtil::GetMomentumBin(fP);
   // fill data
   fData->Fill();
+  // fill monitor
+  ((TH2*)fContainer->At(0))->Fill(fSbin, fPbin);
+  TH2* h2 = (TH2*)((TObjArray*)fContainer->At(1+fPbin))->At(fSbin);
+  h2->Fill(fdEdx[0], fdEdx[1]);
+  //printf("h[%s] : [%f] [%f]\n", h2->GetName(), fdEdx[0], fdEdx[1]);
 }
 
 //________________________________________________________________________
@@ -128,18 +181,34 @@ Bool_t AliTRDpidRefMakerLQ::PostProcess()
   }
   AliDebug(2, Form("Data[%d]", fData->GetEntries()));
 
-  Float_t *data[] = {0x0, 0x0};
-  TPrincipal principal(2, "ND");
+  // save PDF representation
+  TH2 *h2 = 0x0;
+  fResults = new TObjArray(AliTRDCalPID::kNMom);
+  for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
+    TObjArray *arr = new TObjArray(AliPID::kSPECIES);
+    arr->SetName(Form("Pbin%02d", ip));
+    for(Int_t is=AliPID::kSPECIES; is--;) {
+      h2 = new TH2I(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, -6., 6., 50, -6., 6.);
+      h2->GetXaxis()->SetTitle("log(dE/dx^{*}_{am}) [au]");
+      h2->GetYaxis()->SetTitle("log(dE/dx^{*}_{dr}) [au]");
+      h2->GetZaxis()->SetTitle("#");
+      arr->AddAt(h2, is);
+    }
+    fResults->AddAt(arr, ip);
+  }
+
 
   TCanvas *cc = new TCanvas("cc", "", 500, 500);
 
+  Float_t *data[] = {0x0, 0x0};
+  TPrincipal principal(2, "ND");
   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
     for(Int_t is=AliPID::kSPECIES; is--;) {
       principal.Clear();
       Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d]", ip, is, n));
-      if(n<10){
-        AliWarning(Form("Not enough data for %s[%d].", AliPID::ParticleShortName(is), ip));
+      if(n<1000/*Int_t(kMinStat)*Int_t(kMinBuckets)*/){
+        AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
         continue;
       }
       // allocate storage
@@ -165,27 +234,44 @@ Bool_t AliTRDpidRefMakerLQ::PostProcess()
       while((xx = principal.GetRow(++n))){
         principal.X2P(xx, rxy);
         data[0][n]=rxy[0]; data[1][n]=rxy[1];
-        //hProj->Fill(rxy[0], rxy[1]);
       }
 
-      // estimate acceptable statistical error per bucket
-
-      // estimate number of buckets
+      // estimate bucket statistics
+      Int_t ns(kMinStat),    //statistics/bucket
+            nb(kMinBuckets); // number of buckets
+      if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
+      else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
 
       // build PDF
-      TKDPDF pdf(n, 2, 10, data);
-      printf("PDF nodes[%d]\n", pdf.GetNTNodes());
+      TKDPDF pdf(n, 2, ns, data);
       pdf.SetStore();
       pdf.SetAlpha(5.);
+      //pdf.GetStatus();
       Float_t *c, v, ve; Double_t r, e;
       for(Int_t in=pdf.GetNTNodes(); in--;){
         pdf.GetCOGPoint(in, c, v, ve);
         rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
-        printf("%2d x[%f] y[%f]\n", in, rxy[0], rxy[1]);
         pdf.Eval(rxy, r, e, kTRUE);
       }
       pdf.DrawBins(0,1,-6,6,-6,6);cc->Modified(); cc->Update();
-      AliDebug(2, Form("n[%d]", n));
+
+
+      // save a discretization of the PDF for monitoring
+      TH2 *h2s = (TH2D*)((TObjArray*)fResults->At(ip))->At(is);
+      TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
+      for(int ix=1; ix<=ax->GetNbins(); ix++){
+        rxy[0] = ax->GetBinCenter(ix);
+        for(int iy=1; iy<=ay->GetNbins(); iy++){
+          rxy[1] = ay->GetBinCenter(iy);
+      
+          Double_t r,e;
+          pdf.Eval(rxy, r, e);
+          if(r<0. || e/r>.15) continue; // 15% relative error
+          //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, r, e);
+          h2s->SetBinContent(ix, iy, r);
+        }
+      }
+
 
       //pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
       
@@ -198,48 +284,3 @@ Bool_t AliTRDpidRefMakerLQ::PostProcess()
 }
 
 
-//__________________________________________________________________
-void   AliTRDpidRefMakerLQ::Reset()
-{
-  //
-  // Reset reference histograms
-  //
-
-  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
-    if(fH2dEdx[ispec]) fH2dEdx[ispec]->Reset();
-  }    
-}
-
-//__________________________________________________________________
-void  AliTRDpidRefMakerLQ::SaveReferences(const Int_t mom, const char *fn)
-{
-  //
-  // Save the reference histograms
-  //
-
-  TFile *fSave = 0x0;
-  TListIter it((TList*)gROOT->GetListOfFiles());
-  Bool_t kFOUND = kFALSE;
-  TDirectory *pwd = gDirectory;
-  while((fSave=(TFile*)it.Next()))
-    if(strcmp(fn, fSave->GetName())==0){
-      kFOUND = kTRUE;
-      break;
-    }
-  if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
-  fSave->cd();
-
-  // save dE/dx references
-  TH2 *h2 = 0x0;
-  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
-    h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom));
-    h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom));
-    h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
-    h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
-    h2->GetZaxis()->SetTitle("Entries");
-    h2->Write();
-  }
-
-  pwd->cd();
-}
-

diff --git a/TRD/qaRec/AliTRDpidRefMakerLQ.h b/TRD/qaRec/AliTRDpidRefMakerLQ.h
index 2cfc679c52ce09e8345d3fc3f2ff35520125cd5e..a3f7f2d68e45a543e5facbda96eb33be486f176d 100644 (file)
--- a/TRD/qaRec/AliTRDpidRefMakerLQ.h
+++ b/TRD/qaRec/AliTRDpidRefMakerLQ.h
@@ -23,16 +23,19 @@
 #include "AliTRDpidUtil.h"
 #endif
 
-class TH2;
+class TObjArray;
 class AliTRDpidRefMakerLQ : public AliTRDpidRefMaker {
 public:
   enum ETRDpidRefMakerLQsteer{
-    kMaxStat = 100000 // maximum statistics/species/momentum 
+    kMinStat = 100    // minimum statistics/bucket 10%
+   ,kMinBuckets = 450 // minimum number of buckets [lambda(6)*alpha(1.5)*regions(50)]
   };
   AliTRDpidRefMakerLQ();
   ~AliTRDpidRefMakerLQ();
  
   void      CreateOutputObjects();
+  Bool_t    GenerateOCDBEntry(Option_t *);
+  Bool_t    GetRefFigure(Int_t ifig);
   Bool_t    PostProcess();
 
 protected:
@@ -43,13 +46,11 @@ protected:
 private:
   AliTRDpidRefMakerLQ(const AliTRDpidRefMakerLQ &ref);
   AliTRDpidRefMakerLQ& operator=(const AliTRDpidRefMakerLQ &ref);
-  void   Reset();
-  void   SaveReferences(const Int_t mom, const char *fn);
  
 private:
   UChar_t   fPbin;        //! momentum bin
   UChar_t   fSbin;        //! species bin
-  TH2       *fH2dEdx[5];  //! dE/dx data holders
+  TObjArray *fResults;    //! array to store PDF representation
 
   ClassDef(AliTRDpidRefMakerLQ, 4)  // Reference builder for Multidim-LQ TRD-PID