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