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 "AliTRDCalPIDLQ.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                 h2dEdx[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.h2dEdx[ispec]){
  95                         h2dEdx[ispec] = new  TH2D((TH2D&)*(ref.h2dEdx[ispec]));
  96                 } else h2dEdx[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(h2dEdx[ispec]) delete h2dEdx[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.h2dEdx[ispec]) (ref.h2dEdx[ispec])->Copy(*h2dEdx[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 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
 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         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");
 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[%%] ", AliTRDCalPIDLQ::fpartSymb[ispec]);
 182         printf("\n-----------------------------------------------\n");
 183
 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;
 189         AliESD *esd = 0x0;
 190         AliESDtrack *esdTrack = 0x0;
 191
 192         //
 193         // Momentum loop
 194         for (Int_t imom = 0; imom < AliTRDCalPIDLQ::kNMom; imom++) {
 195                 Reset();
 196
 197                 // init statistics for momentum
 198                 nTotPart = 0;
 199                 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
 200
 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));
 205                         if (!fRunLoader) {
 206                                 AliError(Form("Getting run loader for momentum %3.1f GeV/c batch %03d failed.", fTrackMomentum[imom], ibatch));
 207                                 return kFALSE;
 208                         }
 209                         TString s; s.Form("%s/%3.1fGeV/%03d/", dir, fTrackMomentum[imom], ibatch);
 210                         fRunLoader->SetDirName(s);
 211                         fRunLoader->LoadgAlice();
 212                         gAlice = fRunLoader->GetAliRun();
 213                         if (!gAlice) {
 214                                 AliError(Form("galice object not found for momentum %3.1f GeV/c batch %03d.", fTrackMomentum[imom], ibatch));
 215                                 return kFALSE;
 216                         }
 217                         fRunLoader->LoadKinematics();
 218                         fRunLoader->LoadHeader();
 219
 220                         // open the ESD file
 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));
 224                                 return kFALSE;
 225                         }
 226                         esdTree = (TTree*)esdFile->Get("esdTree");
 227                         if (!esdTree) {
 228                                 AliError(Form("ESD tree not found for momentum %3.1f GeV/c batch %03d.", fTrackMomentum[imom], ibatch));
 229                                 return kFALSE;
 230                         }
 231                         esd = new AliESD;
 232                         esdTree->SetBranchAddress("ESD", &esd);
 233
 234
 235                         // Event loop
 236                         for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 237
 238                                 // Load MC info
 239                                 fRunLoader->GetEvent(iEvent);
 240                                 AliStack* stack = gAlice->Stack();
 241                                 TArrayF vertex(3);
 242                                 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 243
 244                                 // Load event summary data
 245                                 esdTree->GetEvent(iEvent);
 246                                 if (!esd) {
 247                                         AliWarning(Form("ESD object not found for event %d. (@ momentum %3.1f GeV/c, batch %03d)", iEvent, fTrackMomentum[imom], ibatch));
 248                                         continue;
 249                                 }
 250
 251                                 // Track loop
 252                                 for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
 253                                         esdTrack = esd->GetTrack(iTrack);
 254
 255                                         //if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
 256
 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;
 261
 262                                         // read MC info
 263                                         Int_t label = esdTrack->GetLabel();
 264                                         if(label<0) continue;
 265                                         if (label > stack->GetNtrack()) continue;     // background
 266                                         TParticle* particle = stack->Particle(label);
 267                                         if(!particle){
 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));
 269                                                 continue;
 270                                         }
 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));
 278                                                 //particle->Print();
 279                                                 continue;
 280                                         }
 281                                         Int_t iGen = -1;
 282                                         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
 283                                                 if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
 284                                                         iGen = ispec;
 285                                                         break;
 286                                                 }
 287                                         if(iGen<0) continue;
 288
 289                                         nPart[iGen]++; nTotPart++;
 290
 291                                         Float_t mom, length;
 292                                         Double_t dedx[AliTRDtrack::kNslice], dEdx;
 293                                         Int_t timebin;
 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);
 300
 301                                                 // check data
 302                                                 if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 303
 304                                                 // retrive kinematic info for this track segment
 305                                                 //if(!AliTRDpidESD::RecalculateTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
 306                                                 mom = esdTrack->GetOuterParam()->GetP();
 307
 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;
 313                                                 if(refMom<0){
 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));
 315                                                         continue;
 316                                                 }
 317                                                 /*while(jleng<AliTRDCalPIDLQ::kNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
 318
 319                                                 // this track segment has fulfilled all requierments
 320                                                 //nPlanePID++;
 321
 322                                                 if(dedx[0] > 0. && dedx[1] > 0.){
 323                                                         dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
 324                                                         fPrinc[iGen]->AddRow(dedx);
 325                                                 }
 326                                                 h1TB[(iGen>0)?1:0]->Fill(timebin);
 327                                         } // end plane loop
 328                                 } // end track loop
 329                         } // end events loop
 330
 331                         delete esd; esd = 0x0;
 332                         esdFile->Close();
 333                         delete esdFile; esdFile = 0x0;
 334
 335                         fRunLoader->UnloadHeader();
 336                         fRunLoader->UnloadKinematics();
 337                         delete fRunLoader; fRunLoader = 0x0;
 338                 } // end directory loop
 339
 340                 // use data to prepare references
 341                 Prepare2D();
 342                 // save references
 343                 SaveReferences(imom, File);
 344
 345
 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);
 349                 printf("\n");
 350         } // end momentum loop
 351
 352         TFile *fSave = 0x0;
 353         TListIter it((TList*)gROOT->GetListOfFiles());
 354         while((fSave=(TFile*)it.Next()))
 355                 if(strcmp(File, fSave->GetName())==0) break;
 356
 357         fSave->cd();
 358         fSave->Close();
 359         delete fSave;
 360
 361         return kTRUE;
 362 }
 363
 364 //__________________________________________________________________
 365 Bool_t AliTRDCalPIDRefMaker::BuildNNReferences(Char_t *File, Char_t *dir)
 366 {
 367         return kTRUE;
 368 }
 369
 370 //__________________________________________________________________
 371 Double_t AliTRDCalPIDRefMaker::Estimate2D2(TH2 *h, Float_t &x, Float_t &y)
 372 {
 373   //
 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 !
 377   //
 378   // Observation:
 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().
 382   //
 383   // Author : Alex Bercuci (A.Bercuci@gsi.de)
 384   //
 385
 386         if(!h){
 387                 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D2()", "No histogram defined.");
 388                 return 0.;
 389         }
 390
 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();
 396         Double_t p[2];
 397         Double_t entries;
 398
 399         // late construction of fitter
 400         if(!fFitter2D2) fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y");
 401
 402         fFitter2D2->ClearPoints();
 403         Int_t npoints=0;
 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));
 417                                 npoints++;
 418                         }
 419                 }
 420                 if(npoints>=25) break;
 421         }
 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);
 426
 427                 return 0.;
 428         }
 429         TVectorD vec(6);
 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];
 432         return exp(result);
 433
 434 }
 435
 436 //__________________________________________________________________
 437 Double_t AliTRDCalPIDRefMaker::Estimate2D1(TH2 *h, Float_t &x, Float_t &y
 438                                       , Float_t &dCT, Float_t &rmin
 439                                       , Float_t &rmax)
 440 {
 441   //
 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 !
 446   //
 447
 448         if(!h){
 449                 AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D1()", "No histogram defined.");
 450                 return 0.;
 451         }
 452
 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();
 458         Double_t p[2];
 459         Double_t entries;
 460         Double_t rxy = sqrt(x*x + y*y), rpxy;
 461
 462         // late construction of fitter
 463         if(!fFitter2D1) fFitter2D1 = new TLinearFitter(3, "1++x++y");
 464
 465         fFitter2D1->ClearPoints();
 466         Int_t npoints=0;
 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;
 475
 476                         fFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries));
 477                         npoints++;
 478                 }
 479         }
 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);
 484                 return 0.;
 485         }
 486         TVectorD vec(3);
 487         fFitter2D1->GetParameters(vec);
 488         Double_t result = vec[0] + x*vec[1] + y*vec[2];
 489         return exp(result);
 490 }
 491
 492 //__________________________________________________________________
 493 // Double_t     AliTRDCalPIDRefMaker::Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z)
 494 // {
 495 //      // Author Alex Bercuci (A.Bercuci@gsi.de)
 496 //      return 0.;
 497 // }
 498
 499
 500
 501 /////////////  Private functions ///////////////////////////////////
 502
 503 //__________________________________________________________________
 504 Int_t AliTRDCalPIDRefMaker::CheckProdDirTree(Char_t *dir)
 505 {
 506         // Scan directory tree for momenta. Returns the smallest number of
 507         // batches found in all directories or 0 if one momentum is missing.
 508
 509         const char *pwd = gSystem->pwd();
 510
 511         if(!gSystem->ChangeDirectory(dir)){
 512                 AliError(Form("Couldn't access production root directory %s.", dir));
 513                 return 0;
 514         }
 515
 516         Int_t iDir;
 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)));
 521                         return 0;
 522                 }
 523
 524                 iDir = 0;
 525                 while(gSystem->ChangeDirectory(Form("%03d", iDir))){
 526                         iDir++;
 527                         gSystem->ChangeDirectory("..");
 528                 }
 529                 if(iDir < nDir) nDir = iDir;
 530                 gSystem->ChangeDirectory(dir);
 531         }
 532
 533         gSystem->ChangeDirectory(pwd);
 534
 535         return nDir;
 536 }
 537
 538
 539 //__________________________________________________________________
 540 void  AliTRDCalPIDRefMaker::Prepare2D()
 541 {
 542   //
 543   // Prepares the 2-dimensional reference histograms
 544   //
 545
 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++){
 551                 // check PCA data
 552                 if(!fPrinc[ispec]){
 553                         AliError(Form("No data defined for %s.", AliTRDCalPIDLQ::fpartName[ispec]));
 554                         return;
 555                 }
 556                 // build reference histograms
 557                 nBinsSector = 10;
 558                 nSectors = 20;
 559                 nbinsx = nbinsy = nBinsSector * nSectors;
 560                 xmin = ymin = 0.;
 561                 xmax = 8000.; ymax = 6000.;
 562                 if(!h2dEdx[ispec]){
 563                         h2dEdx[ispec] = new  TH2D(Form("h2%s", AliTRDCalPIDLQ::fpartSymb[ispec]), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
 564                         h2dEdx[ispec]->SetLineColor(color[ispec]);
 565                 }
 566         }
 567
 568         // build transformed rotated histograms
 569         nbinsx = nbinsy = 100;
 570         xmin = ymin = -6.;
 571         xmax = 10.; ymax = 6.;
 572         TH2I *hProj = 0x0;
 573         if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new  TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
 574         else hProj->Reset();
 575
 576         printf("Doing interpolation and invertion ... "); fflush(stdout);
 577         Bool_t kDebugPlot = kTRUE;
 578         //Bool_t kDebugPrint = kFALSE;
 579
 580         TCanvas *c = 0x0;
 581         TEllipse *ellipse = 0x0;
 582         TMarker *mark = 0x0;
 583         if(kDebugPlot){
 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);
 589         }
 590
 591         // define observable variables
 592         TAxis *ax, *ay;
 593         Double_t xy[2], lxy[2], rxy[2];
 594         Double_t estimate, position;
 595         const TVectorD *eValues;
 596
 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;
 610         }
 611         std::vector<UShort_t> storeX[nPhi], storeY[nPhi];
 612
 613         // define working variables
 614         Float_t x0, y0, rx, ry, rc, rmin, rmax, dr, dCT;
 615         Double_t Phi, r;
 616         Int_t iPhi;
 617         Double_t entries;
 618         for(int ispec=0; ispec<5; ispec++){
 619                 hProj->Reset();
 620                 for(int iphi=0; iphi<nPhi; iphi++){
 621                         nPoints[iphi] = 0;
 622                         refsFitter[iphi].ClearPoints();
 623                         fFitterRange[iphi] = .8;
 624                         kFitterStatus[iphi] = kFALSE;
 625                         storeX[iphi].clear();
 626                         storeY[iphi].clear();
 627                 }
 628
 629                 // calculate covariance ellipse
 630                 fPrinc[ispec]->MakePrincipals();
 631                 eValues  = fPrinc[ispec]->GetEigenValues();
 632                 x0  = 0.;
 633                 y0  = 0.;
 634                 rx  = 3.5*sqrt((*eValues)[0]);
 635                 ry  = 3.5*sqrt((*eValues)[1]);
 636
 637                 // rotate to principal axis
 638                 Int_t irow = 0;
 639                 const Double_t *xx;
 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]);
 644                 }
 645
 646
 647 //      // debug plot
 648 //              if(kDebugPlot){
 649 //                      hProj->Draw();
 650 //                      ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.);
 651 //                      mark->DrawMarker(x0, y0);
 652 //                      gPad->Modified(); gPad->Update();
 653 //              }
 654
 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);
 662
 663                                 if((entries = hProj->GetBinContent(ibin, jbin)) == 0) continue;
 664
 665                                 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
 666                                 if(position < 1.) continue;
 667
 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);
 671
 672                                 refsFitter[iPhi].AddPoint(rxy, log(entries), 1./sqrt(entries));
 673                                 nPoints[iPhi]++;
 674                         }
 675                 }
 676
 677                 // do interpolation
 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);
 682                         lxy[0] = log(xy[0]);
 683                         for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
 684                                 xy[1]  = ay->GetBinCenter(jbin);
 685                                 lxy[1] = log(xy[1]);
 686
 687                                 // rotate to PCA
 688                                 fPrinc[ispec]->X2P(lxy, rxy);
 689
 690                                 // calculate border of covariance ellipse
 691                                 position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry;
 692
 693                                 // interpolation inside the covariance ellipse
 694                                 if(position < 1.){
 695                                   Float_t xTemp = rxy[0];
 696                                   Float_t yTemp = rxy[1];
 697                                         estimate = Estimate2D2((TH2 *) hProj, xTemp, yTemp);
 698                                         rxy[0] = xTemp;
 699                                         rxy[1] = yTemp;
 700                                         //                                      estimate = Estimate2D2((TH2 *) hProj, (Float_t)rxy[0], (Float_t)rxy[1]);
 701                                         h2dEdx[ispec]->SetBinContent(ibin, jbin, estimate/xy[0]/xy[1]);
 702                                 } else { // interpolation outside the covariance ellipse
 703                                         r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]);
 704                                         Phi   = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi]
 705                                         iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1);
 706
 707                                         storeX[iPhi].push_back(ibin);
 708                                         storeY[iPhi].push_back(jbin);
 709                                 }
 710                         }
 711                 }
 712
 713                 // Fill outliers
 714                 // Radial fit on transformed rotated
 715                 TVectorD vec(3);
 716                 Int_t xbin, ybin;
 717                 for(int iphi=0; iphi<nPhi; iphi++){
 718                         Phi = iphi * dPhi - TMath::Pi();
 719                         if(fabs(fabs(Phi)-TMath::Pi()) < 100.*TMath::DegToRad()) continue;
 720
 721
 722                         refsFitter[iphi].Eval();
 723                         refsFitter[iphi].GetParameters(vec);
 724                         for(int ipoint=0; ipoint<storeX[iphi].size(); ipoint++){
 725                                 xbin = storeX[iphi].at(ipoint);
 726                                 ybin = storeY[iphi].at(ipoint);
 727                                 xy[0] = ax->GetBinCenter(xbin); lxy[0] = log(xy[0]);
 728                                 xy[1] = ay->GetBinCenter(ybin); lxy[1] = log(xy[1]);
 729
 730                                 // rotate to PCA
 731                                 fPrinc[ispec]->X2P(lxy, rxy);
 732                                 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]);
 733                                 h2dEdx[ispec]->SetBinContent(xbin, ybin, estimate/xy[0]/xy[1]);
 734                         }
 735                 }
 736
 737                 // Longitudinal fit on ref histo in y direction
 738                 for(int isector=1; isector<nSectors; isector++){
 739                         // define sectors
 740                         binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
 741                         binStop  = (isector+1)*nBinsSector;
 742
 743                         nPoints[0] = 0;
 744                         refsLongFitter.ClearPoints();
 745                         storeX[0].clear();
 746                         storeY[0].clear();
 747
 748                         for(int ibin=1; ibin<=ax->GetNbins(); ibin++){
 749                                 xy[0] = ax->GetBinCenter(ibin);
 750                                 for(int jbin=binStart; jbin<=binStop; jbin++){
 751                                         xy[1] = ay->GetBinCenter(jbin);
 752                                         if((entries = h2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
 753                                                 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
 754                                                 nPoints[0]++;
 755                                         } else {
 756                                                 storeX[0].push_back(ibin);
 757                                                 storeY[0].push_back(jbin);
 758                                         }
 759                                 }
 760                                 if(nPoints[0] >= ((isector==0)?60:420)) break;
 761                         }
 762
 763
 764                         if(refsLongFitter.Eval() == 1){
 765                                 printf("Error in Y sector %d\n", isector);
 766                                 continue;
 767                         }
 768                         refsLongFitter.GetParameters(vec);
 769                         for(int ipoint=0; ipoint<storeX[0].size(); ipoint++){
 770                                 xbin = storeX[0].at(ipoint);
 771                                 ybin = storeY[0].at(ipoint);
 772                                 xy[0] = ax->GetBinCenter(xbin);
 773                                 xy[1] = ay->GetBinCenter(ybin);
 774
 775                                 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];
 776                                 h2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
 777                         }
 778                 }
 779
 780                 // Longitudinal fit on ref histo in x direction
 781                 for(int isector=1; isector<nSectors; isector++){
 782                         binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1;
 783                         binStop  = (isector+1)*nBinsSector;
 784
 785                         nPoints[0] = 0;
 786                         nFitPoints = 0;
 787                         refsLongFitter.ClearPoints();
 788                         storeX[0].clear();
 789                         storeY[0].clear();
 790
 791                         for(int jbin=1; jbin<=ay->GetNbins(); jbin++){
 792                                 xy[1] = ay->GetBinCenter(jbin);
 793                                 for(int ibin=binStart; ibin<=binStop; ibin++){
 794                                         xy[0] = ax->GetBinCenter(ibin);
 795                                         if((entries = h2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){
 796                                                 refsLongFitter.AddPoint(xy, log(entries), 1.e-7);
 797                                                 nPoints[0]++;
 798                                         } else {
 799                                                 storeX[0].push_back(ibin);
 800                                                 storeY[0].push_back(jbin);
 801                                                 nFitPoints++;
 802                                         }
 803                                 }
 804                                 if(nPoints[0] >= ((isector==0)?60:420)) break;
 805                         }
 806                         if(nFitPoints == 0) continue;
 807
 808
 809                         if(refsLongFitter.Eval() == 1){
 810                                 printf("Error in X sector %d\n", isector);
 811                                 continue;
 812                         }
 813                         refsLongFitter.GetParameters(vec);
 814                         for(int ipoint=0; ipoint<storeX[0].size(); ipoint++){
 815                                 xbin = storeX[0].at(ipoint);
 816                                 ybin = storeY[0].at(ipoint);
 817                                 xy[0] = ax->GetBinCenter(xbin);
 818                                 xy[1] = ay->GetBinCenter(ybin);
 819
 820                                 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];
 821                                 h2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate));
 822                         }
 823                 }
 824
 825
 826                 h2dEdx[ispec]->Scale(1./h2dEdx[ispec]->Integral());
 827
 828                 // debug plot
 829                 if(kDebugPlot){
 830                         h2dEdx[ispec]->Draw("cont1"); gPad->SetLogz();
 831                         gPad->Modified(); gPad->Update();
 832                 }
 833         }
 834         AliInfo("Finish interpolation.");
 835 }
 836
 837 //__________________________________________________________________
 838 void    AliTRDCalPIDRefMaker::Reset()
 839 {
 840   //
 841   // Reset reference histograms
 842   //
 843
 844         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 845                 if(h2dEdx[ispec]) h2dEdx[ispec]->Reset();
 846                 fPrinc[ispec]->Clear();
 847         }
 848         if(h1TB[0]) h1TB[0]->Reset();
 849         if(h1TB[1]) h1TB[1]->Reset();
 850 }
 851
 852 //__________________________________________________________________
 853 void  AliTRDCalPIDRefMaker::SaveReferences(const Int_t mom, const char *fn)
 854 {
 855   //
 856   // Save the reference histograms
 857   //
 858
 859         TFile *fSave = 0x0;
 860         TListIter it((TList*)gROOT->GetListOfFiles());
 861         Bool_t kFOUND = kFALSE;
 862         TDirectory *pwd = gDirectory;
 863         while((fSave=(TFile*)it.Next()))
 864                 if(strcmp(fn, fSave->GetName())==0){
 865                         kFOUND = kTRUE;
 866                         break;
 867                 }
 868         if(!kFOUND) fSave = TFile::Open(fn, "RECREATE");
 869         fSave->cd();
 870
 871         Float_t fmom = AliTRDCalPIDLQ::GetMomentum(mom);
 872
 873         // save dE/dx references
 874         TH2 *h2 = 0x0;
 875         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 876                 h2 = (TH2D*)h2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPIDLQ::fpartSymb[ispec], mom));
 877                 h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPIDLQ::fpartName[ispec], mom));
 878                 h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]");
 879                 h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]");
 880                 h2->GetZaxis()->SetTitle("Entries");
 881                 h2->Write();
 882         }
 883
 884         // save maximum time bin references
 885         TH1 *h1 = 0x0;
 886         h1 = (TH1F*)h1TB[0]->Clone(Form("h1MaxTBEL%02d", mom));
 887         h1->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fmom));
 888         h1->GetXaxis()->SetTitle("time [100 ns]");
 889         h1->GetYaxis()->SetTitle("Probability");
 890         h1->Write();
 891
 892         h1 = (TH1F*)h1TB[1]->Clone(Form("h1MaxTBPI%02d", mom));
 893         h1->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fmom));
 894         h1->GetXaxis()->SetTitle("time [100 ns]");
 895         h1->GetYaxis()->SetTitle("Probability");
 896         h1->Write();
 897
 898
 899         pwd->cd();
 900 }
 901