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()
  56   ,fPDF(NULL)
  57 {
  58   //
  59   // AliTRDpidRefMakerLQ default constructor
  60   //
  61   SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker");
  62 }
  63
  64 //__________________________________________________________________
  65 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name)
  66   : AliTRDpidRefMaker(name, "PID(LQ) Reference Maker")
  67   ,fPDF(NULL)
  68 {
  69   //
  70   // AliTRDpidRefMakerLQ default constructor
  71   //
  72   DefineOutput(3, TObjArray::Class());
  73 }
  74
  75 //__________________________________________________________________
  76 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
  77 {
  78   //
  79   // AliTRDCalPIDQRef destructor
  80   //
  81   if(fPDF){
  82     //fPDF->Write("PDF_LQ", TObject::kSingleKey);
  83     fPDF->Delete();
  84     delete fPDF;
  85   }
  86 }
  87
  88 // //________________________________________________________________________
  89 void AliTRDpidRefMakerLQ::UserCreateOutputObjects()
  90 {
  91   // Create histograms
  92   // Called once
  93
  94   fContainer = Histos();
  95   PostData(1, fContainer);
  96
  97   //OpenFile(2, "RECREATE");
  98   fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
  99   fPDF->SetOwner();fPDF->SetName("PDF_LQ");
 100   PostData(3, fPDF);
 101 }
 102
 103
 104 //__________________________________________________________________
 105 TObjArray* AliTRDpidRefMakerLQ::Histos()
 106 {
 107   // Create histograms
 108
 109   if(fContainer) return fContainer;
 110   fContainer  = new TObjArray(AliTRDCalPID::kNMom);
 111
 112   // save dE/dx references
 113   TH1 *h(NULL);
 114   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
 115     TObjArray *arr = new TObjArray(2*AliPID::kSPECIES);
 116     arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
 117     for(Int_t is=AliPID::kSPECIES; is--;) {
 118       h = new TH1D(Form("h1%s%d", AliPID::ParticleShortName(is), ip), Form("1D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12.);
 119       h->GetXaxis()->SetTitle("log(dE/dx) [au]");
 120       h->GetYaxis()->SetTitle("#");
 121       h->SetLineColor(AliTRDCalPIDLQ::GetPartColor(is));
 122       arr->AddAt(h, is);
 123     }
 124     for(Int_t is=AliPID::kSPECIES; is--;) {
 125       h = new TH2D(Form("h2%s%d", AliPID::ParticleShortName(is), ip), Form("2D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
 126       h->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
 127       h->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
 128       h->GetZaxis()->SetTitle("#");
 129       arr->AddAt(h, AliPID::kSPECIES+is);
 130     }
 131     fContainer->AddAt(arr, ip);
 132   }
 133   fNRefFigures=AliTRDCalPID::kNMom;
 134   return fContainer;
 135 }
 136
 137
 138 //__________________________________________________________________
 139 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/)
 140 {
 141 // Steer loading of OCDB LQ PID
 142
 143   if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
 144     AliError(Form("File TRD.Calib%s.root not readable", GetName()));
 145     return NULL;
 146   }
 147   AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
 148   cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
 149   return cal;
 150 }
 151
 152 //__________________________________________________________________
 153 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
 154 {
 155 // Steering reference picture
 156
 157   if(ifig<0 || ifig>=fNRefFigures){
 158     AliError("Ref fig requested outside definition.");
 159     return kFALSE;
 160   }
 161   if(!gPad){
 162     AliWarning("Please provide a canvas to draw results.");
 163     return kFALSE;
 164   }
 165
 166   TObjArray *arr(NULL);TList *l(NULL);TH1 *h(NULL);
 167   if(!(arr = (TObjArray*)fContainer->At(ifig))){
 168     AliError(Form("PDF container @ pBin[%d] missing.", ifig));
 169     return kFALSE;
 170   }
 171   gPad->Divide(5, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
 172   for(Int_t is=0; is<AliPID::kSPECIES; is++){
 173     ((TVirtualPad*)l->At(is))->cd();
 174     if(!(h=(TH1*)arr->At(is))){
 175       AliError(Form("1D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
 176       return kFALSE;
 177     }
 178     h->GetYaxis()->SetRangeUser(0., 1.2);
 179     h->Draw("l");
 180
 181     ((TVirtualPad*)l->At(AliPID::kSPECIES+is))->cd();
 182     if(!(h=(TH1*)arr->At(AliPID::kSPECIES+is))){
 183       AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
 184       return kFALSE;
 185     }
 186     h->Draw("cont4z");
 187   }
 188
 189   return kTRUE;
 190 }
 191
 192
 193 //________________________________________________________________________
 194 Bool_t AliTRDpidRefMakerLQ::Load(const Char_t *file, const Char_t *dir)
 195 {
 196 // Load tree with data in case of detached PostProcess processing.
 197
 198   if(gSystem->AccessPathName(file, kReadPermission)){
 199     AliError(Form("File %s not readable", file));
 200     return kFALSE;
 201   }
 202   if(!TFile::Open(file)) {
 203     AliError(Form("File %s corrupted", file));
 204     return kFALSE;
 205   }
 206   if(!gFile->cd(dir)){
 207     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
 208     return kFALSE;
 209   }
 210   if (!(fData = dynamic_cast<TTree*>(gDirectory->Get("PIDrefMaker")))) {
 211     AliError("PIDref Tree not available");
 212     return kFALSE;
 213   }
 214   LinkPIDdata();
 215
 216   TObjArray *o(NULL);
 217 /*  if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
 218     AliWarning(Form("Monitor container Moni%s not available.", name));
 219     return kFALSE;
 220   }
 221   fContainer = (TObjArray*)o->Clone("monitor");
 222   fNRefFigures=AliTRDCalPID::kNMom;
 223 */
 224   // temporary until new calibration data are being produced
 225   Histos();
 226
 227
 228   if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
 229     AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
 230     return kFALSE;
 231   }
 232   if(!(o = (TObjArray*)gFile->Get(Form("PDF_LQ")))) {
 233     AliInfo("PDF container not available. Create.");
 234     fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
 235     fPDF->SetOwner();fPDF->SetName("PDF_LQ");
 236   } else fPDF = (TObjArray*)o->Clone("PDF_LQ");
 237
 238   return kTRUE;
 239 }
 240
 241 #include "AliAnalysisManager.h"
 242 //________________________________________________________________________
 243 void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/)
 244 {
 245 // Process pid info data array
 246 // The tasks have also access to the following containers which, for the time being, are not used
 247 // fTracks = dynamic_cast<TObjArray*>(GetInputData(1))
 248 // fEvent  = dynamic_cast<AliTRDeventInfo*>(GetInputData(2))
 249 // fV0s    = dynamic_cast<TObjArray*>(GetInputData(3))
 250
 251   if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(4)))){
 252     Int_t ev((Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry());
 253     AliDebug(3, Form("Missing pid info container in ev %d", ev));
 254     return;
 255   }
 256
 257   AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast()));
 258   AliTRDpidInfo *pid(NULL);
 259   const AliTRDpidInfo::AliTRDpidData *data(NULL);
 260   Char_t s(-1);
 261   for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
 262     if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
 263     if((s=pid->GetPID())<0) continue;
 264     for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
 265       data=pid->GetData(itrklt);
 266       Int_t ip(data->Momentum());
 267
 268       Double_t dedx[] = {0., 0.};
 269       if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
 270       ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
 271       AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
 272       ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]);
 273     }
 274   }
 275 }
 276
 277
 278 //________________________________________________________________________
 279 Bool_t AliTRDpidRefMakerLQ::PostProcess()
 280 {
 281 // Analyse merged dedx = f(p) distributions.
 282 //   - select momentum - species bins
 283 //   - rotate to principal components
 284 //   - locally interpolate with TKDPDF
 285 //   - save interpolation to monitoring histograms
 286 //   - write pdf to file for loading to OCDB
 287 //
 288
 289   TCanvas *fMonitor(NULL);
 290   // allocate working storage
 291   const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
 292   Float_t *data[2*kWS], *data1D[kWS];
 293   for(Int_t i(0); i<kWS; i++){
 294     data1D[i]   = new Float_t[kMaxStat];
 295     data[i]     = new Float_t[kMaxStat];
 296     data[kWS+i] = new Float_t[kMaxStat];
 297   }
 298   Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
 299
 300   AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries()));
 301   for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
 302     if(!(fData->GetEntry(itrk))) continue;
 303     Int_t sbin(fPIDdataArray->GetPID());
 304     for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
 305       Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum());
 306
 307       Double_t dedx[] = {0., 0.},
 308                dedx1D[] = {0., 0.};
 309       if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue;
 310       AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE);
 311       Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
 312       if(ndata[idx]==kMaxStat) continue;
 313
 314       // store data
 315       data1D[idx][ndata[idx]]   = dedx1D[0];
 316       data[idx][ndata[idx]]     = dedx[0];
 317       data[idx+kWS][ndata[idx]] = dedx[1];
 318       ndata[idx]++;
 319     }
 320   }
 321
 322   TKDPDF *pdf(NULL);
 323   Int_t in(0); Float_t par[6], *pp(NULL);
 324   for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
 325     for(Int_t is=AliPID::kSPECIES; is--;) {
 326       // estimate bucket statistics
 327       Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
 328             nb(kMinBuckets), // number of buckets
 329             ns((Int_t)(((Float_t)(ndata[idx]))/nb));    //statistics/bucket
 330
 331       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
 332       if(ns<Int_t(kMinStat)){
 333         AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
 334         continue;
 335       }
 336
 337       // build helper 1D PDF
 338       pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
 339       pdf->SetCOG(kFALSE);
 340       pdf->SetWeights();
 341       //pdf->SetStore();
 342       pdf->SetAlpha(15.);
 343       //pdf.GetStatus();
 344       fPDF->AddAt(pdf, idx);
 345       in=pdf->GetNTNodes(); pp=&par[0];
 346       while(in--){
 347         const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
 348         nn->GetCOG(pp);
 349         Double_t p(par[0]),r,e;
 350         pdf->Eval(&p,r,e,1);
 351       }
 352
 353       // build helper 2D PDF
 354       Float_t *ldata[2]={data[idx], data[kWS+idx]};
 355       pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
 356       pdf->SetCOG(kFALSE);
 357       pdf->SetWeights();
 358       //pdf->SetStore();
 359       pdf->SetAlpha(5.);
 360       //pdf.GetStatus();
 361       fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
 362       in=pdf->GetNTNodes(); pp=&par[0];
 363       while(in--){
 364         const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
 365         nn->GetCOG(pp);
 366         Double_t p[] = {par[0], par[1]}, r,e;
 367         pdf->Eval(p,r,e,1);
 368       }
 369 /*
 370       Int_t nnodes = pdf.GetNTNodes(),
 371             nside = Int_t(0.05*nnodes),
 372             nzeros = 4*(nside+1);
 373       printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
 374
 375
 376       // Build interpolator on the pdf skeleton
 377       TKDInterpolator interpolator(2, nnodes+nzeros);
 378       for(Int_t in=nnodes; in--;)
 379         interpolator.SetNode(in, *pdf.GetNodeInfo(in));
 380       TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
 381       Float_t ax0min, ax0max, ax1min, ax1max;
 382       pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
 383       pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
 384       printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
 385
 386       Int_t jn = nnodes;
 387       SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
 388       SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
 389       SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
 390       SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
 391       delete [] nodes;
 392       Int_t in=nnodes; Float_t par[6], *pp=&par[0];
 393       while(in--){
 394         const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
 395         nn->GetCOG(pp);
 396         //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
 397         Double_t p[] = {par[0], par[1]}, r,e;
 398         interpolator.Eval(p,r,e,1);
 399       }
 400 */
 401
 402       // visual on-line monitoring
 403       if(HasOnlineMonitor()){
 404         if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
 405         pdf->DrawProjection();
 406         fMonitor->Modified(); fMonitor->Update();
 407         fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
 408       }
 409
 410       //fContainer->ls();
 411       // save a discretization of the PDF for result monitoring
 412       Double_t rxy[]={0.,0.};
 413       TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is);
 414       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
 415       h2s->Clear();
 416       for(int ix=1; ix<=ax->GetNbins(); ix++){
 417         rxy[0] = ax->GetBinCenter(ix);
 418         for(int iy=1; iy<=ay->GetNbins(); iy++){
 419           rxy[1] = ay->GetBinCenter(iy);
 420
 421           Double_t rr,ee;
 422           pdf->Eval(rxy, rr, ee, kFALSE);
 423           if(rr<0. || ee/rr>.15) continue; // 15% relative error
 424           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
 425           h2s->SetBinContent(ix, iy, rr);
 426         }
 427       }
 428
 429       if(!(pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx)))){
 430         AliWarning(Form("Missing pdf for model id[%d]", idx));
 431         continue;
 432       }
 433       TH1 *h1 = (TH1D*)((TObjArray*)fContainer->At(ip))->At(is);
 434       ax = h1->GetXaxis();
 435       h1->Clear();
 436       for(int ix=1; ix<=ax->GetNbins(); ix++){
 437         rxy[0] = ax->GetBinCenter(ix);
 438
 439         Double_t rr,ee;
 440         pdf->Eval(rxy, rr, ee, kFALSE);
 441         if(rr<0. || ee/rr>.15) continue; // 15% relative error
 442         //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
 443         h1->SetBinContent(ix, rr);
 444       }
 445     }
 446   }
 447   for(Int_t i(0); i<kWS; i++){
 448     delete [] data1D[i];
 449     delete [] data[i];
 450     delete [] data[i+kWS];
 451   }
 452   return kTRUE;
 453 }
 454
 455
 456 //__________________________________________________________________
 457 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)
 458 {
 459 // Set extra nodes to ensure boundary conditions
 460
 461   printf("SetZeroes(%c)\n", opt);
 462   Float_t par[6], val[] = {0., 1.};
 463   Int_t a[6];
 464   if(opt=='x'){
 465     a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
 466   } else if(opt=='y'){
 467     a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
 468   } else return;
 469   Float_t tmp;
 470   par[a[1]] = y;
 471   par[a[4]] = y; par[a[5]] = y+dy;
 472   if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
 473   for(Int_t in=n; in--; node++, idx++, x+=dx){
 474     par[a[0]] = x+.5*dx;
 475     par[a[2]] = x;  par[a[3]] = x+dx;
 476     if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
 477     node->SetNode(2, par, val);
 478     printf("\n\tnode[%d]\n", idx); node->Print();
 479     interpolator->SetNode(idx, *node);
 480   }
 481   par[a[0]] = x;
 482   par[a[2]] = x;  par[a[3]] = x+dx;
 483   if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
 484   node->SetNode(2, par, val);
 485   printf("\n\tnode[%d]\n", idx); node->Print();
 486   interpolator->SetNode(idx, *node);node++;idx++;
 487 }
 488