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