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.3  2006/03/28 14:57:02  arcelli
  19 updates to handle new V5 geometry & some re-arrangements
  20
  21 Revision 1.2  2006/02/13 17:22:26  arcelli
  22 just Fixing Log info
  23
  24 Revision 1.1  2006/02/13 16:10:48  arcelli
  25 Add classes for TOF Calibration (C.Zampolli)
  26
  27 author: Chiara Zampolli, zampolli@bo.infn.it
  28 */
  29
  30 ///////////////////////////////////////////////////////////////////////////////
  31 //                                                                           //
  32 // class for TOF calibration                                                 //
  33 //                                                                           //
  34 ///////////////////////////////////////////////////////////////////////////////
  35
  36 #include "AliTOFcalib.h"
  37 #include "AliRun.h"
  38 #include <TTask.h>
  39 #include <TFile.h>
  40 #include <TROOT.h>
  41 #include <TSystem.h>
  42 #include "AliTOF.h"
  43 #include "AliTOFcalibESD.h"
  44 #include "AliESD.h"
  45 #include <TObject.h>
  46 #include "TF1.h"
  47 #include "TH1F.h"
  48 #include "TH2F.h"
  49 #include "AliESDtrack.h"
  50 #include "AliTOFChannel.h"
  51 #include "AliTOFChSim.h"
  52 #include "AliTOFGeometryV5.h"
  53 #include "TClonesArray.h"
  54 #include "AliTOFCal.h"
  55 #include "TRandom.h"
  56 #include "AliTOFcluster.h"
  57 #include "TList.h"
  58 #include "AliCDBManager.h"
  59 #include "AliCDBMetaData.h"
  60 #include "AliCDBStorage.h"
  61 #include "AliCDBId.h"
  62 #include "AliCDBEntry.h"
  63
  64 extern TROOT *gROOT;
  65 extern TStyle *gStyle;
  66
  67 ClassImp(AliTOFcalib)
  68
  69 const Int_t AliTOFcalib::fgkchannel = 5000;
  70 //_______________________________________________________________________
  71 AliTOFcalib::AliTOFcalib():TTask("AliTOFcalib",""){
  72
  73   fArrayToT = 0x0;
  74   fArrayTime = 0x0;
  75   fESDsel = 0x0;
  76   AliTOFGeometry *geom=new AliTOFGeometryV5();
  77   AliInfo("V5 TOF Geometry is taken as the default");
  78   fNSector = geom->NSectors();
  79   fNPlate  = geom->NPlates();
  80   fNStripA = geom->NStripA();
  81   fNStripB = geom->NStripB();
  82   fNStripC = geom->NStripC();
  83   fNpadZ = geom->NpadZ();
  84   fNpadX = geom->NpadX();
  85   fNChannels = fNSector*(2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX; //generalized version
  86   fTOFCal = new AliTOFCal(geom);
  87   fTOFSimCal = new AliTOFCal(geom);
  88   fTOFCal->CreateArray();
  89   fTOFSimCal->CreateArray();
  90   delete geom;
  91 }
  92 //_______________________________________________________________________
  93 AliTOFcalib::AliTOFcalib(AliTOFGeometry *geom):TTask("AliTOFcalib",""){
  94   fArrayToT = 0x0;
  95   fArrayTime = 0x0;
  96   fESDsel = 0x0;
  97   fNSector = geom->NSectors();
  98   fNPlate  = geom->NPlates();
  99   fNStripA = geom->NStripA();
 100   fNStripB = geom->NStripB();
 101   fNStripC = geom->NStripC();
 102   fNpadZ = geom->NpadZ();
 103   fNpadX = geom->NpadX();
 104   fNChannels = fNSector*(2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX; //generalized version
 105   fTOFCal = new AliTOFCal(geom);
 106   fTOFSimCal = new AliTOFCal(geom);
 107   fTOFCal->CreateArray();
 108   fTOFSimCal->CreateArray();
 109 }
 110 //____________________________________________________________________________
 111
 112 AliTOFcalib::AliTOFcalib(const AliTOFcalib & calib):TTask("AliTOFcalib","")
 113 {
 114   fNSector = calib.fNSector;
 115   fNPlate = calib.fNPlate;
 116   fNStripA = calib.fNStripA;
 117   fNStripB = calib.fNStripB;
 118   fNStripC = calib.fNStripC;
 119   fNpadZ = calib.fNpadZ;
 120   fNpadX = calib.fNpadX;
 121   fNChannels = calib.fNChannels;
 122   fArrayToT = calib.fArrayToT;
 123   fArrayTime = calib.fArrayTime;
 124   fTOFCal=calib.fTOFCal;
 125   fTOFSimCal = calib.fTOFSimCal;
 126 }
 127
 128 //____________________________________________________________________________
 129
 130 AliTOFcalib::~AliTOFcalib()
 131 {
 132   delete fArrayToT;
 133   delete fArrayTime;
 134   delete fTOFCal;
 135   delete fTOFSimCal;
 136   delete fESDsel;
 137 }
 138 //__________________________________________________________________________
 139
 140 TF1* AliTOFcalib::SetFitFunctions(TH1F *histo){
 141   TF1 * fpol[3];
 142   const Int_t nbins = histo->GetNbinsX();
 143   Float_t Delta = histo->GetBinWidth(1);  //all the bins have the same width
 144   Double_t max = histo->GetBinLowEdge(nbins)+Delta;
 145   max = 15;
 146   fpol[0]=new TF1("poly3","pol3",5,max);
 147   fpol[1]=new TF1("poly4","pol4",5,max);
 148   fpol[2]=new TF1("poly5","pol5",5,max);
 149   char npoly[10];
 150   Double_t chi[3]={1E6,1E6,1E6};
 151   Int_t ndf[3]={-1,-1,-1};
 152   Double_t Nchi[3]={1E6,1E6,1E6};
 153   Double_t bestchi=1E6;
 154   TF1 * fGold=0x0;
 155   Int_t nonzero =0;
 156   Int_t numberOfpar =0;
 157   for (Int_t j=0; j<nbins; j++){
 158     if (histo->GetBinContent(j)!=0) {
 159       nonzero++;
 160     }
 161   }
 162   Int_t norderfit = 0;
 163   if (nonzero<=4) {
 164     AliError(" Too few points in the histo. No fit performed.");
 165     return 0x0;
 166   }
 167   else if (nonzero<=5) {
 168     norderfit = 3;
 169     AliInfo(" Only 3rd order polynomial fit possible.");
 170   }
 171   else if (nonzero<=6) {
 172     norderfit = 4;
 173     AliInfo(" Only 3rd and 4th order polynomial fit possible.");
 174   }
 175   else {
 176     norderfit = 5;
 177     AliInfo(" All 3rd, 4th and 5th order polynomial fit possible.");
 178   }
 179   for (Int_t ifun=norderfit-3;ifun<norderfit-2;ifun++){
 180     sprintf(npoly,"poly%i",ifun+3);
 181     histo->Fit(npoly, "ERN", " ", 5.,14.);
 182     chi[ifun] = fpol[ifun]->GetChisquare();
 183     ndf[ifun] = fpol[ifun]->GetNDF();
 184     Nchi[ifun] = (Double_t)chi[ifun]/ndf[ifun];
 185     if (Nchi[ifun]<bestchi) {
 186       bestchi=Nchi[ifun];
 187       fGold = fpol[ifun];
 188       numberOfpar = fGold->GetNpar();
 189     }
 190   }
 191   fGold=fpol[2];  //Gold fit function set to pol5 in any case
 192   histo->Fit(fGold,"ER " ," ",5.,15.);
 193   return fGold;
 194 }
 195 //____________________________________________________________________________
 196
 197 void AliTOFcalib::SelectESD(AliESD *event)
 198 {
 199   Float_t LowerMomBound=0.8; // [GeV/c] default value Pb-Pb
 200   Float_t UpperMomBound=1.8 ; // [GeV/c] default value Pb-Pb
 201   Int_t ntrk =0;
 202   Int_t ngoodtrkfinalToT = 0;
 203   ntrk=event->GetNumberOfTracks();
 204   fESDsel = new TObjArray(ntrk);
 205   fESDsel->SetOwner();
 206   TObjArray  UCdatatemp(ntrk);
 207   Int_t ngoodtrk = 0;
 208   Int_t ngoodtrkfinal = 0;
 209   Float_t mintime =1E6;
 210   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 211     AliESDtrack *t=event->GetTrack(itrk);
 212     //track selection: reconstrution to TOF:
 213     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
 214       continue;
 215     }
 216     //IsStartedTimeIntegral
 217     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
 218       continue;
 219     }
 220     Double_t time=t->GetTOFsignal();
 221     time*=1.E-3; // tof given in nanoseconds
 222     if(time <= mintime)mintime=time;
 223     Double_t mom=t->GetP();
 224     if (!(mom<=UpperMomBound && mom>=LowerMomBound))continue;
 225     UInt_t AssignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
 226     if(AssignedTOFcluster==0){ // not matched
 227       continue;
 228     }
 229     AliTOFcalibESD *unc = new AliTOFcalibESD;
 230     unc->CopyFromAliESD(t);
 231     Double_t c1[15];
 232     unc->GetExternalCovariance(c1);
 233     UCdatatemp.Add(unc);
 234     ngoodtrk++;
 235   }
 236   for (Int_t i = 0; i < ngoodtrk ; i ++){
 237     AliTOFcalibESD *unc = (AliTOFcalibESD*)UCdatatemp.At(i);
 238     if((unc->GetTOFsignal()-mintime*1.E3)<5.E3){
 239       fESDsel->Add(unc);
 240       ngoodtrkfinal++;
 241       ngoodtrkfinalToT++;
 242     }
 243   }
 244   fESDsel->Sort();
 245 }
 246 //_____________________________________________________________________________
 247
 248 void AliTOFcalib::CombESDId(){
 249   Float_t t0offset=0;
 250   Float_t loffset=0;
 251   Int_t   ntracksinset=6;
 252   Float_t exptof[6][3];
 253   Float_t momentum[6]={0.,0.,0.,0.,0.,0.};
 254   Int_t   assparticle[6]={3,3,3,3,3,3};
 255   Float_t massarray[3]={0.13957,0.493677,0.9382723};
 256   Float_t timeofflight[6]={0.,0.,0.,0.,0.,0.};
 257   Float_t beta[6]={0.,0.,0.,0.,0.,0.};
 258   Float_t texp[6]={0.,0.,0.,0.,0.,0.};
 259   Float_t sqMomError[6]={0.,0.,0.,0.,0.,0.};
 260   Float_t sqTrackError[6]={0.,0.,0.,0.,0.,0.};
 261   Float_t tracktoflen[6]={0.,0.,0.,0.,0.,0.};
 262   Float_t TimeResolution   = 0.90e-10; // 90 ps by default
 263   Float_t timeresolutioninns=TimeResolution*(1.e+9); // convert in [ns]
 264   Float_t timezero[6]={0.,0.,0.,0.,0.,0.};
 265   Float_t weightedtimezero[6]={0.,0.,0.,0.,0.,0.};
 266   Float_t besttimezero[6]={0.,0.,0.,0.,0.,0.};
 267   Float_t bestchisquare[6]={0.,0.,0.,0.,0.,0.};
 268   Float_t bestweightedtimezero[6]={0.,0.,0.,0.,0.,0.};
 269   Float_t bestsqTrackError[6]={0.,0.,0.,0.,0.,0.};
 270
 271   Int_t nelements = fESDsel->GetEntries();
 272   Int_t nset= (Int_t)(nelements/ntracksinset);
 273   for (Int_t i=0; i< nset; i++) {
 274
 275     AliTOFcalibESD **unc=new AliTOFcalibESD*[ntracksinset];
 276     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 277       Int_t index = itrk+i*ntracksinset;
 278       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
 279       unc[itrk]=element;
 280     }
 281
 282     for (Int_t j=0; j<ntracksinset; j++) {
 283       AliTOFcalibESD *element=unc[j];
 284       Double_t mom=element->GetP();
 285       Double_t time=element->GetTOFsignal()*1.E-3; // in ns
 286       Double_t exptime[10];
 287       element->GetIntegratedTimes(exptime);
 288       Double_t toflen=element->GetIntegratedLength()/100.;  // in m
 289       timeofflight[j]=time+t0offset;
 290       tracktoflen[j]=toflen+loffset;
 291       exptof[j][0]=exptime[2]*1.E-3+0.005;
 292       exptof[j][1]=exptime[3]*1.E-3+0.005;
 293       exptof[j][2]=exptime[4]*1.E-3+0.005;
 294       momentum[j]=mom;
 295     }
 296     Float_t t0best=999.;
 297     Float_t Et0best=999.;
 298     Float_t chisquarebest=999.;
 299     for (Int_t i1=0; i1<3;i1++) {
 300       beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
 301       texp[0]=exptof[0][i1];
 302       for (Int_t i2=0; i2<3;i2++) {
 303         beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
 304         texp[1]=exptof[1][i2];
 305         for (Int_t i3=0; i3<3;i3++) {
 306           beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
 307           texp[2]=exptof[2][i3];
 308           for (Int_t i4=0; i4<3;i4++) {
 309             beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
 310             texp[3]=exptof[3][i4];
 311
 312             for (Int_t i5=0; i5<3;i5++) {
 313               beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
 314               texp[4]=exptof[4][i5];
 315               for (Int_t i6=0; i6<3;i6++) {
 316                 beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
 317                 texp[5]=exptof[5][i6];
 318
 319                 Float_t sumAllweights=0.;
 320                 Float_t meantzero=0.;
 321                 Float_t Emeantzero=0.;
 322
 323                 for (Int_t itz=0; itz<ntracksinset;itz++) {
 324                   sqMomError[itz]=
 325                     ((1.-beta[itz]*beta[itz])*0.025)*
 326                     ((1.-beta[itz]*beta[itz])*0.025)*
 327                     (tracktoflen[itz]/
 328                      (0.299792*beta[itz]))*
 329                     (tracktoflen[itz]/
 330                      (0.299792*beta[itz]));
 331                   sqTrackError[itz]=
 332                     (timeresolutioninns*
 333                      timeresolutioninns
 334                      +sqMomError[itz]);
 335
 336                   timezero[itz]=texp[itz]-timeofflight[itz];
 337                   weightedtimezero[itz]=timezero[itz]/sqTrackError[itz];
 338                   sumAllweights+=1./sqTrackError[itz];
 339                   meantzero+=weightedtimezero[itz];
 340
 341                 } // end loop for (Int_t itz=0; itz<15;itz++)
 342
 343                 meantzero=meantzero/sumAllweights; // it is given in [ns]
 344                 Emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
 345
 346                 // calculate chisquare
 347
 348                 Float_t chisquare=0.;
 349                 for (Int_t icsq=0; icsq<ntracksinset;icsq++) {
 350                   chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/sqTrackError[icsq];
 351                 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
 352                 //              cout << " chisquare " << chisquare << endl;
 353
 354                 Int_t npion=0;
 355                 if(i1==0)npion++;
 356                 if(i2==0)npion++;
 357                 if(i3==0)npion++;
 358                 if(i4==0)npion++;
 359                 if(i5==0)npion++;
 360                 if(i6==0)npion++;
 361
 362                 if(chisquare<=chisquarebest  && ((Float_t) npion/ ((Float_t) ntracksinset)>0.3)){
 363                   //  if(chisquare<=chisquarebest){
 364
 365                   for(Int_t iqsq = 0; iqsq<ntracksinset; iqsq++) {
 366                     bestsqTrackError[iqsq]=sqTrackError[iqsq];
 367                     besttimezero[iqsq]=timezero[iqsq];
 368                     bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
 369                     bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/sqTrackError[iqsq];
 370                   }
 371
 372                   assparticle[0]=i1;
 373                   assparticle[1]=i2;
 374                   assparticle[2]=i3;
 375                   assparticle[3]=i4;
 376                   assparticle[4]=i5;
 377                   assparticle[5]=i6;
 378
 379                   chisquarebest=chisquare;
 380                   t0best=meantzero;
 381                   Et0best=Emeantzero;
 382                 } // close if(dummychisquare<=chisquare)
 383               } // end loop on i6
 384             } // end loop on i5
 385           } // end loop on i4
 386         } // end loop on i3
 387       } // end loop on i2
 388     } // end loop on i1
 389
 390
 391     Float_t confLevel=999;
 392     if(chisquarebest<999.){
 393       Double_t dblechisquare=(Double_t)chisquarebest;
 394       confLevel=(Float_t)TMath::Prob(dblechisquare,ntracksinset-1);
 395     }
 396     // assume they are all pions for fake sets
 397     if(confLevel<0.01 || confLevel==999. ){
 398       for (Int_t itrk=0; itrk<ntracksinset; itrk++)assparticle[itrk]=0;
 399     }
 400     for (Int_t itrk=0; itrk<ntracksinset; itrk++) {
 401       Int_t index = itrk+i*ntracksinset;
 402       AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(index);
 403       element->SetCombID(assparticle[itrk]);
 404     }
 405   }
 406 }
 407
 408 //_____________________________________________________________________________
 409
 410 void AliTOFcalib::CalibrateESD(){
 411   Int_t nelements = fESDsel->GetEntries();
 412   Int_t *number=new Int_t[fNChannels];
 413   fArrayToT = new AliTOFArray(fNChannels);
 414   fArrayTime = new AliTOFArray(fNChannels);
 415   for (Int_t i=0; i<fNChannels; i++){
 416     number[i]=0;
 417     fArrayToT->AddArray(i, new TArrayF(fgkchannel));
 418     TArrayF * parrToT = fArrayToT->GetArray(i);
 419     TArrayF & refaToT = * parrToT;
 420     fArrayTime->AddArray(i, new TArrayF(fgkchannel));
 421     TArrayF * parrTime = fArrayToT->GetArray(i);
 422     TArrayF & refaTime = * parrTime;
 423     for (Int_t j = 0;j<AliTOFcalib::fgkchannel;j++){
 424       refaToT[j]=0.;      //ToT[i][j]=j;
 425       refaTime[j]=0.;      //Time[i][j]=j;
 426     }
 427   }
 428
 429   for (Int_t i=0; i< nelements; i++) {
 430     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
 431     Int_t ipid = element->GetCombID();
 432     Double_t etime = 0;   //expected time
 433     Double_t expTime[10];
 434     element->GetIntegratedTimes(expTime);
 435     if (ipid == 0) etime = expTime[2]*1E-3; //ns
 436     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
 437     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
 438     else AliError("No pid from combinatorial algo for this track");
 439     Double_t mtime = (Double_t)element->GetTOFsignal()*1E-3;  //measured time
 440     Double_t mToT = (Double_t) element->GetToT();  //measured ToT, ns
 441     //select the correspondent channel with its simulated ToT spectrum
 442     //summing up everything, index = 0 for all channels:
 443     Int_t index = element->GetTOFCalChannel();
 444     Int_t index2 = number[index];
 445     TArrayF * parrToT = fArrayToT->GetArray(index);
 446     TArrayF & refaToT = * parrToT;
 447     refaToT[index2] = (Float_t)mToT;
 448     TArrayF * parrTime = fArrayTime->GetArray(index);
 449     TArrayF & refaTime = * parrTime;
 450     refaTime[index2] = (Float_t)(mtime-etime);
 451     number[index]++;
 452   }
 453
 454   for (Int_t i=0;i<1;i++){
 455     TH1F * hProf = Profile(i);
 456     TF1* fGold = SetFitFunctions(hProf);
 457     Int_t nfpar = fGold->GetNpar();
 458     Float_t par[6];
 459     for(Int_t kk=0;kk<6;kk++){
 460       par[kk]=0;
 461     }
 462     for (Int_t kk = 0; kk< nfpar; kk++){
 463       par[kk]=fGold->GetParameter(kk);
 464     }
 465     AliTOFChannel * CalChannel = fTOFCal->GetChannel(i);
 466     CalChannel->SetSlewPar(par);
 467   }
 468   delete[] number;
 469 }
 470
 471 //___________________________________________________________________________
 472
 473 TH1F* AliTOFcalib::Profile(Int_t ich){
 474   const Int_t nbinToT = 650;
 475   Int_t nbinTime = 400;
 476   Float_t minTime = -10.5; //ns
 477   Float_t maxTime = 10.5; //ns
 478   Float_t minToT = 7.5; //ns
 479   Float_t maxToT = 40.; //ns
 480   Float_t DeltaToT = (maxToT-minToT)/nbinToT;
 481   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];
 482   Int_t n[nbinToT+1], nentrx[nbinToT+1];
 483   Double_t sigmaToT[nbinToT+1];
 484   for (Int_t i = 0; i < nbinToT+1 ; i++){
 485     mTime[i]=0;
 486     mToT[i]=0;
 487     n[i]=0;
 488     meanTime[i]=0;
 489     meanTime2[i]=0;
 490     ToT[i]=0;
 491     ToT2[i]=0;
 492     meanToT[i]=0;
 493     meanToT2[i]=0;
 494     Time[i]=0;
 495     Time2[i]=0;
 496     xlow[i]=0;
 497     sigmaTime[i]=0;
 498     sigmaToT[i]=0;
 499     n[i]=0;
 500     nentrx[i]=0;
 501   }
 502   TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", nbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
 503   TArrayF * parrToT = fArrayToT->GetArray(ich);
 504   TArrayF & refaToT = * parrToT;
 505   TArrayF * parrTime = fArrayTime->GetArray(ich);
 506   TArrayF & refaTime = * parrTime;
 507   for (Int_t j = 0; j < AliTOFcalib::fgkchannel; j++){
 508     if (refaToT[j] == 0) continue;
 509     Int_t nx = (Int_t)((refaToT[j]-minToT)/DeltaToT)+1;
 510     if ((refaToT[j] != 0) && (refaTime[j] != 0)){
 511       Time[nx]+=refaTime[j];
 512       Time2[nx]+=(refaTime[j])*(refaTime[j]);
 513       ToT[nx]+=refaToT[j];
 514       ToT2[nx]+=refaToT[j]*refaToT[j];
 515       nentrx[nx]++;
 516       hSlewing->Fill(refaToT[j],refaTime[j]);
 517     }
 518   }
 519   Int_t nbinsToT=hSlewing->GetNbinsX();
 520   if (nbinsToT != nbinToT) {
 521     AliError("Profile :: incompatible numbers of bins");
 522     return 0x0;
 523   }
 524
 525   Int_t usefulBins=0;
 526   TH1F *histo = new TH1F("histo", "1D Time vs ToT", nbinsToT, minToT, maxToT);
 527   for (Int_t i=1;i<=nbinsToT;i++){
 528     if (nentrx[i]!=0){
 529     n[usefulBins]+=nentrx[i];
 530     if (n[usefulBins]==0 && i == nbinsToT) {
 531       break;
 532     }
 533     meanTime[usefulBins]+=Time[i];
 534     meanTime2[usefulBins]+=Time2[i];
 535     meanToT[usefulBins]+=ToT[i];
 536     meanToT2[usefulBins]+=ToT2[i];
 537     if (n[usefulBins]<20 && i!=nbinsToT) continue;
 538     mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
 539     mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
 540     sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
 541                                    *(meanTime2[usefulBins]-meanTime[usefulBins]
 542                                      *meanTime[usefulBins]/n[usefulBins]));
 543     if ((1./n[usefulBins]/n[usefulBins]
 544          *(meanToT2[usefulBins]-meanToT[usefulBins]
 545            *meanToT[usefulBins]/n[usefulBins]))< 0) {
 546       AliError(" too small radical" );
 547       sigmaToT[usefulBins]=0;
 548     }
 549     else{
 550       sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
 551                                      *(meanToT2[usefulBins]-meanToT[usefulBins]
 552                                        *meanToT[usefulBins]/n[usefulBins]));
 553     }
 554     usefulBins++;
 555     }
 556   }
 557   for (Int_t i=0;i<usefulBins;i++){
 558     Int_t binN = (Int_t)((mToT[i]-minToT)/DeltaToT)+1;
 559     histo->Fill(mToT[i],mTime[i]);
 560     histo->SetBinError(binN,sigmaTime[i]);
 561   }
 562   return histo;
 563 }
 564 //_____________________________________________________________________________
 565
 566 void AliTOFcalib::CorrectESDTime(){
 567   Int_t nelements = fESDsel->GetEntries();
 568   for (Int_t i=0; i< nelements; i++) {
 569     AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
 570     Int_t index = element->GetTOFCalChannel();
 571     Float_t ToT = element->GetToT();
 572     //select the correspondent channel with its simulated ToT spectrum
 573     //summing up everything, index = 0 for all channels:
 574     Int_t ipid = element->GetCombID();
 575     Double_t etime = 0;   //expected time
 576     Double_t expTime[10];
 577     element->GetIntegratedTimes(expTime);
 578     if (ipid == 0) etime = expTime[2]*1E-3; //ns
 579     else if (ipid == 1) etime = expTime[3]*1E-3; //ns
 580     else if (ipid == 2) etime = expTime[4]*1E-3; //ns
 581     Float_t par[6];
 582     AliTOFChannel * CalChannel = fTOFCal->GetChannel(index);
 583     for (Int_t j = 0; j<6; j++){
 584       par[j]=CalChannel->GetSlewPar(j);
 585     }
 586     Float_t TimeCorr=0;
 587     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;
 588   }
 589 }
 590 //_____________________________________________________________________________
 591
 592 void AliTOFcalib::CorrectESDTime(AliESD *event){
 593
 594   Int_t ntrk =0;
 595   ntrk=event->GetNumberOfTracks();
 596   for (Int_t itrk=0; itrk<ntrk; itrk++) {
 597     AliESDtrack *t=event->GetTrack(itrk);
 598     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
 599       continue;
 600     }
 601     //IsStartedTimeIntegral
 602     if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
 603       continue;
 604     }
 605     UInt_t AssignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
 606     if(AssignedTOFcluster==0){ // not matched
 607       continue;
 608     }
 609     Int_t index = t->GetTOFCalChannel();
 610     AliTOFChannel * CalChannel = fTOFCal->GetChannel(index);
 611     Float_t par[6];
 612     for (Int_t j = 0; j<6; j++){
 613       par[j]=CalChannel->GetSlewPar(j);
 614     }
 615     Float_t ToT = t->GetTOFsignalToT();
 616     Float_t TimeCorr =0;
 617     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;
 618   }
 619 }
 620 //_____________________________________________________________________________
 621
 622 void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun){
 623   AliCDBManager *man = AliCDBManager::Instance();
 624   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 625   Char_t *sel1 = "Par" ;
 626   Char_t  out[100];
 627   sprintf(out,"%s/%s",sel,sel1);
 628   AliCDBId id(out,minrun,maxrun);
 629   AliCDBMetaData *md = new AliCDBMetaData();
 630   md->SetResponsible("Chiara Zampolli");
 631   man->Put(fTOFCal,id,md);
 632 }
 633 //_____________________________________________________________________________
 634
 635 void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal){
 636   AliCDBManager *man = AliCDBManager::Instance();
 637   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 638   Char_t *sel1 = "Par" ;
 639   Char_t  out[100];
 640   sprintf(out,"%s/%s",sel,sel1);
 641   AliCDBId id(out,minrun,maxrun);
 642   AliCDBMetaData *md = new AliCDBMetaData();
 643   md->SetResponsible("Chiara Zampolli");
 644   man->Put(cal,id,md);
 645 }
 646 //_____________________________________________________________________________
 647
 648 void AliTOFcalib::ReadParFromCDB(Char_t *sel, Int_t nrun){
 649   AliCDBManager *man = AliCDBManager::Instance();
 650   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 651   Char_t *sel1 = "Par" ;
 652   Char_t  out[100];
 653   sprintf(out,"%s/%s",sel,sel1);
 654   AliCDBEntry *entry = man->Get(out,nrun);
 655   AliTOFCal *cal =(AliTOFCal*)entry->GetObject();
 656   fTOFCal = cal;
 657 }
 658 //_____________________________________________________________________________
 659 void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun){
 660
 661
 662   //for the time being, only one spectrum is used
 663   TFile *spFile = new TFile("$ALICE_ROOT/TOF/data/spectrum.root","read");
 664   TH1F * hToT;
 665   // Retrieve ToT Spectrum
 666   spFile->GetObject("ToT",hToT);
 667
 668   fTOFSimToT=hToT;
 669
 670   // Retrieve Time over TOT dependence
 671
 672   TH1F * h = (TH1F*)spFile->Get("TimeToTFit");
 673   TList * list = (TList*)h->GetListOfFunctions();
 674   TF1* f = (TF1*)list->At(0);
 675   Float_t par[6] = {0,0,0,0,0,0};
 676   Int_t npar=f->GetNpar();
 677   for (Int_t ipar=0;ipar<npar;ipar++){
 678     par[ipar]=f->GetParameter(ipar);
 679   }
 680
 681   for(Int_t iTOFch=0; iTOFch<fTOFSimCal->NPads();iTOFch++){
 682     AliTOFChannel * CalChannel = fTOFSimCal->GetChannel(iTOFch);
 683     CalChannel->SetSlewPar(par);
 684   }
 685
 686   // Store them in the CDB
 687
 688   AliCDBManager *man = AliCDBManager::Instance();
 689   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 690   AliCDBMetaData *md = new AliCDBMetaData();
 691   md->SetResponsible("Chiara Zampolli");
 692   Char_t *sel1 = "SimPar" ;
 693   Char_t  out[100];
 694   sprintf(out,"%s/%s",sel,sel1);
 695   AliCDBId id1(out,minrun,maxrun);
 696   man->Put(fTOFSimCal,id1,md);
 697   Char_t *sel2 = "SimHisto" ;
 698   sprintf(out,"%s/%s",sel,sel2);
 699   AliCDBId id2(out,minrun,maxrun);
 700   man->Put(fTOFSimToT,id2,md);
 701 }
 702
 703 //_____________________________________________________________________________
 704 void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal, TH1F * histo){
 705
 706   // Retrieve ToT Spectrum
 707   fTOFSimToT=histo;
 708   fTOFSimCal=cal;
 709   AliCDBManager *man = AliCDBManager::Instance();
 710   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 711   AliCDBMetaData *md = new AliCDBMetaData();
 712   md->SetResponsible("Chiara Zampolli");
 713   Char_t *sel1 = "SimPar" ;
 714   Char_t  out[100];
 715   sprintf(out,"%s/%s",sel,sel1);
 716   AliCDBId id1(out,minrun,maxrun);
 717   man->Put(fTOFSimCal,id1,md);
 718   Char_t *sel2 = "SimHisto" ;
 719   sprintf(out,"%s/%s",sel,sel2);
 720   AliCDBId id2(out,minrun,maxrun);
 721   man->Put(fTOFSimToT,id2,md);
 722 }
 723 //_____________________________________________________________________________
 724 void AliTOFcalib::ReadSimParFromCDB(Char_t *sel, Int_t nrun){
 725   AliCDBManager *man = AliCDBManager::Instance();
 726   if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
 727   AliCDBMetaData *md = new AliCDBMetaData();
 728   md->SetResponsible("Chiara Zampolli");
 729   Char_t *sel1 = "SimPar" ;
 730   Char_t  out[100];
 731   sprintf(out,"%s/%s",sel,sel1);
 732   AliCDBEntry *entry1 = man->Get(out,nrun);
 733   AliTOFCal *cal =(AliTOFCal*)entry1->GetObject();
 734   fTOFSimCal=cal;
 735   Char_t *sel2 = "SimHisto" ;
 736   sprintf(out,"%s/%s",sel,sel2);
 737   AliCDBEntry *entry2 = man->Get(out,nrun);
 738   TH1F *histo =(TH1F*)entry2->GetObject();
 739   fTOFSimToT=histo;
 740 }
 741 //_____________________________________________________________________________
 742
 743 Int_t AliTOFcalib::GetIndex(Int_t *detId){
 744   Int_t isector = detId[0];
 745   if (isector >= fNSector)
 746     AliError(Form("Wrong sector number in TOF (%d) !",isector));
 747   Int_t iplate = detId[1];
 748   if (iplate >= fNPlate)
 749     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 750   Int_t istrip = detId[2];
 751   Int_t ipadz = detId[3];
 752   Int_t ipadx = detId[4];
 753   Int_t stripOffset = 0;
 754   switch (iplate) {
 755   case 0:
 756     stripOffset = 0;
 757     break;
 758   case 1:
 759     stripOffset = fNStripC;
 760     break;
 761   case 2:
 762     stripOffset = fNStripC+fNStripB;
 763     break;
 764   case 3:
 765     stripOffset = fNStripC+fNStripB+fNStripA;
 766     break;
 767   case 4:
 768     stripOffset = fNStripC+fNStripB+fNStripA+fNStripB;
 769     break;
 770   default:
 771     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 772     break;
 773   };
 774
 775   Int_t idet = ((2*(fNStripC+fNStripB)+fNStripA)*fNpadZ*fNpadX)*isector +
 776                (stripOffset*fNpadZ*fNpadX)+
 777                (fNpadZ*fNpadX)*istrip+
 778                (fNpadX)*ipadz+
 779                 ipadx;
 780   return idet;
 781 }
 782