TRD/Cal/AliTRDCalPIDRefMaker.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$ */
  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 //  - Neural Network (responsible A.Wilk)
  25 //
  26 //   Origin
  27 //   Alex Bercuci  (A.Bercuci@gsi.de)
  28 //
  29 ///////////////////////////////////////////////////////////////////////////////
  30
  31 #include <TROOT.h>
  32 #include <TFile.h>
  33 #include <TTree.h>
  34 #include <TH2D.h>
  35 #include <TH2I.h>
  36 #include <TH3D.h>
  37 #include <TParticle.h>
  38 #include <TParticle.h>
  39 #include <TPrincipal.h>
  40 #include <TVector3.h>
  41 #include <TLinearFitter.h>
  42 #include <TVectorT.h>
  43 #include <TCanvas.h>
  44 #include <TEllipse.h>
  45 #include <TMarker.h>
  46
  47 #include "AliLog.h"
  48 #include "AliPID.h"
  49 #include "AliESD.h"
  50 #include "AliRun.h"
  51 #include "AliRunLoader.h"
  52 #include "AliStack.h"
  53 #include "AliHeader.h"
  54 #include "AliGenEventHeader.h"
  55 #include "AliESDtrack.h"
  56
  57 #include "AliTRDCalPIDRefMaker.h"
  58 #include "AliTRDCalPID.h"
  59 #include "AliTRDcalibDB.h"
  60 #include "AliTRDgeometry.h"
  61 #include "AliTRDtrack.h"
  62
  63 #include <vector>
  64
  65 ClassImp(AliTRDCalPIDRefMaker)
  66
  67 TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D2 = 0x0;
  68 TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D1 = 0x0;
  69
  70 //__________________________________________________________________
  71 AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker()
  72   :TObject()
  73 {
  74   //
  75   // AliTRDCalPIDRefMaker default constructor
  76   //
  77
  78         // histogram settings
  79         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
  80                 fH2dEdx[ispec] = 0x0;
  81                 fPrinc[ispec] = 0x0;
  82         }
  83 }
  84
  85 //__________________________________________________________________
  86 AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker(const AliTRDCalPIDRefMaker &ref)
  87   :TObject()
  88 {
  89   //
  90   // AliTRDCalPIDRefMaker copy constructor
  91   //
  92
  93         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
  94                 if(ref.fH2dEdx[ispec]){
  95                        fH2dEdx[ispec] = new  TH2D((TH2D&)*(ref.fH2dEdx[ispec]));
  96                 } else fH2dEdx[ispec] = 0x0;
  97                 fPrinc[ispec] = 0x0;
  98         }
  99 }
 100
 101 //__________________________________________________________________
 102 AliTRDCalPIDRefMaker::~AliTRDCalPIDRefMaker()
 103 {
 104   //
 105   // AliTRDCalPIDQRef destructor
 106   //
 107
 108         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 109                 if(fH2dEdx[ispec]) delete fH2dEdx[ispec];
 110                 if(fPrinc[ispec]) delete fPrinc[ispec];
 111         }
 112         if(fFitter2D1){ delete fFitter2D1; fFitter2D1 = 0x0;}
 113         if(fFitter2D2){ delete fFitter2D2; fFitter2D2 = 0x0;}
 114
 115 }
 116
 117 //__________________________________________________________________
 118 AliTRDCalPIDRefMaker& AliTRDCalPIDRefMaker::operator=(const AliTRDCalPIDRefMaker &ref)
 119 {
 120   //
 121   // Assignment operator
 122   //
 123
 124         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 125                 if(ref.fH2dEdx[ispec]) (ref.fH2dEdx[ispec])->Copy(*fH2dEdx[ispec]);
 126                 fPrinc[ispec] = 0x0;
 127         }
 128         return *this;
 129 }
 130
 131
 132 //__________________________________________________________________
 133 Bool_t AliTRDCalPIDRefMaker::BuildLQReferences(Char_t *File, Char_t *dir)
 134 {
 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 trackMomentum).
 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
 144         //
 145         // Author:
 146         // Alex Bercuci (A.Bercuci@gsi.de)
 147
 148   Int_t partCode[AliPID::kSPECIES] =
 149     {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
 150
 151         // check and retrive number of directories in the production
 152         Int_t nBatches;
 153         if(!(nBatches = CheckProdDirTree(dir))) return kFALSE;
 154
 155
 156         // Number of Time bins
 157         Int_t nTimeBins;
 158         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 159         if(!calibration){
 160           AliError("No AliTRDcalibDB available.");
 161           return kFALSE;
 162         } else {
 163                 nTimeBins = calibration->GetNumberOfTimeBins();
 164 //              if (calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
 165 //              else {
 166 //                      AliError("No run number set.");
 167 //          return kFALSE;
 168 //        }
 169         }
 170
 171         // Build PID reference/working objects
 172         fH1TB[0] = new TH1F("h1TB_0", "", nTimeBins, -.5, nTimeBins-.5);
 173         fH1TB[0]->SetLineColor(4);
 174         fH1TB[1] = new TH1F("h1TB_1", "", nTimeBins, -.5, nTimeBins-.5);
 175         fH1TB[1]->SetLineColor(2);
 176         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) if(!fPrinc[ispec]) fPrinc[ispec] = new TPrincipal(2, "ND");
 177
 178         // Print statistics header
 179         Int_t nPart[AliPID::kSPECIES], nTotPart;
 180         printf("P[GeV/c] ");
 181         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::GetPartSymb(ispec));
 182         printf("\n-----------------------------------------------\n");
 183
 184         Float_t trackMomentum[AliTRDCalPID::kNMom];
 185         //Float_t trackSegLength[AliTRDCalPID::kNLength];
 186         for(int i=0; i<AliTRDCalPID::kNMom; i++) trackMomentum[i] = AliTRDCalPID::GetMomentum(i);
 187         AliRunLoader *fRunLoader = 0x0;
 188         TFile *esdFile = 0x0;
 189         TTree *esdTree = 0x0;
 190         AliESD *esd = 0x0;
 191         AliESDtrack *esdTrack = 0x0;
 192
 193         //
 194         // Momentum loop
 195         for (Int_t imom = 0; imom < AliTRDCalPID::kNMom; imom++) {
 196                 Reset();
 197
 198                 // init statistics for momentum
 199                 nTotPart = 0;
 200                 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
 201
 202                 // loop over production directories
 203                 for(Int_t ibatch = 0; ibatch<nBatches; ibatch++){
 204                         // open run loader and load gAlice, kinematics and header
 205                         fRunLoader = AliRunLoader::Open(Form("%s/%3.1fGeV/%03d/galice.root", dir, trackMomentum[imom], ibatch));
 206                         if (!fRunLoader) {
 207                                 AliError(Form("Getting run loader for momentum %3.1f GeV/c batch %03d failed.", trackMomentum[imom], ibatch));
 208                                 return kFALSE;
 209                         }
 210                         TString s; s.Form("%s/%3.1fGeV/%03d/", dir, trackMomentum[imom], ibatch);
 211                         fRunLoader->SetDirName(s);
 212                         fRunLoader->LoadgAlice();
 213                         gAlice = fRunLoader->GetAliRun();
 214                         if (!gAlice) {
 215                                 AliError(Form("galice object not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
 216                                 return kFALSE;
 217                         }
 218                         fRunLoader->LoadKinematics();
 219                         fRunLoader->LoadHeader();
 220
 221                         // open the ESD file
 222                         esdFile = TFile::Open(Form("%s/%3.1fGeV/%03d/AliESDs.root", dir, trackMomentum[imom], ibatch));
 223                         if (!esdFile || esdFile->IsZombie()) {
 224                                 AliError(Form("Opening ESD file failed for momentum  %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
 225                                 return kFALSE;
 226                         }
 227                         esdTree = (TTree*)esdFile->Get("esdTree");
 228                         if (!esdTree) {
 229                                 AliError(Form("ESD tree not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch));
 230                                 return kFALSE;
 231                         }
 232                         esd = new AliESD;
 233                         esdTree->SetBranchAddress("ESD", &esd);
 234
 235
 236                         // Event loop
 237                         for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 238
 239                                 // Load MC info
 240                                 fRunLoader->GetEvent(iEvent);
 241                                 AliStack* stack = gAlice->Stack();
 242                                 TArrayF vertex(3);
 243                                 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 244
 245                                 // Load event summary data
 246                                 esdTree->GetEvent(iEvent);
 247                                 if (!esd) {
 248                                         AliWarning(Form("ESD object not found for event %d. (@ momentum %3.1f GeV/c, batch %03d)", iEvent, trackMomentum[imom], ibatch));
 249                                         continue;
 250                                 }
 251
 252                                 // Track loop
 253                                 for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
 254                                         esdTrack = esd->GetTrack(iTrack);
 255
 256                                         //if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
 257
 258                                         if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
 259                                         if(esdTrack->GetConstrainedChi2() > 1E9) continue;
 260                                         if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
 261                                         if (esdTrack->GetTRDsignal() == 0.) continue;
 262
 263                                         // read MC info
 264                                         Int_t label = esdTrack->GetLabel();
 265                                         if(label<0) continue;
 266                                         if (label > stack->GetNtrack()) continue;     // background
 267                                         TParticle* particle = stack->Particle(label);
 268                                         if(!particle){
 269                                                 AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ momentum %3.1f batch %03d event %d track %d]", label, trackMomentum[imom], ibatch, iEvent, iTrack));
 270                                                 continue;
 271                                         }
 272                                         if(particle->Pt() < 1.E-3) continue;
 273                                         //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
 274                                         TVector3 dVertex(particle->Vx() - vertex[0],
 275                                                                                 particle->Vy() - vertex[1],
 276                                                                                 particle->Vz() - vertex[2]);
 277                                         if (dVertex.Mag() > 1.E-4){
 278                                                 //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
 279                                                 //particle->Print();
 280                                                 continue;
 281                                         }
 282                                         Int_t iGen = -1;
 283                                         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
 284                                                 if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
 285                                                         iGen = ispec;
 286                                                         break;
 287                                                 }
 288                                         if(iGen<0) continue;
 289
 290                                         nPart[iGen]++; nTotPart++;
 291
 292                                         Float_t mom;
 293                                         //Float_t length;
 294                                         Double_t dedx[AliTRDtrack::kNslice], dEdx;
 295                                         Int_t timebin;
 296                                         for (Int_t iLayer=0; iLayer<AliTRDgeometry::kNlayer; iLayer++){
 297                                                 // read data for track segment
 298                                                 for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
 299                                                         dedx[iSlice] = esdTrack->GetTRDslice(iLayer, iSlice);
 300                                                 dEdx    = esdTrack->GetTRDslice(iLayer, -1);
 301                                                 timebin = esdTrack->GetTRDTimBin(iLayer);
 302
 303                                                 // check data
 304                                                 if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 305
 306                                                 // retrive kinematic info for this track segment
 307                                                 //if(!AliTRDpidESD::RecalculateTrackSegmentKine(esdTrack, iLayer, mom, length)) continue;
 308                                                 mom = esdTrack->GetOuterParam()->GetP();
 309
 310                                                 // find segment length and momentum bin
 311                                                 Int_t jmom = 1, refMom = -1;
 312                                                 while(jmom<AliTRDCalPID::kNMom-1 && mom>trackMomentum[jmom]) jmom++;
 313                                                 if(TMath::Abs(trackMomentum[jmom-1] - mom) < trackMomentum[jmom-1] * .2) refMom = jmom-1;
 314                                                 else if(TMath::Abs(trackMomentum[jmom] - mom) < trackMomentum[jmom] * .2) refMom = jmom;
 315                                                 if(refMom<0){
 316                                                         AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ momentum %3.1f batch %03d event %d track %d]", iLayer, trackMomentum[imom], ibatch, iEvent, iTrack));
 317                                                         continue;
 318                                                 }
 319                                                 /*while(jleng<AliTRDCalPID::kNLength-1 && length>trackSegLength[jleng]) jleng++;*/
 320
 321                                                 // this track segment has fulfilled all requierments
 322                                                 //nPlanePID++;
 323
 324                                                 if(dedx[0] > 0. && dedx[1] > 0.){
 325                                                         dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
 326                                                         fPrinc[iGen]->AddRow(dedx);
 327                                                 }
 328                                                 fH1TB[(iGen>0)?1:0]->Fill(timebin);
 329                                         } // end plane loop
 330                                 } // end track loop
 331                         } // end events loop
 332
 333                         delete esd; esd = 0x0;
 334                         esdFile->Close();
 335                         delete esdFile; esdFile = 0x0;
 336
 337                         fRunLoader->UnloadHeader();
 338                         fRunLoader->UnloadKinematics();
 339                         delete fRunLoader; fRunLoader = 0x0;
 340                 } // end directory loop
 341
 342                 // use data to prepare references
 343                 Prepare2D();
 344                 // save references
 345                 SaveReferences(imom, File);
 346
 347
 348                 // print momentum statistics
 349                 printf("  %3.1f  ", trackMomentum[imom]);
 350                 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
 351                 printf("\n");
 352         } // end momentum loop
 353
 354         TFile *fSave = 0x0;
 355         TListIter it((TList*)gROOT->GetListOfFiles());
 356         while((fSave=(TFile*)it.Next()))
 357                 if(strcmp(File, fSave->GetName())==0) break;
 358
 359         fSave->cd();
 360         fSave->Close();
 361         delete fSave;
 362
 363         return kTRUE;
 364 }
 365
 366 //__________________________________________________________________
 367 Bool_t AliTRDCalPIDRefMaker::BuildNNReferences(Char_t* /*File*/, Char_t* /*dir*/)
 368 {
 369         return kTRUE;
 370 }
 371
 372 //__________________________________________________________________
 373 Double_t AliTRDCalPIDRefMaker::Estimate2D2(TH2 *h, Float_t &x, Float_t &y)
 374 {
 375   //
 376   // Linear interpolation of data point with a parabolic expresion using
 377   // the logarithm of the bin content in a close neighbourhood. It is
 378   // assumed that the bin content of h is in number of events !
 379   //
 380   // Observation:
 381   // This function have to be used when the statistics of each bin is
 382   // sufficient and all bins are populated. For cases of sparse data
 383   // please refere to Estimate2D1().
 384   //
 385   // Author : Alex Bercuci (A.Bercuci@gsi.de)
 386   //
 387
 388         if(!h){
 389                 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D2()", "No histogram defined.");
 390                 return 0.;
 391         }
 392
 393         TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
 394         Int_t binx   = ax->FindBin(x);
 395         Int_t biny   = ay->FindBin(y);
 396         Int_t nbinsx = ax->GetNbins();
 397         Int_t nbinsy = ay->GetNbins();
 398         Double_t p[2];
 399         Double_t entries;
 400
 401         // late construction of fitter
 402         if(!fFitter2D2) fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
 403
 404         fFitter2D2->ClearPoints();
 405         Int_t npoints=0;
 406         Int_t binx0, binx1, biny0, biny1;
 407         for(int bin=0; bin<5; bin++){
 408                 binx0 = TMath::Max(1, binx-bin);
 409                 binx1 = TMath::Min(nbinsx, binx+bin);
 410                 for(int ibin=binx0; ibin<=binx1; ibin++){
 411                         biny0 = TMath::Max(1, biny-bin);
 412                         biny1 = TMath::Min(nbinsy, biny+bin);
 413                         for(int jbin=biny0; jbin<=biny1; jbin++){
 414                                 if(ibin != binx0 && ibin != binx1 && jbin != biny0 && jbin != biny1) continue;
 415                                 if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
 416                                 p[0] = ax->GetBinCenter(ibin);
 417                                 p[1] = ay->GetBinCenter(jbin);
 418                                 fFitter2D2->AddPoint(p, log(entries), 1./sqrt(entries));
 419                                 npoints++;
 420                         }
 421                 }
 422                 if(npoints>=25) break;
 423         }
 424         if(fFitter2D2->Eval() == 1){
 425                 printf("<I2> x = %9.4f y = %9.4f\n", x, y);
 426                 printf("\tbinx %d biny %d\n", binx, biny);
 427                 printf("\tpoints %d\n", npoints);
 428
 429                 return 0.;
 430         }
 431         TVectorD vec(6);
 432         fFitter2D2->GetParameters(vec);
 433         Double_t result = vec[0] + x*vec[1] + y*vec[2] + x*x*vec[3] + y*y*vec[4] + x*y*vec[5];
 434         return exp(result);
 435
 436 }
 437
 438 //__________________________________________________________________
 439 Double_t AliTRDCalPIDRefMaker::Estimate2D1(TH2 *h, Float_t &x, Float_t &y
 440                                       , Float_t &dCT, Float_t &rmin
 441                                       , Float_t &rmax)
 442 {
 443   //
 444   // Linear interpolation of data point with a plane using
 445   // the logarithm of the bin content in the area defined by the
 446   // d(cos(phi)) and dr=(rmin, rmax). It is assumed that the bin content
 447   // of h is number of events !
 448   //
 449
 450         if(!h){
 451                 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D1()", "No histogram defined.");
 452                 return 0.;
 453         }
 454
 455         TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis();
 456 //      Int_t binx   = ax->FindBin(x);
 457 //      Int_t biny   = ay->FindBin(y);
 458         Int_t nbinsx = ax->GetNbins();
 459         Int_t nbinsy = ay->GetNbins();
 460         Double_t p[2];
 461         Double_t entries;
 462         Double_t rxy = sqrt(x*x + y*y), rpxy;
 463
 464         // late construction of fitter
 465         if(!fFitter2D1) fFitter2D1 = new TLinearFitter(3, "1++x++y");
 466
 467         fFitter2D1->ClearPoints();
 468         Int_t npoints=0;
 469         for(int ibin=1; ibin<=nbinsx; ibin++){
 470                 for(int jbin=1; jbin<=nbinsy; jbin++){
 471                         if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue;
 472                         p[0] = ax->GetBinCenter(ibin);
 473                         p[1] = ay->GetBinCenter(jbin);
 474                         rpxy = sqrt(p[0]*p[0] + p[1]*p[1]);
 475                         if((x*p[0] + y*p[1])/rxy/rpxy < dCT) continue;
 476                         if(rpxy<rmin || rpxy > rmax) continue;
 477
 478                         fFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries));
 479                         npoints++;
 480                 }
 481         }
 482         if(npoints<15) return 0.;
 483         if(fFitter2D1->Eval() == 1){
 484                 printf("<O2> x = %9.4f y = %9.4f\n", x, y);
 485                 printf("\tpoints %d\n", npoints);
 486                 return 0.;
 487         }
 488         TVectorD vec(3);
 489         fFitter2D1->GetParameters(vec);
 490         Double_t result = vec[0] + x*vec[1] + y*vec[2];
 491         return exp(result);
 492 }
 493
 494 //__________________________________________________________________
 495 // Double_t     AliTRDCalPIDRefMaker::Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z)
 496 // {
 497 //      // Author Alex Bercuci (A.Bercuci@gsi.de)
 498 //      return 0.;
 499 // }
 500
 501
 502
 503 /////////////  Private functions ///////////////////////////////////
 504
 505 //__________________________________________________________________
 506 Int_t AliTRDCalPIDRefMaker::CheckProdDirTree(Char_t *dir)
 507 {
 508         // Scan directory tree for momenta. Returns the smallest number of
 509         // batches found in all directories or 0 if one momentum is missing.
 510
 511         const char *pwd = gSystem->pwd();
 512
 513         if(!gSystem->ChangeDirectory(dir)){
 514                 AliError(Form("Couldn't access production root directory %s.", dir));
 515                 return 0;
 516         }
 517
 518         Int_t iDir;
 519         Int_t nDir = Int_t(1.E6);
 520         for(int imom=0; imom<AliTRDCalPID::kNMom; imom++){
 521                 if(!gSystem->ChangeDirectory(Form("%3.1fGeV", AliTRDCalPID::GetMomentum(imom)))){
 522                         AliError(Form("Couldn't find data for momentum %3.1f GeV/c.", AliTRDCalPID::GetMomentum(imom)));
 523                         return 0;
 524                 }
 525
 526                 iDir = 0;
 527                 while(gSystem->ChangeDirectory(Form("%03d", iDir))){
 528                         iDir++;
 529                         gSystem->ChangeDirectory("..");
 530                 }
 531                 if(iDir < nDir) nDir = iDir;
 532                 gSystem->ChangeDirectory(dir);
 533         }
 534
 535         gSystem->ChangeDirectory(pwd);
 536
 537         return nDir;
 538 }
 539
 540
 541 //__________________________________________________________________
 542 void  AliTRDCalPIDRefMaker::Prepare2D()
 543 {
 544   //
 545   // Prepares the 2-dimensional reference histograms
 546   //
 547
 548         // histogram settings
 549         Float_t xmin, xmax, ymin, ymax;
 550         Int_t nbinsx, nbinsy, nBinsSector, nSectors;
 551         const Int_t color[] = {4, 3, 2, 7, 6};
 552         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 553                 // check PCA data
 554                 if(!fPrinc[ispec]){
 555                         AliError(Form("No data defined for %s.", AliTRDCalPID::GetPartName(ispec)));
 556                         return;
 557                 }
 558                 // build reference histograms
 559                 nBinsSector = 10;
 560                 nSectors = 20;
 561                 nbinsx = nbinsy = nBinsSector * nSectors;
 562                 xmin = ymin = 0.;
 563                 xmax = 8000.; ymax = 6000.;
 564                 if(!fH2dEdx[ispec]){
 565                         fH2dEdx[ispec] = new  TH2D(Form("h2%s", AliTRDCalPID::GetPartSymb(ispec)), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
 566                         fH2dEdx[ispec]->SetLineColor(color[ispec]);
 567                 }
 568         }
 569
 570         // build transformed rotated histograms
 571         nbinsx = nbinsy = 100;
 572         xmin = ymin = -6.;
 573         xmax = 10.; ymax = 6.;
 574         TH2I *hProj = 0x0;
 575         if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new  TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
 576         else hProj->Reset();
 577
 578         printf("Doing interpolation and invertion ... "); fflush(stdout);
 579         Bool_t kDebugPlot = kTRUE;
 580         //Bool_t kDebugPrint = kFALSE;
 581
 582         TCanvas *c = 0x0;
 583         TEllipse *ellipse = 0x0;
 584         TMarker *mark = 0x0;
 585         if(kDebugPlot){
 586                 c=new TCanvas("c2", "Interpolation 2D", 10, 10, 500, 500);
 587                 ellipse = new TEllipse();
 588                 ellipse->SetFillStyle(0); ellipse->SetLineColor(2);
 589                 mark = new TMarker();
 590                 mark->SetMarkerColor(2); mark->SetMarkerSize(2); mark->SetMarkerStyle(2);
 591         }
 592
 593         // define observable variables
 594         TAxis *ax, *ay;
 595         Double_t xy[2], lxy[2], rxy[2];
 596         Double_t estimate, position;
 597         const TVectorD *eValues;
 598
 599         // define radial sectors
 600         const Int_t   nPhi      = 36;
 601         const Float_t dPhi      = TMath::TwoPi()/nPhi;
 602         //const Float_t dPhiRange = .1;
 603         Int_t nPoints[nPhi], nFitPoints, binStart, binStop;
 604         TLinearFitter refsFitter[nPhi], refsLongFitter(6, "1++x++y++x*x++y*y++x*y");
 605         Float_t fFitterRange[nPhi];
 606         Bool_t kFitterStatus[nPhi];
 607         for(int iphi=0; iphi<nPhi; iphi++){
 608                 refsFitter[iphi].SetDim(3);
 609                 refsFitter[iphi].SetFormula("1++x++y");//++x*x++y*y++x*y");
 610                 fFitterRange[iphi] = .8;
 611                 kFitterStatus[iphi] = kFALSE;
 612         }
 613         std::vector<UShort_t> storeX[nPhi], storeY[nPhi];
 614
 615         // define working variables
 616         Float_t x0, y0, rx, ry;
 617         //Float_t rc, rmin, rmax, dr, dCT;
 618         Double_t Phi, r;
 619         Int_t iPhi;
 620         Double_t entries;
 621         for(int ispec=0; ispec<5; ispec++){
 622                 hProj->Reset();
 623                 for(int iphi=0; iphi<nPhi; iphi++){
 624                         nPoints[iphi] = 0;
 625                         refsFitter[iphi].ClearPoints();
 626                         fFitterRange[iphi] = .8;
 627                         kFitterStatus[iphi] = kFALSE;
 628                         storeX[iphi].clear();
 629                         storeY[iphi].clear();
 630                 }
 631
 632                 // calculate covariance ellipse
 633                 fPrinc[ispec]->MakePrincipals();
 634                 eValues  = fPrinc[ispec]->GetEigenValues();
 635                 x0  = 0.;
 636                 y0  = 0.;
 637                 rx  = 3.5*sqrt((*eValues)[0]);
 638                 ry  = 3.5*sqrt((*eValues)[1]);
 639
 640                 // rotate to principal axis
 641                 Int_t irow = 0;
 642                 const Double_t *xx;
 643                 printf("filling for spec %d ...\n", ispec);
 644                 while((xx = fPrinc[ispec]->GetRow(irow++))){
 645                         fPrinc[ispec]->X2P(xx, rxy);
 646                         hProj->Fill(rxy[0], rxy[1]);
 647                 }
 648
 649
 650 //      // debug plot
 651 //              if(kDebugPlot){
 652 //                      hProj->Draw();
 653 //                      ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.);
 654 //                      mark->DrawMarker(x0, y0);
 655 //                      gPad->Modified(); gPad->Update();
 656 //              }
 657
 658                 // define radial sectors
 659                 ax=hProj->GetXaxis();
 660                 ay=hProj->GetYaxis();
 661                 for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
 662                         rxy[0] = ax->GetBinCenter(ibin);
 663                         for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
 664                                 rxy[1] = ay->GetBinCenter(jbin);
 665
 666                                 if((entries = hProj->GetBinContent(ibin, jbin)) == 0) continue;
 667
 668                                 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
 669                                 if(position < 1.) continue;
 670
 671                                 r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
 672                                 Phi   = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
 673                                 iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1);
 674
 675                                 refsFitter[iPhi].AddPoint(rxy, log(entries), 1./sqrt(entries));
 676                                 nPoints[iPhi]++;
 677                         }
 678                 }
 679
 680                 // do interpolation
 681                 ax=fH2dEdx[ispec]->GetXaxis();
 682                 ay=fH2dEdx[ispec]->GetYaxis();
 683                 for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
 684                         xy[0]  = ax->GetBinCenter(ibin);
 685                         lxy[0] = log(xy[0]);
 686                         for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
 687                                 xy[1]  = ay->GetBinCenter(jbin);
 688                                 lxy[1] = log(xy[1]);
 689
 690                                 // rotate to PCA
 691                                 fPrinc[ispec]->X2P(lxy, rxy);
 692
 693                                 // calculate border of covariance ellipse
 694                                 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
 695
 696                                 // interpolation inside the covariance ellipse
 697                                 if(position < 1.){
 698                                   Float_t xTemp = rxy[0];
 699                                   Float_t yTemp = rxy[1];
 700                                         estimate = Estimate2D2((TH2 *) hProj, xTemp, yTemp);
 701                                         rxy[0] = xTemp;
 702                                         rxy[1] = yTemp;
 703                                         //                                      estimate = Estimate2D2((TH2 *) hProj, (Float_t)rxy[0], (Float_t)rxy[1]);
 704                                         fH2dEdx[ispec]->SetBinContent(ibin, jbin, estimate/xy[0]/xy[1]);
 705                                 } else { // interpolation outside the covariance ellipse
 706                                         r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
 707                                         Phi   = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
 708                                         iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1);
 709
 710                                         storeX[iPhi].push_back(ibin);
 711                                         storeY[iPhi].push_back(jbin);
 712                                 }
 713                         }
 714                 }
 715
 716                 // Fill outliers
 717                 // Radial fit on transformed rotated
 718                 TVectorD vec(3);
 719                 Int_t xbin, ybin;
 720                 for(int iphi=0; iphi<nPhi; iphi++){
 721                         Phi = iphi * dPhi - TMath::Pi();
 722                         if(TMath::Abs(TMath::Abs(Phi)-TMath::Pi()) < 100.*TMath::DegToRad()) continue;
 723
 724
 725                         refsFitter[iphi].Eval();
 726                         refsFitter[iphi].GetParameters(vec);
 727                         for(UInt_t ipoint=0; ipoint<storeX[iphi].size(); ipoint++){
 728                                 xbin = storeX[iphi].at(ipoint);
 729                                 ybin = storeY[iphi].at(ipoint);
 730                                 xy[0] = ax->GetBinCenter(xbin); lxy[0] = log(xy[0]);
 731                                 xy[1] = ay->GetBinCenter(ybin); lxy[1] = log(xy[1]);
 732
 733                                 // rotate to PCA
 734                                 fPrinc[ispec]->X2P(lxy, rxy);
 735                                 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]);
 736                                 fH2dEdx[ispec]->SetBinContent(xbin, ybin, estimate/xy[0]/xy[1]);
 737                         }
 738                 }
 739
 740                 // Longitudinal fit on ref histo in y direction
 741                 for(int isector=1; isector<nSectors; isector++){
 742                         // define sectors
 743                         binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
 744                         binStop  = (isector+1)*nBinsSector;
 745
 746                         nPoints[0] = 0;
 747                         refsLongFitter.ClearPoints();
 748                         storeX[0].clear();
 749                         storeY[0].clear();
 750
 751                         for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
 752                                 xy[0] = ax->GetBinCenter(ibin);
 753                                 for(int jbin=binStart; jbin<=binStop; jbin++){
 754                                         xy[1] = ay->GetBinCenter(jbin);
 755                                         if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
 756                                                 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
 757                                                 nPoints[0]++;
 758                                         } else {
 759                                                 storeX[0].push_back(ibin);
 760                                                 storeY[0].push_back(jbin);
 761                                         }
 762                                 }
 763                                 if(nPoints[0] >= ((isector==0)?60:420)) break;
 764                         }
 765
 766
 767                         if(refsLongFitter.Eval() == 1){
 768                                 printf("Error in Y sector %d\n", isector);
 769                                 continue;
 770                         }
 771                         refsLongFitter.GetParameters(vec);
 772                         for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){
 773                                 xbin = storeX[0].at(ipoint);
 774                                 ybin = storeY[0].at(ipoint);
 775                                 xy[0] = ax->GetBinCenter(xbin);
 776                                 xy[1] = ay->GetBinCenter(ybin);
 777
 778                                 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];
 779                                 fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
 780                         }
 781                 }
 782
 783                 // Longitudinal fit on ref histo in x direction
 784                 for(int isector=1; isector<nSectors; isector++){
 785                         binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
 786                         binStop  = (isector+1)*nBinsSector;
 787
 788                         nPoints[0] = 0;
 789                         nFitPoints = 0;
 790                         refsLongFitter.ClearPoints();
 791                         storeX[0].clear();
 792                         storeY[0].clear();
 793
 794                         for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
 795                                 xy[1] = ay->GetBinCenter(jbin);
 796                                 for(int ibin=binStart; ibin<=binStop; ibin++){
 797                                         xy[0] = ax->GetBinCenter(ibin);
 798                                         if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
 799                                                 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
 800                                                 nPoints[0]++;
 801                                         } else {
 802                                                 storeX[0].push_back(ibin);
 803                                                 storeY[0].push_back(jbin);
 804                                                 nFitPoints++;
 805                                         }
 806                                 }
 807                                 if(nPoints[0] >= ((isector==0)?60:420)) break;
 808                         }
 809                         if(nFitPoints == 0) continue;
 810
 811
 812                         if(refsLongFitter.Eval() == 1){
 813                                 printf("Error in X sector %d\n", isector);
 814                                 continue;
 815                         }
 816                         refsLongFitter.GetParameters(vec);
 817                         for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){
 818                                 xbin = storeX[0].at(ipoint);
 819                                 ybin = storeY[0].at(ipoint);
 820                                 xy[0] = ax->GetBinCenter(xbin);
 821                                 xy[1] = ay->GetBinCenter(ybin);
 822
 823                                 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];
 824                                 fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
 825                         }
 826                 }
 827
 828
 829                 fH2dEdx[ispec]->Scale(1./fH2dEdx[ispec]->Integral());
 830
 831                 // debug plot
 832                 if(kDebugPlot){
 833                         fH2dEdx[ispec]->Draw("cont1"); gPad->SetLogz();
 834                         gPad->Modified(); gPad->Update();
 835                 }
 836         }
 837         AliInfo("Finish interpolation.");
 838 }
 839
 840 //__________________________________________________________________
 841 void    AliTRDCalPIDRefMaker::Reset()
 842 {
 843   //
 844   // Reset reference histograms
 845   //
 846
 847         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 848                 if(fH2dEdx[ispec]) fH2dEdx[ispec]->Reset();
 849                 fPrinc[ispec]->Clear();
 850         }
 851         if(fH1TB[0]) fH1TB[0]->Reset();
 852         if(fH1TB[1]) fH1TB[1]->Reset();
 853 }
 854
 855 //__________________________________________________________________
 856 void  AliTRDCalPIDRefMaker::SaveReferences(const Int_t mom, const char *fn)
 857 {
 858   //
 859   // Save the reference histograms
 860   //
 861
 862         TFile *fSave = 0x0;
 863         TListIter it((TList*)gROOT->GetListOfFiles());
 864         Bool_t kFOUND = kFALSE;
 865         TDirectory *pwd = gDirectory;
 866         while((fSave=(TFile*)it.Next()))
 867                 if(strcmp(fn, fSave->GetName())==0){
 868                         kFOUND = kTRUE;
 869                         break;
 870                 }
 871         if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
 872         fSave->cd();
 873
 874         Float_t fmom = AliTRDCalPID::GetMomentum(mom);
 875
 876         // save dE/dx references
 877         TH2 *h2 = 0x0;
 878         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 879                 h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom));
 880                 h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom));
 881                 h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
 882                 h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
 883                 h2->GetZaxis()->SetTitle("Entries");
 884                 h2->Write();
 885         }
 886
 887         // save maximum time bin references
 888         TH1 *h1 = 0x0;
 889         h1 = (TH1F*)fH1TB[0]->Clone(Form("h1MaxTBEL%02d", mom));
 890         h1->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fmom));
 891         h1->GetXaxis()->SetTitle("time [100 ns]");
 892         h1->GetYaxis()->SetTitle("Probability");
 893         h1->Write();
 894
 895         h1 = (TH1F*)fH1TB[1]->Clone(Form("h1MaxTBPI%02d", mom));
 896         h1->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fmom));
 897         h1->GetXaxis()->SetTitle("time [100 ns]");
 898         h1->GetYaxis()->SetTitle("Probability");
 899         h1->Write();
 900
 901
 902         pwd->cd();
 903 }
 904