1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
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 // - Neural Network (responsible A.Wilk)
27 // Alex Bercuci (A.Bercuci@gsi.de)
29 ///////////////////////////////////////////////////////////////////////////////
37 #include <TParticle.h>
38 #include <TParticle.h>
39 #include <TPrincipal.h>
41 #include <TLinearFitter.h>
51 #include "AliRunLoader.h"
53 #include "AliHeader.h"
54 #include "AliGenEventHeader.h"
55 #include "AliESDtrack.h"
57 #include "AliTRDCalPIDRefMaker.h"
58 #include "AliTRDCalPIDLQ.h"
59 #include "AliTRDcalibDB.h"
60 #include "AliTRDgeometry.h"
61 #include "AliTRDtrack.h"
65 ClassImp(AliTRDCalPIDRefMaker)
67 TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D2 = 0x0;
68 TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D1 = 0x0;
70 //__________________________________________________________________
71 AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker()
75 // AliTRDCalPIDRefMaker default constructor
79 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
85 //__________________________________________________________________
86 AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker(const AliTRDCalPIDRefMaker &ref)
90 // AliTRDCalPIDRefMaker copy constructor
93 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
94 if(ref.h2dEdx[ispec]){
95 h2dEdx[ispec] = new TH2D((TH2D&)*(ref.h2dEdx[ispec]));
96 } else h2dEdx[ispec] = 0x0;
101 //__________________________________________________________________
102 AliTRDCalPIDRefMaker::~AliTRDCalPIDRefMaker()
105 // AliTRDCalPIDQRef destructor
108 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
109 if(h2dEdx[ispec]) delete h2dEdx[ispec];
110 if(fPrinc[ispec]) delete fPrinc[ispec];
112 if(fFitter2D1){ delete fFitter2D1; fFitter2D1 = 0x0;}
113 if(fFitter2D2){ delete fFitter2D2; fFitter2D2 = 0x0;}
117 //__________________________________________________________________
118 AliTRDCalPIDRefMaker& AliTRDCalPIDRefMaker::operator=(const AliTRDCalPIDRefMaker &ref)
121 // Assignment operator
124 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
125 if(ref.h2dEdx[ispec]) (ref.h2dEdx[ispec])->Copy(*h2dEdx[ispec]);
132 //__________________________________________________________________
133 Bool_t AliTRDCalPIDRefMaker::BuildLQReferences(Char_t *File, Char_t *dir)
135 // Build, Fill and write to file the histograms used for PID.
136 // The simulations are looked in the
137 // directories with the general form Form("p%3.1f", momentum)
138 // starting from dir (default .). Here momentum belongs to the list
139 // of known momentum to PID (see fTrackMomentum).
140 // The output histograms are
141 // written to the file "File" in cwd (default
142 // TRDPIDHistograms.root). In order to build a DB entry
143 // consider running $ALICE_ROOT/Cal/AliTRDpidCDB.C
146 // Alex Bercuci (A.Bercuci@gsi.de)
148 Int_t partCode[AliPID::kSPECIES] =
149 {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
151 // check and retrive number of directories in the production
153 if(!(nBatches = CheckProdDirTree(dir))) return kFALSE;
156 // Number of Time bins
158 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
160 AliError("No AliTRDcalibDB available.");
163 nTimeBins = calibration->GetNumberOfTimeBins();
164 // if (calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
166 // AliError("No run number set.");
171 // Build PID reference/working objects
172 h1TB[0] = new TH1F("h1TB_0", "", nTimeBins, -.5, nTimeBins-.5);
173 h1TB[0]->SetLineColor(4);
174 h1TB[1] = new TH1F("h1TB_1", "", nTimeBins, -.5, nTimeBins-.5);
175 h1TB[1]->SetLineColor(2);
176 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) if(!fPrinc[ispec]) fPrinc[ispec] = new TPrincipal(2, "ND");
178 // Print statistics header
179 Int_t nPart[AliPID::kSPECIES], nTotPart;
181 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPIDLQ::fpartSymb[ispec]);
182 printf("\n-----------------------------------------------\n");
184 Float_t fTrackMomentum[AliTRDCalPIDLQ::kNMom], fTrackSegLength[AliTRDCalPIDLQ::kNLength];
185 for(int i=0; i<AliTRDCalPIDLQ::kNMom; i++) fTrackMomentum[i] = AliTRDCalPIDLQ::GetMomentum(i);
186 AliRunLoader *fRunLoader = 0x0;
187 TFile *esdFile = 0x0;
188 TTree *esdTree = 0x0;
190 AliESDtrack *esdTrack = 0x0;
194 for (Int_t imom = 0; imom < AliTRDCalPIDLQ::kNMom; imom++) {
197 // init statistics for momentum
199 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
201 // loop over production directories
202 for(Int_t ibatch = 0; ibatch<nBatches; ibatch++){
203 // open run loader and load gAlice, kinematics and header
204 fRunLoader = AliRunLoader::Open(Form("%s/%3.1fGeV/%03d/galice.root", dir, fTrackMomentum[imom], ibatch));
206 AliError(Form("Getting run loader for momentum %3.1f GeV/c batch %03d failed.", fTrackMomentum[imom], ibatch));
209 TString s; s.Form("%s/%3.1fGeV/%03d/", dir, fTrackMomentum[imom], ibatch);
210 fRunLoader->SetDirName(s);
211 fRunLoader->LoadgAlice();
212 gAlice = fRunLoader->GetAliRun();
214 AliError(Form("galice object not found for momentum %3.1f GeV/c batch %03d.", fTrackMomentum[imom], ibatch));
217 fRunLoader->LoadKinematics();
218 fRunLoader->LoadHeader();
221 esdFile = TFile::Open(Form("%s/%3.1fGeV/%03d/AliESDs.root", dir, fTrackMomentum[imom], ibatch));
222 if (!esdFile || esdFile->IsZombie()) {
223 AliError(Form("Opening ESD file failed for momentum %3.1f GeV/c batch %03d.", fTrackMomentum[imom], ibatch));
226 esdTree = (TTree*)esdFile->Get("esdTree");
228 AliError(Form("ESD tree not found for momentum %3.1f GeV/c batch %03d.", fTrackMomentum[imom], ibatch));
232 esdTree->SetBranchAddress("ESD", &esd);
236 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
239 fRunLoader->GetEvent(iEvent);
240 AliStack* stack = gAlice->Stack();
242 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
244 // Load event summary data
245 esdTree->GetEvent(iEvent);
247 AliWarning(Form("ESD object not found for event %d. (@ momentum %3.1f GeV/c, batch %03d)", iEvent, fTrackMomentum[imom], ibatch));
252 for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
253 esdTrack = esd->GetTrack(iTrack);
255 //if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
257 if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
258 if(esdTrack->GetConstrainedChi2() > 1E9) continue;
259 if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
260 if (esdTrack->GetTRDsignal() == 0.) continue;
263 Int_t label = esdTrack->GetLabel();
264 if(label<0) continue;
265 if (label > stack->GetNtrack()) continue; // background
266 TParticle* particle = stack->Particle(label);
268 AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ momentum %3.1f batch %03d event %d track %d]", label, fTrackMomentum[imom], ibatch, iEvent, iTrack));
271 if(particle->Pt() < 1.E-3) continue;
272 // if (TMath::Abs(particle->Eta()) > 0.3) continue;
273 TVector3 dVertex(particle->Vx() - vertex[0],
274 particle->Vy() - vertex[1],
275 particle->Vz() - vertex[2]);
276 if (dVertex.Mag() > 1.E-4){
277 //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
282 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
283 if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
289 nPart[iGen]++; nTotPart++;
292 Double_t dedx[AliTRDtrack::kNslice], dEdx;
294 for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
295 // read data for track segment
296 for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
297 dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
298 dEdx = esdTrack->GetTRDsignals(iPlane, -1);
299 timebin = esdTrack->GetTRDTimBin(iPlane);
302 if ((dEdx <= 0.) || (timebin <= -1.)) continue;
304 // retrive kinematic info for this track segment
305 //if(!AliTRDpidESD::RecalculateTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
306 mom = esdTrack->GetOuterParam()->GetP();
308 // find segment length and momentum bin
309 Int_t jmom = 1, refMom = -1;
310 while(jmom<AliTRDCalPIDLQ::kNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
311 if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
312 else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
314 AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ momentum %3.1f batch %03d event %d track %d]", iPlane, fTrackMomentum[imom], ibatch, iEvent, iTrack));
317 /*while(jleng<AliTRDCalPIDLQ::kNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
319 // this track segment has fulfilled all requierments
322 if(dedx[0] > 0. && dedx[1] > 0.){
323 dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
324 fPrinc[iGen]->AddRow(dedx);
326 h1TB[(iGen>0)?1:0]->Fill(timebin);
331 delete esd; esd = 0x0;
333 delete esdFile; esdFile = 0x0;
335 fRunLoader->UnloadHeader();
336 fRunLoader->UnloadKinematics();
337 delete fRunLoader; fRunLoader = 0x0;
338 } // end directory loop
340 // use data to prepare references
343 SaveReferences(imom, File);
346 // print momentum statistics
347 printf(" %3.1f ", fTrackMomentum[imom]);
348 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
350 } // end momentum loop
353 TListIter it((TList*)gROOT->GetListOfFiles());
354 while((fSave=(TFile*)it.Next()))
355 if(strcmp(File, fSave->GetName())==0) break;
364 //__________________________________________________________________
365 Bool_t AliTRDCalPIDRefMaker::BuildNNReferences(Char_t *File, Char_t *dir)
370 //__________________________________________________________________
371 Double_t AliTRDCalPIDRefMaker::Estimate2D2(TH2 *h, Float_t &x, Float_t &y)
374 // Linear interpolation of data point with a parabolic expresion using
375 // the logarithm of the bin content in a close neighbourhood. It is
376 // assumed that the bin content of h is in number of events !
379 // This function have to be used when the statistics of each bin is
380 // sufficient and all bins are populated. For cases of sparse data
381 // please refere to Estimate2D1().
383 // Author : Alex Bercuci (A.Bercuci@gsi.de)
387 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D2()", "No histogram defined.");
391 TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
392 Int_t binx = ax->FindBin(x);
393 Int_t biny = ay->FindBin(y);
394 Int_t nbinsx = ax->GetNbins();
395 Int_t nbinsy = ay->GetNbins();
399 // late construction of fitter
400 if(!fFitter2D2) fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
402 fFitter2D2->ClearPoints();
404 Int_t binx0, binx1, biny0, biny1;
405 for(int bin=0; bin<5; bin++){
406 binx0 = TMath::Max(1, binx-bin);
407 binx1 = TMath::Min(nbinsx, binx+bin);
408 for(int ibin=binx0; ibin<=binx1; ibin++){
409 biny0 = TMath::Max(1, biny-bin);
410 biny1 = TMath::Min(nbinsy, biny+bin);
411 for(int jbin=biny0; jbin<=biny1; jbin++){
412 if(ibin != binx0 && ibin != binx1 && jbin != biny0 && jbin != biny1) continue;
413 if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
414 p[0] = ax->GetBinCenter(ibin);
415 p[1] = ay->GetBinCenter(jbin);
416 fFitter2D2->AddPoint(p, log(entries), 1./sqrt(entries));
420 if(npoints>=25) break;
422 if(fFitter2D2->Eval() == 1){
423 printf("<I2> x = %9.4f y = %9.4f\n", x, y);
424 printf("\tbinx %d biny %d\n", binx, biny);
425 printf("\tpoints %d\n", npoints);
430 fFitter2D2->GetParameters(vec);
431 Double_t result = vec[0] + x*vec[1] + y*vec[2] + x*x*vec[3] + y*y*vec[4] + x*y*vec[5];
436 //__________________________________________________________________
437 Double_t AliTRDCalPIDRefMaker::Estimate2D1(TH2 *h, Float_t &x, Float_t &y
438 , Float_t &dCT, Float_t &rmin
442 // Linear interpolation of data point with a plane using
443 // the logarithm of the bin content in the area defined by the
444 // d(cos(phi)) and dr=(rmin, rmax). It is assumed that the bin content
445 // of h is number of events !
449 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D1()", "No histogram defined.");
453 TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
454 // Int_t binx = ax->FindBin(x);
455 // Int_t biny = ay->FindBin(y);
456 Int_t nbinsx = ax->GetNbins();
457 Int_t nbinsy = ay->GetNbins();
460 Double_t rxy = sqrt(x*x + y*y), rpxy;
462 // late construction of fitter
463 if(!fFitter2D1) fFitter2D1 = new TLinearFitter(3, "1++x++y");
465 fFitter2D1->ClearPoints();
467 for(int ibin=1; ibin<=nbinsx; ibin++){
468 for(int jbin=1; jbin<=nbinsy; jbin++){
469 if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
470 p[0] = ax->GetBinCenter(ibin);
471 p[1] = ay->GetBinCenter(jbin);
472 rpxy = sqrt(p[0]*p[0] + p[1]*p[1]);
473 if((x*p[0] + y*p[1])/rxy/rpxy < dCT) continue;
474 if(rpxy<rmin || rpxy > rmax) continue;
476 fFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries));
480 if(npoints<15) return 0.;
481 if(fFitter2D1->Eval() == 1){
482 printf("<O2> x = %9.4f y = %9.4f\n", x, y);
483 printf("\tpoints %d\n", npoints);
487 fFitter2D1->GetParameters(vec);
488 Double_t result = vec[0] + x*vec[1] + y*vec[2];
492 //__________________________________________________________________
493 // Double_t AliTRDCalPIDRefMaker::Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z)
495 // // Author Alex Bercuci (A.Bercuci@gsi.de)
501 ///////////// Private functions ///////////////////////////////////
503 //__________________________________________________________________
504 Int_t AliTRDCalPIDRefMaker::CheckProdDirTree(Char_t *dir)
506 // Scan directory tree for momenta. Returns the smallest number of
507 // batches found in all directories or 0 if one momentum is missing.
509 const char *pwd = gSystem->pwd();
511 if(!gSystem->ChangeDirectory(dir)){
512 AliError(Form("Couldn't access production root directory %s.", dir));
517 Int_t nDir = Int_t(1.E6);
518 for(int imom=0; imom<AliTRDCalPIDLQ::kNMom; imom++){
519 if(!gSystem->ChangeDirectory(Form("%3.1fGeV", AliTRDCalPIDLQ::GetMomentum(imom)))){
520 AliError(Form("Couldn't find data for momentum %3.1f GeV/c.", AliTRDCalPIDLQ::GetMomentum(imom)));
525 while(gSystem->ChangeDirectory(Form("%03d", iDir))){
527 gSystem->ChangeDirectory("..");
529 if(iDir < nDir) nDir = iDir;
530 gSystem->ChangeDirectory(dir);
533 gSystem->ChangeDirectory(pwd);
539 //__________________________________________________________________
540 void AliTRDCalPIDRefMaker::Prepare2D()
543 // Prepares the 2-dimensional reference histograms
546 // histogram settings
547 Float_t xmin, xmax, ymin, ymax;
548 Int_t nbinsx, nbinsy, nBinsSector, nSectors;
549 const Int_t color[] = {4, 3, 2, 7, 6};
550 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
553 AliError(Form("No data defined for %s.", AliTRDCalPIDLQ::fpartName[ispec]));
556 // build reference histograms
559 nbinsx = nbinsy = nBinsSector * nSectors;
561 xmax = 8000.; ymax = 6000.;
563 h2dEdx[ispec] = new TH2D(Form("h2%s", AliTRDCalPIDLQ::fpartSymb[ispec]), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
564 h2dEdx[ispec]->SetLineColor(color[ispec]);
568 // build transformed rotated histograms
569 nbinsx = nbinsy = 100;
571 xmax = 10.; ymax = 6.;
573 if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
576 printf("Doing interpolation and invertion ... "); fflush(stdout);
577 Bool_t kDebugPlot = kTRUE;
578 //Bool_t kDebugPrint = kFALSE;
581 TEllipse *ellipse = 0x0;
584 c=new TCanvas("c2", "Interpolation 2D", 10, 10, 500, 500);
585 ellipse = new TEllipse();
586 ellipse->SetFillStyle(0); ellipse->SetLineColor(2);
587 mark = new TMarker();
588 mark->SetMarkerColor(2); mark->SetMarkerSize(2); mark->SetMarkerStyle(2);
591 // define observable variables
593 Double_t xy[2], lxy[2], rxy[2];
594 Double_t estimate, position;
595 const TVectorD *eValues;
597 // define radial sectors
598 const Int_t nPhi = 36;
599 const Float_t dPhi = TMath::TwoPi()/nPhi;
600 const Float_t dPhiRange = .1;
601 Int_t nPoints[nPhi], nFitPoints, binStart, binStop;
602 TLinearFitter refsFitter[nPhi], refsLongFitter(6, "1++x++y++x*x++y*y++x*y");
603 Float_t fFitterRange[nPhi];
604 Bool_t kFitterStatus[nPhi];
605 for(int iphi=0; iphi<nPhi; iphi++){
606 refsFitter[iphi].SetDim(3);
607 refsFitter[iphi].SetFormula("1++x++y");//++x*x++y*y++x*y");
608 fFitterRange[iphi] = .8;
609 kFitterStatus[iphi] = kFALSE;
611 std::vector<UShort_t> storeX[nPhi], storeY[nPhi];
613 // define working variables
614 Float_t x0, y0, rx, ry, rc, rmin, rmax, dr, dCT;
618 for(int ispec=0; ispec<5; ispec++){
620 for(int iphi=0; iphi<nPhi; iphi++){
622 refsFitter[iphi].ClearPoints();
623 fFitterRange[iphi] = .8;
624 kFitterStatus[iphi] = kFALSE;
625 storeX[iphi].clear();
626 storeY[iphi].clear();
629 // calculate covariance ellipse
630 fPrinc[ispec]->MakePrincipals();
631 eValues = fPrinc[ispec]->GetEigenValues();
634 rx = 3.5*sqrt((*eValues)[0]);
635 ry = 3.5*sqrt((*eValues)[1]);
637 // rotate to principal axis
640 printf("filling for spec %d ...\n", ispec);
641 while((xx = fPrinc[ispec]->GetRow(irow++))){
642 fPrinc[ispec]->X2P(xx, rxy);
643 hProj->Fill(rxy[0], rxy[1]);
650 // ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.);
651 // mark->DrawMarker(x0, y0);
652 // gPad->Modified(); gPad->Update();
655 // define radial sectors
656 ax=hProj->GetXaxis();
657 ay=hProj->GetYaxis();
658 for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
659 rxy[0] = ax->GetBinCenter(ibin);
660 for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
661 rxy[1] = ay->GetBinCenter(jbin);
663 if((entries = hProj->GetBinContent(ibin, jbin)) == 0) continue;
665 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
666 if(position < 1.) continue;
668 r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
669 Phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
670 iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1);
672 refsFitter[iPhi].AddPoint(rxy, log(entries), 1./sqrt(entries));
678 ax=h2dEdx[ispec]->GetXaxis();
679 ay=h2dEdx[ispec]->GetYaxis();
680 for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
681 xy[0] = ax->GetBinCenter(ibin);
683 for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
684 xy[1] = ay->GetBinCenter(jbin);
688 fPrinc[ispec]->X2P(lxy, rxy);
690 // calculate border of covariance ellipse
691 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
693 // interpolation inside the covariance ellipse
695 estimate = Estimate2D2(hProj, (Float_t)rxy[0], (Float_t)rxy[1]);
696 h2dEdx[ispec]->SetBinContent(ibin, jbin, estimate/xy[0]/xy[1]);
697 } else { // interpolation outside the covariance ellipse
698 r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
699 Phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
700 iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1);
702 storeX[iPhi].push_back(ibin);
703 storeY[iPhi].push_back(jbin);
709 // Radial fit on transformed rotated
712 for(int iphi=0; iphi<nPhi; iphi++){
713 Phi = iphi * dPhi - TMath::Pi();
714 if(fabs(fabs(Phi)-TMath::Pi()) < 100.*TMath::DegToRad()) continue;
717 refsFitter[iphi].Eval();
718 refsFitter[iphi].GetParameters(vec);
719 for(int ipoint=0; ipoint<storeX[iphi].size(); ipoint++){
720 xbin = storeX[iphi].at(ipoint);
721 ybin = storeY[iphi].at(ipoint);
722 xy[0] = ax->GetBinCenter(xbin); lxy[0] = log(xy[0]);
723 xy[1] = ay->GetBinCenter(ybin); lxy[1] = log(xy[1]);
726 fPrinc[ispec]->X2P(lxy, rxy);
727 estimate = exp(vec[0] + rxy[0]*vec[1] + rxy[1]*vec[2]);//+ rxy[0]*rxy[0]*vec[3]+ rxy[1]*rxy[1]*vec[4]+ rxy[0]*rxy[1]*vec[5]);
728 h2dEdx[ispec]->SetBinContent(xbin, ybin, estimate/xy[0]/xy[1]);
732 // Longitudinal fit on ref histo in y direction
733 for(int isector=1; isector<nSectors; isector++){
735 binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
736 binStop = (isector+1)*nBinsSector;
739 refsLongFitter.ClearPoints();
743 for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
744 xy[0] = ax->GetBinCenter(ibin);
745 for(int jbin=binStart; jbin<=binStop; jbin++){
746 xy[1] = ay->GetBinCenter(jbin);
747 if((entries = h2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
748 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
751 storeX[0].push_back(ibin);
752 storeY[0].push_back(jbin);
755 if(nPoints[0] >= ((isector==0)?60:420)) break;
759 if(refsLongFitter.Eval() == 1){
760 printf("Error in Y sector %d\n", isector);
763 refsLongFitter.GetParameters(vec);
764 for(int ipoint=0; ipoint<storeX[0].size(); ipoint++){
765 xbin = storeX[0].at(ipoint);
766 ybin = storeY[0].at(ipoint);
767 xy[0] = ax->GetBinCenter(xbin);
768 xy[1] = ay->GetBinCenter(ybin);
770 estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5];
771 h2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
775 // Longitudinal fit on ref histo in x direction
776 for(int isector=1; isector<nSectors; isector++){
777 binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
778 binStop = (isector+1)*nBinsSector;
782 refsLongFitter.ClearPoints();
786 for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
787 xy[1] = ay->GetBinCenter(jbin);
788 for(int ibin=binStart; ibin<=binStop; ibin++){
789 xy[0] = ax->GetBinCenter(ibin);
790 if((entries = h2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
791 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
794 storeX[0].push_back(ibin);
795 storeY[0].push_back(jbin);
799 if(nPoints[0] >= ((isector==0)?60:420)) break;
801 if(nFitPoints == 0) continue;
804 if(refsLongFitter.Eval() == 1){
805 printf("Error in X sector %d\n", isector);
808 refsLongFitter.GetParameters(vec);
809 for(int ipoint=0; ipoint<storeX[0].size(); ipoint++){
810 xbin = storeX[0].at(ipoint);
811 ybin = storeY[0].at(ipoint);
812 xy[0] = ax->GetBinCenter(xbin);
813 xy[1] = ay->GetBinCenter(ybin);
815 estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5];
816 h2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
821 h2dEdx[ispec]->Scale(1./h2dEdx[ispec]->Integral());
825 h2dEdx[ispec]->Draw("cont1"); gPad->SetLogz();
826 gPad->Modified(); gPad->Update();
829 AliInfo("Finish interpolation.");
832 //__________________________________________________________________
833 void AliTRDCalPIDRefMaker::Reset()
836 // Reset reference histograms
839 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
840 if(h2dEdx[ispec]) h2dEdx[ispec]->Reset();
841 fPrinc[ispec]->Clear();
843 if(h1TB[0]) h1TB[0]->Reset();
844 if(h1TB[1]) h1TB[1]->Reset();
847 //__________________________________________________________________
848 void AliTRDCalPIDRefMaker::SaveReferences(const Int_t mom, const char *fn)
851 // Save the reference histograms
855 TListIter it((TList*)gROOT->GetListOfFiles());
856 Bool_t kFOUND = kFALSE;
857 TDirectory *pwd = gDirectory;
858 while((fSave=(TFile*)it.Next()))
859 if(strcmp(fn, fSave->GetName())==0){
863 if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
866 Float_t fmom = AliTRDCalPIDLQ::GetMomentum(mom);
868 // save dE/dx references
870 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
871 h2 = (TH2D*)h2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPIDLQ::fpartSymb[ispec], mom));
872 h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPIDLQ::fpartName[ispec], mom));
873 h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
874 h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
875 h2->GetZaxis()->SetTitle("Entries");
879 // save maximum time bin references
881 h1 = (TH1F*)h1TB[0]->Clone(Form("h1MaxTBEL%02d", mom));
882 h1->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fmom));
883 h1->GetXaxis()->SetTitle("time [100 ns]");
884 h1->GetYaxis()->SetTitle("Probability");
887 h1 = (TH1F*)h1TB[1]->Clone(Form("h1MaxTBPI%02d", mom));
888 h1->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fmom));
889 h1->GetXaxis()->SetTitle("time [100 ns]");
890 h1->GetYaxis()->SetTitle("Probability");