PWG1/TRD/AliTRDpidRefMakerLQ.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 //
  21 //  TRD calibration class for building reference data for PID
  22 //  - 2D reference histograms (responsible A.Bercuci)
  23 //  - 3D reference histograms (not yet implemented) (responsible A.Bercuci)
  24 //
  25 //   Origin
  26 //   Alex Bercuci  (A.Bercuci@gsi.de)
  27 //
  28 ///////////////////////////////////////////////////////////////////////////////
  29
  30 #include <TSystem.h>
  31 #include <TROOT.h>
  32 #include <TFile.h>
  33 #include <TTree.h>
  34 #include <TMath.h>
  35 #include <TEventList.h>
  36 #include <TH2D.h>
  37 #include <TH2I.h>
  38 #include <TCanvas.h>
  39
  40 #include "AliLog.h"
  41 #include "../STAT/TKDPDF.h"
  42 #include "../STAT/TKDInterpolator.h"
  43 #include "AliTRDpidRefMakerLQ.h"
  44 #include "Cal/AliTRDCalPID.h"
  45 #include "Cal/AliTRDCalPIDLQ.h"
  46 #include "AliTRDseedV1.h"
  47 #include "AliTRDcalibDB.h"
  48 #include "AliTRDgeometry.h"
  49 #include "info/AliTRDpidInfo.h"
  50
  51 ClassImp(AliTRDpidRefMakerLQ)
  52
  53 //__________________________________________________________________
  54 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
  55   : AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
  56   , fPDF(NULL)
  57 {
  58   //
  59   // AliTRDpidRefMakerLQ default constructor
  60   //
  61   DefineOutput(2, TObjArray::Class());
  62 }
  63
  64 //__________________________________________________________________
  65 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
  66 {
  67   //
  68   // AliTRDCalPIDQRef destructor
  69   //
  70   if(fPDF){
  71     fPDF->Write("PDF_2DLQ", TObject::kSingleKey);
  72     fPDF->Delete();
  73     delete fPDF;
  74   }
  75 }
  76
  77 // //________________________________________________________________________
  78 void AliTRDpidRefMakerLQ::CreateOutputObjects()
  79 {
  80   // Create histograms
  81   // Called once
  82
  83   //OpenFile(0, "RECREATE");
  84   fContainer = Histos();
  85
  86   OpenFile(2, "RECREATE");
  87   fPDF = new TObjArray(AliTRDCalPID::kNMom*AliPID::kSPECIES);
  88   fPDF->SetOwner();fPDF->SetName("PDF_2DLQ");
  89 }
  90
  91
  92 //__________________________________________________________________
  93 TObjArray* AliTRDpidRefMakerLQ::Histos()
  94 {
  95   // Create histograms
  96
  97   if(fContainer) return fContainer;
  98
  99   fContainer = new TObjArray(AliTRDCalPID::kNMom);
 100   //fContainer->SetOwner(kTRUE);
 101   //fContainer->SetName(Form("Moni%s", GetName()));
 102
 103   // save dE/dx references
 104   TH2 *h2 = 0x0;
 105   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
 106     TObjArray *arr = new TObjArray(AliPID::kSPECIES);
 107     arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
 108     for(Int_t is=AliPID::kSPECIES; is--;) {
 109       h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
 110       h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
 111       h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
 112       h2->GetZaxis()->SetTitle("#");
 113       arr->AddAt(h2, is);
 114     }
 115     fContainer->AddAt(arr, ip);
 116   }
 117   fNRefFigures=AliTRDCalPID::kNMom;
 118   return fContainer;
 119 }
 120
 121
 122 //__________________________________________________________________
 123 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/)
 124 {
 125 // Steer loading of OCDB LQ PID
 126
 127   if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
 128     AliError(Form("File TRD.Calib%s.root not readable", GetName()));
 129     return NULL;
 130   }
 131   AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
 132   cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
 133   return cal;
 134 }
 135
 136 //__________________________________________________________________
 137 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
 138 {
 139 // Steering reference picture
 140
 141   if(ifig<0 || ifig>=fNRefFigures){
 142     AliError("Ref fig requested outside definition.");
 143     return kFALSE;
 144   }
 145   if(!gPad){
 146     AliWarning("Please provide a canvas to draw results.");
 147     return kFALSE;
 148   }
 149
 150   TObjArray *arr(NULL);TList *l(NULL);TH2 *h2(NULL);
 151   switch(ifig){
 152   case 0: // PDG plot
 153     if(!(h2=(TH2*)fContainer->At(0))){
 154       AliError("Abundance Plot missing.");
 155       return kFALSE;
 156     }
 157     h2->Draw("cont4z");
 158     break;
 159   default:
 160     if(!(arr = (TObjArray*)fContainer->At(ifig))){
 161       AliError(Form("2D container @ pBin[%d] missing.", ifig-1));
 162       return kFALSE;
 163     }
 164     gPad->Divide(3, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
 165     for(Int_t is=0; is<AliPID::kSPECIES; is++){
 166       ((TVirtualPad*)l->At(is))->cd();
 167       if(!(h2=(TH2*)arr->At(is))){
 168         AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig-1));
 169         return kFALSE;
 170       }
 171       h2->Draw("cont4z");
 172     }
 173   }
 174   return kTRUE;
 175 }
 176
 177
 178 //________________________________________________________________________
 179 Bool_t AliTRDpidRefMakerLQ::Load(const Char_t */*fname*/)
 180 {
 181   const Char_t *name("PIDrefMaker");
 182   if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
 183     AliError(Form("File TRD.Calib%s.root not readable", name));
 184     return kFALSE;
 185   }
 186   if(!TFile::Open(Form("TRD.Calib%s.root", name))){
 187     AliError(Form("File TRD.Calib%s.root corrupted", name));
 188     return kFALSE;
 189   }
 190   if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
 191     AliError(Form("Tree %s not available", name));
 192     return kFALSE;
 193   }
 194   LinkPIDdata();
 195
 196   TObjArray *o(NULL);
 197   if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
 198     AliWarning(Form("Monitor container Moni%s not available.", name));
 199     return kFALSE;
 200   }
 201   fContainer = (TObjArray*)o->Clone("monitor");
 202   fNRefFigures=AliTRDCalPID::kNMom;
 203
 204
 205   if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
 206     AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
 207     return kFALSE;
 208   }
 209   if(!(o = (TObjArray*)gFile->Get(Form("PDF_2DLQ")))) {
 210     AliInfo("PDF container PDF_2DLQ not available. Create.");
 211     fPDF = new TObjArray(AliTRDCalPID::kNMom*AliPID::kSPECIES);
 212     fPDF->SetOwner();fPDF->SetName("PDF_2DLQ");
 213   } else fPDF = (TObjArray*)o->Clone("PDF_2DLQ");
 214
 215   return kTRUE;
 216 }
 217
 218
 219 //________________________________________________________________________
 220 void AliTRDpidRefMakerLQ::Exec(Option_t */*opt*/)
 221 {
 222 // Load PID data into local data storage
 223
 224   AliTRDpidInfo *pid(NULL);
 225   const AliTRDpidInfo::AliTRDpidData *data(NULL);
 226   Char_t s(-1);
 227   for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
 228     if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
 229     if((s=pid->GetPID())<0) continue;
 230     for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
 231       data=pid->GetData(itrklt);
 232       Int_t ip(data->Momentum());
 233
 234
 235       Double_t dedx[] = {0., 0.};
 236       if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
 237       ((TH2*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0], dedx[1]);
 238     }
 239   }
 240
 241   PostData(0, fContainer);
 242   PostData(2, fPDF);
 243 }
 244
 245
 246 //________________________________________________________________________
 247 Bool_t AliTRDpidRefMakerLQ::PostProcess()
 248 {
 249 // Analyse merged dedx = f(p) distributions.
 250 //   - select momentum - species bins
 251 //   - rotate to principal components
 252 //   - locally interpolate with TKDPDF
 253 //   - save interpolation to monitoring histograms
 254 //   - write pdf to file for loading to OCDB
 255 //
 256
 257   TCanvas *fMonitor(NULL);
 258   // allocate working storage
 259   const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
 260   Float_t *data[2*kWS];
 261   for(Int_t i(0); i<2*kWS; i++) data[i]=new Float_t[kMaxStat];
 262   Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
 263
 264   AliDebug(1, Form("Loading data[%d]", fData->GetEntries()));
 265   for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
 266     if(!(fData->GetEntry(itrk))) continue;
 267     Int_t sbin(fPIDbin);
 268     for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
 269       Int_t pbin(fPIDdataArray->fData[ily].fPLbin & 0xf);
 270
 271       Double_t dedx[] = {0., 0.};
 272       if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->fData[ily].fdEdx, dedx)) continue;
 273       Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
 274       if(ndata[idx]==kMaxStat) continue;
 275
 276       // store data
 277       data[idx][ndata[idx]] = dedx[0];
 278       data[idx+kWS][ndata[idx]] = dedx[1];
 279       ndata[idx]++;
 280     }
 281   }
 282   for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
 283     for(Int_t is=AliPID::kSPECIES; is--;) {
 284       // estimate bucket statistics
 285       Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
 286             nb(kMinBuckets), // number of buckets
 287         ns((Int_t)(((Float_t)(ndata[idx]))/nb));    //statistics/bucket
 288
 289       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
 290       if(ns<Int_t(kMinStat)){
 291         AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
 292         continue;
 293       }
 294
 295       // build helper PDF
 296       Float_t *ldata[2]={data[idx], data[kWS+idx]};
 297       TKDPDF *pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
 298       pdf->SetCOG(kFALSE);
 299       pdf->SetWeights();
 300       //pdf->SetStore();
 301       pdf->SetAlpha(5.);
 302       //pdf.GetStatus();
 303       fPDF->AddAt(pdf, idx);
 304       Int_t in=pdf->GetNTNodes(); Float_t par[6], *pp=&par[0];
 305       while(in--){
 306         const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
 307         nn->GetCOG(pp);
 308         Double_t p[] = {par[0], par[1]}, r,e;
 309         pdf->Eval(p,r,e,1);
 310       }
 311 /*
 312       Int_t nnodes = pdf.GetNTNodes(),
 313             nside = Int_t(0.05*nnodes),
 314             nzeros = 4*(nside+1);
 315       printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
 316
 317
 318       // Build interpolator on the pdf skeleton
 319       TKDInterpolator interpolator(2, nnodes+nzeros);
 320       for(Int_t in=nnodes; in--;)
 321         interpolator.SetNode(in, *pdf.GetNodeInfo(in));
 322       TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
 323       Float_t ax0min, ax0max, ax1min, ax1max;
 324       pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
 325       pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
 326       printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
 327
 328       Int_t jn = nnodes;
 329       SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
 330       SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
 331       SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
 332       SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
 333       delete [] nodes;
 334       Int_t in=nnodes; Float_t par[6], *pp=&par[0];
 335       while(in--){
 336         const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
 337         nn->GetCOG(pp);
 338         //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
 339         Double_t p[] = {par[0], par[1]}, r,e;
 340         interpolator.Eval(p,r,e,1);
 341       }
 342 */
 343       // visual on-line monitoring
 344       if(HasOnlineMonitor()){
 345         if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
 346         pdf->DrawProjection();
 347         fMonitor->Modified(); fMonitor->Update();
 348         fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
 349       }
 350
 351       //fContainer->ls();
 352       // save a discretization of the PDF for result monitoring
 353       Double_t rxy[]={0.,0.};
 354       TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(is);
 355       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
 356       h2s->Clear();
 357       for(int ix=1; ix<=ax->GetNbins(); ix++){
 358         rxy[0] = ax->GetBinCenter(ix);
 359         for(int iy=1; iy<=ay->GetNbins(); iy++){
 360           rxy[1] = ay->GetBinCenter(iy);
 361
 362           Double_t rr,ee;
 363           pdf->Eval(rxy, rr, ee, kFALSE);
 364           if(rr<0. || ee/rr>.15) continue; // 15% relative error
 365           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
 366           h2s->SetBinContent(ix, iy, rr);
 367         }
 368       }
 369     }
 370   }
 371   for(Int_t i(0); i<2*kWS; i++) delete [] data[i];
 372   return kTRUE;
 373 }
 374
 375
 376 //__________________________________________________________________
 377 void AliTRDpidRefMakerLQ::SetZeroes(TKDInterpolator *interpolator, TKDNodeInfo *node, Int_t n, Int_t& idx, Float_t x, Float_t dx, Float_t y, Float_t dy, const Char_t opt)
 378 {
 379 // Set extra nodes to ensure boundary conditions
 380
 381   printf("SetZeroes(%c)\n", opt);
 382   Float_t par[6], val[] = {0., 1.};
 383   Int_t a[6];
 384   if(opt=='x'){
 385     a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
 386   } else if(opt=='y'){
 387     a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
 388   } else return;
 389   Float_t tmp;
 390   par[a[1]] = y;
 391   par[a[4]] = y; par[a[5]] = y+dy;
 392   if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
 393   for(Int_t in=n; in--; node++, idx++, x+=dx){
 394     par[a[0]] = x+.5*dx;
 395     par[a[2]] = x;  par[a[3]] = x+dx;
 396     if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
 397     node->SetNode(2, par, val);
 398     printf("\n\tnode[%d]\n", idx); node->Print();
 399     interpolator->SetNode(idx, *node);
 400   }
 401   par[a[0]] = x;
 402   par[a[2]] = x;  par[a[3]] = x+dx;
 403   if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
 404   node->SetNode(2, par, val);
 405   printf("\n\tnode[%d]\n", idx); node->Print();
 406   interpolator->SetNode(idx, *node);node++;idx++;
 407 }
 408