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