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