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