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