TOF/AliTOFcalib.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 /*$Log$
  17 author: Chiara Zampolli, zampolli@bo.infn.it
  18  */
  19
  20 ///////////////////////////////////////////////////////////////////////////////
  21 //                                                                           //
  22 // class for TOF calibration                                                 //
  23 //                                                                           //
  24 ///////////////////////////////////////////////////////////////////////////////
  25
  26 #include "AliTOFcalib.h"
  27 #include "AliRun.h"
  28 #include <TTask.h>
  29 #include <TFile.h>
  30 #include <TROOT.h>
  31 #include <TSystem.h>
  32 #include "AliTOF.h"
  33 #include "AliTOFcalibESD.h"
  34 #include "AliESD.h"
  35 #include <TObject.h>
  36 #include "TF1.h"
  37 #include "TH1F.h"
  38 #include "TH2F.h"
  39 #include "AliESDtrack.h"
  40 #include "AliTOFChannel.h"
  41 #include "AliTOFChSim.h"
  42 #include "AliTOFGeometry.h"
  43 #include "AliTOFdigit.h"
  44 #include "TClonesArray.h"
  45 #include "AliTOFCal.h"
  46 #include "TRandom.h"
  47 #include "AliTOFcluster.h"
  48 #include "TList.h"
  49 #include "AliCDBManager.h"
  50 #include "AliCDBMetaData.h"
  51 #include "AliCDBStorage.h"
  52 #include "AliCDBId.h"
  53 #include "AliCDBEntry.h"
  54
  55 extern TROOT *gROOT;
  56 extern TStyle *gStyle;
  57
  58 ClassImp(AliTOFcalib)
  59
  60 const Int_t AliTOFcalib::fgkchannel = 5000;
  61 const Char_t* AliTOFcalib::ffile[6]={"$ALICE_ROOT/TOF/Spectra/spectrum0_1.root","$ALICE_ROOT/TOF/Spectra/spectrum0_2.root","$ALICE_ROOT/TOF/Spectra/spectrum1_1.root","$ALICE_ROOT/TOF/Spectra/spectrum1_2.root","$ALICE_ROOT/TOF/Spectra/spectrum2_1.root","$ALICE_ROOT/TOF/Spectra/spectrum2_2.root"};
  62 //_______________________________________________________________________
  63 AliTOFcalib::AliTOFcalib():TTask("AliTOFcalib",""){
  64   fNSector = 0;
  65   fNPlate  = 0;
  66   fNStripA = 0;
  67   fNStripB = 0;
  68   fNStripC = 0;
  69   fNpadZ = 0;
  70   fNpadX = 0;
  71   fsize = 0;
  72   fArrayToT = 0x0;
  73   fArrayTime = 0x0;
  74   flistFunc = 0x0;
  75   fTOFCal = 0x0;
  76   fESDsel = 0x0;
  77   for (Int_t i = 0;i<6;i++){
  78     fhToT[i]=0x0;
  79   }
  80   fGeom=0x0;
  81 }
  82 //____________________________________________________________________________
  83
  84 AliTOFcalib::AliTOFcalib(char* headerFile, Int_t nEvents):TTask("AliTOFcalib","")
  85 {
  86   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
  87   rl->CdGAFile();
  88   TFile *in=(TFile*)gFile;
  89   in->cd();
  90   fGeom = (AliTOFGeometry*)in->Get("TOFgeometry");
  91   fNSector = fGeom->NSectors();
  92   fNPlate  = fGeom->NPlates();
  93   fNStripA = fGeom->NStripA();
  94   fNStripB = fGeom->NStripB();
  95   fNStripC = fGeom->NStripC();
  96   fNpadZ = fGeom->NpadZ();
  97   fNpadX = fGeom->NpadX();
  98   fsize = 2*(fNStripC+fNStripB) + fNStripA;
  99   for (Int_t i = 0;i<6;i++){
 100     fhToT[i]=0x0;
 101   }
 102   fArrayToT = 0x0;
 103   fArrayTime = 0x0;
 104   flistFunc = 0x0;
 105   fNevents=nEvents;
 106   fHeadersFile=headerFile;
 107   fTOFCal = 0x0;
 108   fESDsel = 0x0;
 109   TFile* file = (TFile*) gROOT->GetFile(fHeadersFile.Data()) ;
 110   if(file == 0){
 111     if(fHeadersFile.Contains("rfio"))
 112       file = TFile::Open(fHeadersFile,"update") ;
 113     else
 114       file = new TFile(fHeadersFile.Data(),"update") ;
 115     gAlice = (AliRun*)file->Get("gAlice") ;
 116   }
 117
 118   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
 119   roottasks->Add(this) ;
 120 }
 121 //____________________________________________________________________________
 122
 123 AliTOFcalib::AliTOFcalib(const AliTOFcalib & calib):TTask("AliTOFcalib","")
 124 {
 125   fNSector = calib.fNSector;
 126   fNPlate = calib.fNPlate;
 127   fNStripA = calib.fNStripA;
 128   fNStripB = calib.fNStripB;
 129   fNStripC = calib.fNStripC;
 130   fNpadZ = calib.fNpadZ;
 131   fNpadX = calib.fNpadX;
 132   fsize = calib.fsize;
 133   fArrayToT = calib.fArrayToT;
 134   fArrayTime = calib.fArrayTime;
 135   flistFunc = calib.flistFunc;
 136   for (Int_t i = 0;i<6;i++){
 137     fhToT[i]=calib.fhToT[i];
 138   }
 139   fTOFCal=calib.fTOFCal;
 140   fESDsel = calib.fESDsel;
 141   fGeom = calib.fGeom;
 142 }
 143
 144 //____________________________________________________________________________
 145
 146 AliTOFcalib::~AliTOFcalib()
 147 {
 148   delete fArrayToT;
 149   delete fArrayTime;
 150   delete flistFunc;
 151   delete[] fhToT;
 152   delete fTOFCal;
 153   delete fESDsel;
 154 }
 155 //____________________________________________________________________________
 156
 157 void AliTOFcalib::Init(){
 158   SetHistos();
 159   SetFitFunctions();
 160   fTOFCal = new AliTOFCal();
 161   fTOFCal->CreateArray();
 162   fNSector = 18;
 163   fNPlate  = 5;
 164   fNStripA = 15;
 165   fNStripB = 19;
 166   fNStripC = 20;
 167   fNpadZ = 2;
 168   fNpadX = 96;
 169   fsize = 1;
 170   //fsize = fNSector*(2*(fNStripC+fNStripB)+fNStripA())*fNpadZ*fNpadX; //generalized version
 171 }
 172 //____________________________________________________________________________
 173
 174 void AliTOFcalib::SetHistos(){
 175   TFile * spFile;
 176   TH1F * hToT;
 177   for (Int_t i =0;i<6;i++){
 178     //for the time being, only one spectrum is used
 179     spFile = new TFile("$ALICE_ROOT/TOF/Spectra/spectrumtest0_1.root","read");
 180     //otherwise
 181     //spFile = new TFile(ffile[i],"read");
 182     spFile->GetObject("ToT",hToT);
 183     fhToT[i]=hToT;
 184     Int_t nbins = hToT->GetNbinsX();
 185     Float_t Delta = hToT->GetBinWidth(1);
 186     Float_t maxch = hToT->GetBinLowEdge(nbins)+Delta;
 187     Float_t minch = hToT->GetBinLowEdge(1);
 188     Float_t max=0,min=0; //maximum and minimum value of the distribution
 189     Int_t maxbin=0,minbin=0; //maximum and minimum bin of the distribution
 190     for (Int_t ii=nbins; ii>0; ii--){
 191       if (hToT->GetBinContent(ii)!= 0) {
 192         max = maxch - (nbins-ii-1)*Delta;
 193         maxbin = ii;
 194         break;}
 195     }
 196     for (Int_t j=1; j<nbins; j++){
 197       if (hToT->GetBinContent(j)!= 0) {
 198         min = minch + (j-1)*Delta;
 199         minbin = j;
 200         break;}
 201     }
 202     fMaxToT[i]=max;
 203     fMinToT[i]=min;
 204     fMaxToTDistr[i]=hToT->GetMaximum();
 205   }
 206 }
 207 //__________________________________________________________________________
 208
 209 void AliTOFcalib::SetFitFunctions(){
 210   TFile * spFile;
 211   flistFunc = new TList();
 212   for (Int_t i =0;i<6;i++){
 213     //only one type of spectrum used for the time being
 214     spFile = new TFile("$ALICE_ROOT/TOF/Spectra/spectrumtest0_1.root","read");
 215     //otherwise
 216     //spFile = new TFile(ffile[i],"read");
 217     TH1F * h = (TH1F*)spFile->Get("TimeToTFit");
 218     TList * list = (TList*)h->GetListOfFunctions();
 219     TF1* f = (TF1*)list->At(0);
 220     Double_t par[6] = {0,0,0,0,0,0};
 221     Int_t npar=f->GetNpar();
 222     for (Int_t ipar=0;ipar<npar;ipar++){
 223       par[ipar]=f->GetParameter(ipar);
 224     }
 225     flistFunc->AddLast(f);
 226   }
 227   return;
 228 }
 229 //__________________________________________________________________________
 230
 231 TF1* AliTOFcalib::SetFitFunctions(TH1F *histo){
 232   TF1 * fpol[3];
 233   const Int_t nbins = histo->GetNbinsX();
 234   Float_t Delta = histo->GetBinWidth(1);  //all the bins have the same width
 235   Double_t max = histo->GetBinLowEdge(nbins)+Delta;
 236   max = 15;
 237   fpol[0]=new TF1("poly3","pol3",5,max);
 238   fpol[1]=new TF1("poly4","pol4",5,max);
 239   fpol[2]=new TF1("poly5","pol5",5,max);
 240   char npoly[10];
 241   Double_t chi[3]={1E6,1E6,1E6};
 242   Int_t ndf[3]={-1,-1,-1};
 243   Double_t Nchi[3]={1E6,1E6,1E6};
 244   Double_t bestchi=1E6;
 245   TF1 * fGold=0x0;
 246   Int_t nonzero =0;
 247   Int_t numberOfpar =0;
 248   for (Int_t j=0; j<nbins; j++){
 249     if (histo->GetBinContent(j)!=0) {
 250       nonzero++;
 251     }
 252   }
 253   Int_t norderfit = 0;
 254   if (nonzero<=4) {
 255     AliError(" Too few points in the histo. No fit performed.");
 256     return 0x0;
 257   }
 258   else if (nonzero<=5) {
 259     norderfit = 3;
 260     AliInfo(" Only 3rd order polynomial fit possible.");
 261   }
 262   else if (nonzero<=6) {
 263     norderfit = 4;
 264     AliInfo(" Only 3rd and 4th order polynomial fit possible.");
 265   }
 266   else {
 267     norderfit = 5;
 268     AliInfo(" All 3rd, 4th and 5th order polynomial fit possible.");
 269   }
 270   for (Int_t ifun=norderfit-3;ifun<norderfit-2;ifun++){
 271     sprintf(npoly,"poly%i",ifun+3);
 272     histo->Fit(npoly, "ERN", " ", 5.,14.);
 273     chi[ifun] = fpol[ifun]->GetChisquare();
 274     ndf[ifun] = fpol[ifun]->GetNDF();
 275     Nchi[ifun] = (Double_t)chi[ifun]/ndf[ifun];
 276     if (Nchi[ifun]<bestchi) {
 277       bestchi=Nchi[ifun];
 278       fGold = fpol[ifun];
 279       numberOfpar = fGold->GetNpar();
 280     }
 281   }
 282   fGold=fpol[2];  //Gold fit function set to pol5 in any case
 283   histo->Fit(fGold,"ER " ," ",5.,15.);
 284   return fGold;
 285 }
 286 //____________________________________________________________________________
 287
 288 TClonesArray* AliTOFcalib::DecalibrateDigits(TClonesArray * digits){
 289
 290   TObjArray ChArray(fsize);
 291   ChArray.SetOwner();
 292   for (Int_t kk = 0 ; kk < fsize; kk++){
 293     AliTOFChSim * channel = new AliTOFChSim();
 294     ChArray.Add(channel);
 295   }
 296   Int_t ndigits = digits->GetEntriesFast();
 297   for (Int_t i=0;i<ndigits;i++){
 298     Float_t trix = 0;
 299     Float_t triy = 0;
 300     Float_t SimToT = 0;
 301     AliTOFdigit * element = (AliTOFdigit*)digits->At(i);
 302     /*
 303     Int_t *detId[5];
 304     detId[0] = element->GetSector();
 305     detId[1] = element->GetPlate();
 306     detId[2] = element->GetStrip();
 307     detId[3] = element->GetPadz();
 308     detId[4] = element->GetPadx();
 309     Int_t index = GetIndex(detId);
 310     */
 311     //select the corresponding channel with its simulated ToT spectrum
 312     //summing up everything, index = 0 for all channels:
 313     Int_t index = 0;
 314     AliTOFChSim *ch = (AliTOFChSim*)ChArray.At(index);
 315     Int_t itype = -1;
 316     if (!ch->IsSlewed()){
 317       Double_t rand = gRandom->Uniform(0,6);
 318       itype = (Int_t)(6-rand);
 319       ch->SetSpectrum(itype);
 320       ch->SetSlewedStatus(kTRUE);
 321     }
 322     else itype = ch->GetSpectrum();
 323     TH1F * hToT = fhToT[itype];
 324     TF1 * f = (TF1*)flistFunc->At(itype);
 325     while (SimToT <= triy){
 326       trix = gRandom->Rndm(i);
 327       triy = gRandom->Rndm(i);
 328       trix = (fMaxToT[itype]-fMinToT[itype])*trix + fMinToT[itype];
 329       triy = fMaxToTDistr[itype]*triy;
 330       Int_t binx=hToT->FindBin(trix);
 331       SimToT=hToT->GetBinContent(binx);
 332     }
 333
 334     Float_t par[6];
 335     for(Int_t kk=0;kk<6;kk++){
 336       par[kk]=0;
 337     }
 338     element->SetToT(trix);
 339     Int_t nfpar = f->GetNpar();
 340     for (Int_t kk = 0; kk< nfpar; kk++){
 341       par[kk]=f->GetParameter(kk);
 342     }
 343     Float_t ToT = element->GetToT();
 344     element->SetTdcND(element->GetTdc());
 345     Float_t tdc = ((element->GetTdc())*AliTOFGeometry::TdcBinWidth()+32)*1.E-3; //tof signal in ns
 346     Float_t timeoffset=par[0] + ToT*par[1] +ToT*ToT*par[2] +ToT*ToT*ToT*par[3] +ToT*ToT*ToT*ToT*par[4] +ToT*ToT*ToT*ToT*ToT*par[5];
 347     timeoffset=0;
 348     Float_t timeSlewed = tdc + timeoffset;
 349     element->SetTdc((timeSlewed*1E3-32)/AliTOFGeometry::TdcBinWidth());
 350   }
 351   TFile * file = new TFile("ToT.root", "recreate");
 352   file->cd();
 353   ChArray.Write("ToT",TObject::kSingleKey);
 354   file->Close();
 355   delete file;
 356   ChArray.Clear();
 357   return digits;
 358 }
 359 //____________________________________________________________________________
 360
 361 void AliTOFcalib::SelectESD(AliESD *event, AliRunLoader *rl)
 362 {
 363   AliLoader *tofl = rl->GetLoader("TOFLoader");
 364   if (tofl == 0x0) {
 365     AliError("Can not get the TOF Loader");
 366     delete rl;
 367   }
 368   tofl->LoadRecPoints("read");
 369   TClonesArray *fClusters =new TClonesArray("AliTOFcluster",  1000);
 370   TTree *rTree=tofl->TreeR();
 371   if (rTree == 0) {
 372     cerr<<"Can't get the Rec Points tree !\n";
 373     delete rl;
 374   }
 375   TBranch *branch=rTree->GetBranch("TOF");
 376   if (!branch) {
 377     cerr<<"Cant' get the branch !\n";
 378     delete rl;
 379   }
 380   branch->SetAddress(&fClusters);
 381   rTree->GetEvent(0);
 382   Float_t LowerMomBound=0.8; // [GeV/c] default value Pb-Pb
 383   Float_t UpperMomBound=1.8 ; // [GeV/c] default value Pb-Pb
 384   Int_t ntrk =0;
 385   Int_t ngoodtrkfinalToT = 0;
 386   ntrk=event->GetNumberOfTracks();
 387   fESDsel = new TObjArray(ntrk);
 388   fESDsel->SetOwner();
 389   TObjArray  UCdatatemp(ntrk);
 390   Int_t ngoodtrk = 0;
 391   Int_t ngoodtrkfinal = 0;
 392   Float_t mintime =1E6;
 393   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 394     AliESDtrack *t=event->GetTrack(itrk);
 395     //track selection: reconstrution to TOF:
 396     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
 397       continue;
 398     }
 399     //IsStartedTimeIntegral
 400     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
 401       continue;
 402     }
 403     Double_t time=t->GetTOFsignal();
 404     time*=1.E-3; // tof given in nanoseconds
 405     if(time <= mintime)mintime=time;
 406     Double_t mom=t->GetP();
 407     if (!(mom<=UpperMomBound && mom>=LowerMomBound))continue;
 408     UInt_t AssignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
 409     if(AssignedTOFcluster==0){ // not matched
 410       continue;
 411     }
 412     AliTOFcluster * tofcl = (AliTOFcluster *)fClusters->UncheckedAt(AssignedTOFcluster);
 413     Int_t isector = tofcl->GetDetInd(0);
 414     Int_t iplate = tofcl->GetDetInd(1);
 415     Int_t istrip = tofcl->GetDetInd(2);
 416     Int_t ipadz = tofcl->GetDetInd(3);
 417     Int_t ipadx = tofcl->GetDetInd(4);
 418     Float_t ToT = tofcl->GetToT();
 419     AliTOFcalibESD *unc = new AliTOFcalibESD;
 420     unc->CopyFromAliESD(t);
 421     unc->SetSector(isector);
 422     unc->SetPlate(iplate);
 423     unc->SetStrip(istrip);
 424     unc->SetPadz(ipadz);
 425     unc->SetPadx(ipadx);
 426     unc->SetToT(ToT);
 427     Double_t c1[15];
 428     unc->GetExternalCovariance(c1);
 429     UCdatatemp.Add(unc);
 430     ngoodtrk++;
 431   }
 432   for (Int_t i = 0; i < ngoodtrk ; i ++){
 433     AliTOFcalibESD *unc = (AliTOFcalibESD*)UCdatatemp.At(i);
 434     if((unc->GetTOFsignal()-mintime*1.E3)<5.E3){
 435       fESDsel->Add(unc);
 436       ngoodtrkfinal++;
 437       ngoodtrkfinalToT++;
 438     }
 439   }
 440   fESDsel->Sort();
 441 }
 442 //_____________________________________________________________________________
 443
 444 void AliTOFcalib::CombESDId(){
 445   Float_t t0offset=0;
 446   Float_t loffset=0;
 447   Int_t   ntracksinset=6;
 448   Float_t exptof[6][3];
 449   Float_t momentum[6]={0.,0.,0.,0.,0.,0.};
 450   Int_t   assparticle[6]={3,3,3,3,3,3};
 451   Float_t massarray[3]={0.13957,0.493677,0.9382723};
 452   Float_t timeofflight[6]={0.,0.,0.,0.,0.,0.};
 453   Float_t beta[6]={0.,0.,0.,0.,0.,0.};
 454   Float_t texp[6]={0.,0.,0.,0.,0.,0.};
 455   Float_t sqMomError[6]={0.,0.,0.,0.,0.,0.};
 456   Float_t sqTrackError[6]={0.,0.,0.,0.,0.,0.};
 457   Float_t tracktoflen[6]={0.,0.,0.,0.,0.,0.};
 458   Float_t TimeResolution   = 0.90e-10; // 90 ps by default
 459   Float_t timeresolutioninns=TimeResolution*(1.e+9); // convert in [ns]
 460   Float_t timezero[6]={0.,0.,0.,0.,0.,0.};
 461   Float_t weightedtimezero[6]={0.,0.,0.,0.,0.,0.};
 462   Float_t besttimezero[6]={0.,0.,0.,0.,0.,0.};
 463   Float_t bestchisquare[6]={0.,0.,0.,0.,0.,0.};
 464   Float_t bestweightedtimezero[6]={0.,0.,0.,0.,0.,0.};
 465   Float_t bestsqTrackError[6]={0.,0.,0.,0.,0.,0.};
 466
 467   Int_t nelements = fESDsel->GetEntries();
 468   Int_t nset= (Int_t)(nelements/ntracksinset);
 469   for (Int_t i=0; i< nset; i++) {
 470
 471     AliTOFcalibESD **unc=new AliTOFcalibESD*[ntracksinset];
 472     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 473       Int_t index = itrk+i*ntracksinset;
 474       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
 475       unc[itrk]=element;
 476     }
 477
 478     for (Int_t j=0; j<ntracksinset; j++) {
 479       AliTOFcalibESD *element=unc[j];
 480       Double_t mom=element->GetP();
 481       Double_t time=element->GetTOFsignal()*1.E-3; // in ns
 482       Double_t exptime[10];
 483       element->GetIntegratedTimes(exptime);
 484       Double_t toflen=element->GetIntegratedLength()/100.;  // in m
 485       timeofflight[j]=time+t0offset;
 486       tracktoflen[j]=toflen+loffset;
 487       exptof[j][0]=exptime[2]*1.E-3+0.005;
 488       exptof[j][1]=exptime[3]*1.E-3+0.005;
 489       exptof[j][2]=exptime[4]*1.E-3+0.005;
 490       momentum[j]=mom;
 491     }
 492     Float_t t0best=999.;
 493     Float_t Et0best=999.;
 494     Float_t chisquarebest=999.;
 495     for (Int_t i1=0; i1<3;i1++) {
 496       beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
 497       texp[0]=exptof[0][i1];
 498       for (Int_t i2=0; i2<3;i2++) {
 499         beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
 500         texp[1]=exptof[1][i2];
 501         for (Int_t i3=0; i3<3;i3++) {
 502           beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
 503           texp[2]=exptof[2][i3];
 504           for (Int_t i4=0; i4<3;i4++) {
 505             beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
 506             texp[3]=exptof[3][i4];
 507
 508             for (Int_t i5=0; i5<3;i5++) {
 509               beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
 510               texp[4]=exptof[4][i5];
 511               for (Int_t i6=0; i6<3;i6++) {
 512                 beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
 513                 texp[5]=exptof[5][i6];
 514
 515                 Float_t sumAllweights=0.;
 516                 Float_t meantzero=0.;
 517                 Float_t Emeantzero=0.;
 518
 519                 for (Int_t itz=0; itz<ntracksinset;itz++) {
 520                   sqMomError[itz]=
 521                     ((1.-beta[itz]*beta[itz])*0.025)*
 522                     ((1.-beta[itz]*beta[itz])*0.025)*
 523                     (tracktoflen[itz]/
 524                      (0.299792*beta[itz]))*
 525                     (tracktoflen[itz]/
 526                      (0.299792*beta[itz]));
 527                   sqTrackError[itz]=
 528                     (timeresolutioninns*
 529                      timeresolutioninns
 530                      +sqMomError[itz]);
 531
 532                   timezero[itz]=texp[itz]-timeofflight[itz];
 533                   weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 534                   sumAllweights+=1./sqTrackError[itz];
 535                   meantzero+=weightedtimezero[itz];
 536
 537                 } // end loop for (Int_t itz=0; itz<15;itz++)
 538
 539                 meantzero=meantzero/sumAllweights; // it is given in [ns]
 540                 Emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
 541
 542                 // calculate chisquare
 543
 544                 Float_t chisquare=0.;
 545                 for (Int_t icsq=0; icsq<ntracksinset;icsq++) {
 546                   chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 547                 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 548                 //              cout << " chisquare " << chisquare << endl;
 549
 550                 Int_t npion=0;
 551                 if(i1==0)npion++;
 552                 if(i2==0)npion++;
 553                 if(i3==0)npion++;
 554                 if(i4==0)npion++;
 555                 if(i5==0)npion++;
 556                 if(i6==0)npion++;
 557
 558                 if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) ntracksinset)>0.3)){
 559                   //  if(chisquare<=chisquarebest){
 560
 561                   for(Int_t iqsq = 0; iqsq<ntracksinset; iqsq++) {
 562                     bestsqTrackError[iqsq]=sqTrackError[iqsq];
 563                     besttimezero[iqsq]=timezero[iqsq];
 564                     bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 565                     bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 566                   }
 567
 568                   assparticle[0]=i1;
 569                   assparticle[1]=i2;
 570                   assparticle[2]=i3;
 571                   assparticle[3]=i4;
 572                   assparticle[4]=i5;
 573                   assparticle[5]=i6;
 574
 575                   chisquarebest=chisquare;
 576                   t0best=meantzero;
 577                   Et0best=Emeantzero;
 578                 } // close if(dummychisquare<=chisquare)
 579               } // end loop on i6
 580             } // end loop on i5
 581           } // end loop on i4
 582         } // end loop on i3
 583       } // end loop on i2
 584     } // end loop on i1
 585
 586
 587     Float_t confLevel=999;
 588     if(chisquarebest<999.){
 589       Double_t dblechisquare=(Double_t)chisquarebest;
 590       confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinset-1);
 591     }
 592     // assume they are all pions for fake sets
 593     if(confLevel<0.01 || confLevel==999. ){
 594       for (Int_t itrk=0; itrk<ntracksinset; itrk++)assparticle[itrk]=0;
 595     }
 596     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 597       Int_t index = itrk+i*ntracksinset;
 598       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
 599       element->SetCombID(assparticle[itrk]);
 600     }
 601   }
 602 }
 603
 604 //_____________________________________________________________________________
 605
 606 void AliTOFcalib::CalibrateESD(){
 607   Int_t nelements = fESDsel->GetEntries();
 608   Int_t *number=new Int_t[fsize];
 609   fArrayToT = new AliTOFArray(fsize);
 610   fArrayTime = new AliTOFArray(fsize);
 611   for (Int_t i=0; i<fsize; i++){
 612     number[i]=0;
 613     fArrayToT->AddArray(i, new TArrayF(fgkchannel));
 614     TArrayF * parrToT = fArrayToT->GetArray(i);
 615     TArrayF & refaToT = * parrToT;
 616     fArrayTime->AddArray(i, new TArrayF(fgkchannel));
 617     TArrayF * parrTime = fArrayToT->GetArray(i);
 618     TArrayF & refaTime = * parrTime;
 619     for (Int_t j = 0;j<AliTOFcalib::fgkchannel;j++){
 620       refaToT[j]=0.;      //ToT[i][j]=j;
 621       refaTime[j]=0.;      //Time[i][j]=j;
 622     }
 623   }
 624
 625   for (Int_t i=0; i< nelements; i++) {
 626     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
 627     Int_t ipid = element->GetCombID();
 628     Double_t etime = 0;   //expected time
 629     Double_t expTime[10];
 630     element->GetIntegratedTimes(expTime);
 631     if (ipid == 0) etime = expTime[2]*1E-3; //ns
 632     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
 633     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
 634     else AliError("No pid from combinatorial algo for this track");
 635     Double_t mtime = (Double_t)element->GetTOFsignal()*1E-3;  //measured time
 636     Double_t mToT = (Double_t) element->GetToT();  //measured ToT, ns
 637     /*
 638     Int_t *detId[5];
 639     detId[0] = element->GetSector();
 640     detId[1] = element->GetPlate();
 641     detId[2] = element->GetStrip();
 642     detId[3] = element->GetPadz();
 643     detId[4] = element->GetPadx();
 644     Int_t index = GetIndex(detId);
 645     */
 646     //select the correspondent channel with its simulated ToT spectrum
 647     //summing up everything, index = 0 for all channels:
 648     Int_t index = 0;
 649     Int_t index2 = number[index];
 650     TArrayF * parrToT = fArrayToT->GetArray(index);
 651     TArrayF & refaToT = * parrToT;
 652     refaToT[index2] = (Float_t)mToT;
 653     TArrayF * parrTime = fArrayTime->GetArray(index);
 654     TArrayF & refaTime = * parrTime;
 655     refaTime[index2] = (Float_t)(mtime-etime);
 656     number[index]++;
 657   }
 658
 659   for (Int_t i=0;i<1;i++){
 660     TH1F * hProf = Profile(i);
 661     TF1* fGold = SetFitFunctions(hProf);
 662     Int_t nfpar = fGold->GetNpar();
 663     Float_t par[6];
 664     for(Int_t kk=0;kk<6;kk++){
 665       par[kk]=0;
 666     }
 667     for (Int_t kk = 0; kk< nfpar; kk++){
 668       par[kk]=fGold->GetParameter(kk);
 669     }
 670     AliTOFChannel * CalChannel = fTOFCal->GetChannel(i);
 671     CalChannel->SetSlewPar(par);
 672   }
 673   delete[] number;
 674 }
 675
 676 //___________________________________________________________________________
 677
 678 TH1F* AliTOFcalib::Profile(Int_t ich){
 679   const Int_t nbinToT = 650;
 680   Int_t nbinTime = 400;
 681   Float_t minTime = -10.5; //ns
 682   Float_t maxTime = 10.5; //ns
 683   Float_t minToT = 7.5; //ns
 684   Float_t maxToT = 40.; //ns
 685   Float_t DeltaToT = (maxToT-minToT)/nbinToT;
 686   Double_t mTime[nbinToT+1],mToT[nbinToT+1],meanTime[nbinToT+1], meanTime2[nbinToT+1],ToT[nbinToT+1], ToT2[nbinToT+1],meanToT[nbinToT+1],meanToT2[nbinToT+1],Time[nbinToT+1],Time2[nbinToT+1],xlow[nbinToT+1],sigmaTime[nbinToT+1];
 687   Int_t n[nbinToT+1], nentrx[nbinToT+1];
 688   Double_t sigmaToT[nbinToT+1];
 689   for (Int_t i = 0; i < nbinToT+1 ; i++){
 690     mTime[i]=0;
 691     mToT[i]=0;
 692     n[i]=0;
 693     meanTime[i]=0;
 694     meanTime2[i]=0;
 695     ToT[i]=0;
 696     ToT2[i]=0;
 697     meanToT[i]=0;
 698     meanToT2[i]=0;
 699     Time[i]=0;
 700     Time2[i]=0;
 701     xlow[i]=0;
 702     sigmaTime[i]=0;
 703     sigmaToT[i]=0;
 704     n[i]=0;
 705     nentrx[i]=0;
 706   }
 707   TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", nbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
 708   TArrayF * parrToT = fArrayToT->GetArray(ich);
 709   TArrayF & refaToT = * parrToT;
 710   TArrayF * parrTime = fArrayTime->GetArray(ich);
 711   TArrayF & refaTime = * parrTime;
 712   for (Int_t j = 0; j < AliTOFcalib::fgkchannel; j++){
 713     if (refaToT[j] == 0) continue;
 714     Int_t nx = (Int_t)((refaToT[j]-minToT)/DeltaToT)+1;
 715     if ((refaToT[j] != 0) && (refaTime[j] != 0)){
 716       Time[nx]+=refaTime[j];
 717       Time2[nx]+=(refaTime[j])*(refaTime[j]);
 718       ToT[nx]+=refaToT[j];
 719       ToT2[nx]+=refaToT[j]*refaToT[j];
 720       nentrx[nx]++;
 721       hSlewing->Fill(refaToT[j],refaTime[j]);
 722     }
 723   }
 724   Int_t nbinsToT=hSlewing->GetNbinsX();
 725   if (nbinsToT != nbinToT) {
 726     AliError("Profile :: incompatible numbers of bins");
 727     return 0x0;
 728   }
 729
 730   Int_t usefulBins=0;
 731   TH1F *histo = new TH1F("histo", "1D Time vs ToT", nbinsToT, minToT, maxToT);
 732   for (Int_t i=1;i<=nbinsToT;i++){
 733     if (nentrx[i]!=0){
 734     n[usefulBins]+=nentrx[i];
 735     if (n[usefulBins]==0 && i == nbinsToT) {
 736       break;
 737     }
 738     meanTime[usefulBins]+=Time[i];
 739     meanTime2[usefulBins]+=Time2[i];
 740     meanToT[usefulBins]+=ToT[i];
 741     meanToT2[usefulBins]+=ToT2[i];
 742     if (n[usefulBins]<20 && i!=nbinsToT) continue;
 743     mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
 744     mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
 745     sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
 746                                    *(meanTime2[usefulBins]-meanTime[usefulBins]
 747                                      *meanTime[usefulBins]/n[usefulBins]));
 748     if ((1./n[usefulBins]/n[usefulBins]
 749          *(meanToT2[usefulBins]-meanToT[usefulBins]
 750            *meanToT[usefulBins]/n[usefulBins]))< 0) {
 751       AliError(" too small radical" );
 752       sigmaToT[usefulBins]=0;
 753     }
 754     else{
 755       sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
 756                                      *(meanToT2[usefulBins]-meanToT[usefulBins]
 757                                        *meanToT[usefulBins]/n[usefulBins]));
 758     }
 759     usefulBins++;
 760     }
 761   }
 762   for (Int_t i=0;i<usefulBins;i++){
 763     Int_t binN = (Int_t)((mToT[i]-minToT)/DeltaToT)+1;
 764     histo->Fill(mToT[i],mTime[i]);
 765     histo->SetBinError(binN,sigmaTime[i]);
 766   }
 767   return histo;
 768 }
 769 //_____________________________________________________________________________
 770
 771 void AliTOFcalib::CorrectESDTime(){
 772   Int_t nelements = fESDsel->GetEntries();
 773   for (Int_t i=0; i< nelements; i++) {
 774     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
 775     /*
 776     Int_t *detId[5];
 777     detId[0] = element->GetSector();
 778     detId[1] = element->GetPlate();
 779     detId[2] = element->GetStrip();
 780     detId[3] = element->GetPadz();
 781     detId[4] = element->GetPadx();
 782     Int_t index = GetIndex(detId);
 783     */
 784     //select the correspondent channel with its simulated ToT spectrum
 785     //summing up everything, index = 0 for all channels:
 786     Int_t index = 0;
 787     Int_t ipid = element->GetCombID();
 788     Double_t etime = 0;   //expected time
 789     Double_t expTime[10];
 790     element->GetIntegratedTimes(expTime);
 791     if (ipid == 0) etime = expTime[2]*1E-3; //ns
 792     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
 793     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
 794     Float_t par[6];
 795     AliTOFChannel * CalChannel = fTOFCal->GetChannel(index);
 796     for (Int_t j = 0; j<6; j++){
 797       par[j]=CalChannel->GetSlewPar(j);
 798     }
 799     //Float_t TimeCorr=par[0]+par[1]*ToT+par[2]*ToT*ToT+par[3]*ToT*ToT*ToT+par[4]*ToT*ToT*ToT*ToT+par[5]*ToT*ToT*ToT*ToT*ToT;
 800   }
 801 }
 802 //_____________________________________________________________________________
 803
 804 void AliTOFcalib::CorrectESDTime(AliESD *event, AliRunLoader *rl ){
 805   AliLoader *tofl = rl->GetLoader("TOFLoader");
 806   if (tofl == 0x0) {
 807     AliError("Can not get the TOF Loader");
 808     delete rl;
 809   }
 810   tofl->LoadRecPoints("read");
 811   TClonesArray *fClusters =new TClonesArray("AliTOFcluster",  1000);
 812   TTree *rTree=tofl->TreeR();
 813   if (rTree == 0) {
 814     cerr<<"Can't get the Rec Points tree !\n";
 815     delete rl;
 816   }
 817   TBranch *branch=rTree->GetBranch("TOF");
 818   if (!branch) {
 819     cerr<<"Cant' get the branch !\n";
 820     delete rl;
 821   }
 822   branch->SetAddress(&fClusters);
 823   rTree->GetEvent(0);
 824   Int_t ntrk =0;
 825   ntrk=event->GetNumberOfTracks();
 826   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 827     AliESDtrack *t=event->GetTrack(itrk);
 828     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
 829       continue;
 830     }
 831     //IsStartedTimeIntegral
 832     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
 833       continue;
 834     }
 835     UInt_t AssignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
 836     if(AssignedTOFcluster==0){ // not matched
 837       continue;
 838     }
 839     AliTOFcluster * tofcl = (AliTOFcluster *)fClusters->UncheckedAt(AssignedTOFcluster);
 840     /*
 841     Int_t *detId[5];
 842     detId[0] = tofcl->GetSector();
 843     detId[1] = tofcl->GetPlate();
 844     detId[2] = tofcl->GetStrip();
 845     detId[3] = tofcl->GetPadz();
 846     detId[4] = tofcl->GetPadx();
 847     Int_t index = GetIndex(detId);
 848     */
 849     //select the correspondent channel with its simulated ToT spectrum
 850     //summing up everything, index = 0 for all channels:
 851     Int_t index = 0;
 852     AliTOFChannel * CalChannel = fTOFCal->GetChannel(index);
 853     Float_t par[6];
 854     for (Int_t j = 0; j<6; j++){
 855       par[j]=CalChannel->GetSlewPar(j);
 856     }
 857     Float_t ToT = tofcl->GetToT();
 858     Float_t TimeCorr =0;
 859     TimeCorr=par[0]+par[1]*ToT+par[2]*ToT*ToT+par[3]*ToT*ToT*ToT+par[4]*ToT*ToT*ToT*ToT+par[5]*ToT*ToT*ToT*ToT*ToT;
 860   }
 861 }
 862 //_____________________________________________________________________________
 863
 864 void AliTOFcalib::WriteOnCDB(){
 865   AliCDBManager *man = AliCDBManager::Instance();
 866   AliCDBId id("TOF/Calib/Constants",1,100);
 867   AliCDBMetaData *md = new AliCDBMetaData();
 868   md->SetResponsible("Shimmize");
 869   man->Put(fTOFCal,id,md);
 870 }
 871 //_____________________________________________________________________________
 872
 873 void AliTOFcalib::ReadFromCDB(Char_t *sel, Int_t nrun){
 874   AliCDBManager *man = AliCDBManager::Instance();
 875   AliCDBEntry *entry = man->Get(sel,nrun);
 876   if (!entry){
 877     AliError("Retrivial failed");
 878     AliCDBStorage *origSto =man->GetDefaultStorage();
 879     man->SetDefaultStorage("local://$ALICE_ROOT");
 880     entry = man->Get("TOF/Calib/Constants",nrun);
 881     man->SetDefaultStorage(origSto);
 882   }
 883   AliTOFCal *cal =(AliTOFCal*)entry->GetObject();
 884   fTOFCal = cal;
 885 }
 886 //_____________________________________________________________________________
 887
 888 Int_t AliTOFcalib::GetIndex(Int_t *detId){
 889   Int_t isector = detId[0];
 890   if (isector >= fNSector)
 891     AliError(Form("Wrong sector number in TOF (%d) !",isector));
 892   Int_t iplate = detId[1];
 893   if (iplate >= fNPlate)
 894     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 895   Int_t istrip = detId[2];
 896   Int_t ipadz = detId[3];
 897   Int_t ipadx = detId[4];
 898   Int_t stripOffset = 0;
 899   switch (iplate) {
 900   case 0:
 901     stripOffset = 0;
 902     break;
 903   case 1:
 904     stripOffset = fNStripC;
 905     break;
 906   case 2:
 907     stripOffset = fNStripC+fNStripB;
 908     break;
 909   case 3:
 910     stripOffset = fNStripC+fNStripB+fNStripA;
 911     break;
 912   case 4:
 913     stripOffset = fNStripC+fNStripB+fNStripA+fNStripB;
 914     break;
 915   default:
 916     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 917     break;
 918   };
 919
 920   Int_t idet = ((2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX)*isector +
 921                (stripOffset*fNpadZ*fNpadX)+
 922                (fNpadZ*fNpadX)*istrip+
 923                (fNpadX)*ipadz+
 924                 ipadx;
 925   return idet;
 926 }
 927