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