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