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