]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/AliTOFcalib.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix in treatment of the vertex finder covariance matrix (Andrea)
[u/mrichter/AliRoot.git] / TOF / AliTOFcalib.cxx
CommitLineData
6dc9348d 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
762446e0 16/*
17$Log$
7a4f634e 18Revision 1.12 2007/03/23 11:31:16 arcelli
19CDB Entry for TOF Reconstruction Parameters
20
7037bd93 21Revision 1.11 2007/02/28 18:08:26 arcelli
22Add protection against failed retrieval of the CDB cal object
23
1b20c168 24Revision 1.10 2006/08/22 13:30:49 arcelli
25removal of effective c++ warnings (C.Zampolli)
26
655e379f 27Revision 1.9 2006/04/20 22:30:50 hristov
28Coding conventions (Annalisa)
29
0e46b9ae 30Revision 1.8 2006/04/16 22:29:05 hristov
31Coding conventions (Annalisa)
32
7aeeaf38 33Revision 1.7 2006/04/16 20:12:46 hristov
34Removing memory leak in case of cached CDB entries
35
85fc78e3 36Revision 1.6 2006/04/11 15:28:32 hristov
37Checks on cache status before deleting calibration objects (A.Colla)
38
5eaae242 39Revision 1.5 2006/04/05 08:35:38 hristov
40Coding conventions (S.Arcelli, C.Zampolli)
41
340693af 42Revision 1.4 2006/03/31 11:26:46 arcelli
43 changing CDB Ids according to standard convention
44
28dd10b6 45Revision 1.3 2006/03/28 14:57:02 arcelli
46updates to handle new V5 geometry & some re-arrangements
47
d4ad0d6b 48Revision 1.2 2006/02/13 17:22:26 arcelli
49just Fixing Log info
50
762446e0 51Revision 1.1 2006/02/13 16:10:48 arcelli
52Add classes for TOF Calibration (C.Zampolli)
53
6dc9348d 54author: Chiara Zampolli, zampolli@bo.infn.it
762446e0 55*/
6dc9348d 56
57///////////////////////////////////////////////////////////////////////////////
58// //
59// class for TOF calibration //
60// //
61///////////////////////////////////////////////////////////////////////////////
62
6dc9348d 63#include "TF1.h"
0e46b9ae 64#include "TFile.h"
6dc9348d 65#include "TH1F.h"
66#include "TH2F.h"
6dc9348d 67#include "TList.h"
0e46b9ae 68#include "TROOT.h"
69#include "TStyle.h"
70
71#include "AliCDBEntry.h"
72#include "AliCDBId.h"
6dc9348d 73#include "AliCDBManager.h"
74#include "AliCDBMetaData.h"
0e46b9ae 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"
7037bd93 85#include "AliTOFRecoParam.h"
6dc9348d 86
87extern TROOT *gROOT;
88extern TStyle *gStyle;
89
90ClassImp(AliTOFcalib)
91
92const Int_t AliTOFcalib::fgkchannel = 5000;
6dc9348d 93//_______________________________________________________________________
655e379f 94AliTOFcalib::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{
340693af 112 //TOF Calibration Class ctor
6dc9348d 113 fArrayToT = 0x0;
114 fArrayTime = 0x0;
6dc9348d 115 fESDsel = 0x0;
d4ad0d6b 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
85fc78e3 126 fTOFCal = 0x0;
127 fTOFSimCal = 0x0;
128 fTOFSimToT = 0x0;
d4ad0d6b 129 delete geom;
6dc9348d 130}
d4ad0d6b 131//_______________________________________________________________________
655e379f 132AliTOFcalib::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{
340693af 150 //TOF Calibration Class ctor, taking the TOF geometry as input
6dc9348d 151 fArrayToT = 0x0;
152 fArrayTime = 0x0;
6dc9348d 153 fESDsel = 0x0;
d4ad0d6b 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
85fc78e3 162 fTOFCal = 0x0;
163 fTOFSimCal = 0x0;
164 fTOFSimToT = 0x0;
6dc9348d 165}
166//____________________________________________________________________________
167
655e379f 168AliTOFcalib::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)
6dc9348d 185{
340693af 186 //TOF Calibration Class copy ctor
6dc9348d 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;
d4ad0d6b 194 fNChannels = calib.fNChannels;
6dc9348d 195 fArrayToT = calib.fArrayToT;
196 fArrayTime = calib.fArrayTime;
6dc9348d 197 fTOFCal=calib.fTOFCal;
d4ad0d6b 198 fTOFSimCal = calib.fTOFSimCal;
340693af 199 fTOFSimToT=calib.fTOFSimToT;
6dc9348d 200}
201
202//____________________________________________________________________________
7aeeaf38 203
204AliTOFcalib& 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//____________________________________________________________________________
6dc9348d 224
225AliTOFcalib::~AliTOFcalib()
226{
340693af 227 //TOF Calibration Class dtor
6dc9348d 228 delete fArrayToT;
229 delete fArrayTime;
5eaae242 230
231 if(!(AliCDBManager::Instance()->GetCacheFlag())){ // CDB objects must NOT be deleted if cache is active!
232 delete fTOFCal;
233 delete fTOFSimCal;
234 }
235
6dc9348d 236 delete fESDsel;
237}
6dc9348d 238//__________________________________________________________________________
239
340693af 240TF1* AliTOFcalib::SetFitFunctions(TH1F *histo)
241{
242 //Define Fit Functions for Slewing Correction
6dc9348d 243 TF1 * fpol[3];
340693af 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;
6dc9348d 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};
340693af 254 Double_t nchi[3]={1E6,1E6,1E6};
6dc9348d 255 Double_t bestchi=1E6;
256 TF1 * fGold=0x0;
257 Int_t nonzero =0;
258 Int_t numberOfpar =0;
340693af 259 for (Int_t j=0; j<knbins; j++){
6dc9348d 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();
340693af 286 nchi[ifun] = (Double_t)chi[ifun]/ndf[ifun];
287 if (nchi[ifun]<bestchi) {
288 bestchi=nchi[ifun];
6dc9348d 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
d4ad0d6b 299void AliTOFcalib::SelectESD(AliESD *event)
6dc9348d 300{
340693af 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
6dc9348d 304 Int_t ntrk =0;
305 Int_t ngoodtrkfinalToT = 0;
306 ntrk=event->GetNumberOfTracks();
307 fESDsel = new TObjArray(ntrk);
308 fESDsel->SetOwner();
340693af 309 TObjArray uCdatatemp(ntrk);
6dc9348d 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();
340693af 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
6dc9348d 330 continue;
331 }
6dc9348d 332 AliTOFcalibESD *unc = new AliTOFcalibESD;
333 unc->CopyFromAliESD(t);
6dc9348d 334 Double_t c1[15];
335 unc->GetExternalCovariance(c1);
340693af 336 uCdatatemp.Add(unc);
6dc9348d 337 ngoodtrk++;
338 }
339 for (Int_t i = 0; i < ngoodtrk ; i ++){
340693af 340 AliTOFcalibESD *unc = (AliTOFcalibESD*)uCdatatemp.At(i);
6dc9348d 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
340693af 351void AliTOFcalib::CombESDId()
352{
353 //track PID for calibration
6dc9348d 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.};
340693af 367 Float_t timeResolution = 0.90e-10; // 90 ps by default
368 Float_t timeresolutioninns=timeResolution*(1.e+9); // convert in [ns]
6dc9348d 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.;
340693af 402 Float_t et0best=999.;
6dc9348d 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.;
340693af 426 Float_t emeantzero=0.;
6dc9348d 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]
340693af 449 emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
6dc9348d 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;
340693af 486 et0best=emeantzero;
6dc9348d 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
515void AliTOFcalib::CalibrateESD(){
340693af 516 //Calibrate selected ESD times
6dc9348d 517 Int_t nelements = fESDsel->GetEntries();
d4ad0d6b 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++){
6dc9348d 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
6dc9348d 547 //select the correspondent channel with its simulated ToT spectrum
548 //summing up everything, index = 0 for all channels:
d4ad0d6b 549 Int_t index = element->GetTOFCalChannel();
6dc9348d 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 }
85fc78e3 571 if (!fTOFCal) {
572 AliTOFGeometry *geom=new AliTOFGeometryV5();
573 fTOFCal = new AliTOFCal(geom);
574 fTOFCal->CreateArray();
575 delete geom;
576 }
340693af 577 AliTOFChannel * calChannel = fTOFCal->GetChannel(i);
578 calChannel->SetSlewPar(par);
6dc9348d 579 }
580 delete[] number;
581}
582
583//___________________________________________________________________________
584
340693af 585TH1F* AliTOFcalib::Profile(Int_t ich)
586{
587 //Prepare histograms for Slewing Correction
588 const Int_t knbinToT = 650;
6dc9348d 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
340693af 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++){
6dc9348d 599 mTime[i]=0;
600 mToT[i]=0;
601 n[i]=0;
602 meanTime[i]=0;
603 meanTime2[i]=0;
340693af 604 vToT[i]=0;
605 vToT2[i]=0;
6dc9348d 606 meanToT[i]=0;
607 meanToT2[i]=0;
340693af 608 vTime[i]=0;
609 vTime2[i]=0;
6dc9348d 610 xlow[i]=0;
611 sigmaTime[i]=0;
612 sigmaToT[i]=0;
613 n[i]=0;
614 nentrx[i]=0;
615 }
340693af 616 TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", knbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
6dc9348d 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;
340693af 623 Int_t nx = (Int_t)((refaToT[j]-minToT)/deltaToT)+1;
6dc9348d 624 if ((refaToT[j] != 0) && (refaTime[j] != 0)){
340693af 625 vTime[nx]+=refaTime[j];
626 vTime2[nx]+=(refaTime[j])*(refaTime[j]);
627 vToT[nx]+=refaToT[j];
628 vToT2[nx]+=refaToT[j]*refaToT[j];
6dc9348d 629 nentrx[nx]++;
630 hSlewing->Fill(refaToT[j],refaTime[j]);
631 }
632 }
633 Int_t nbinsToT=hSlewing->GetNbinsX();
340693af 634 if (nbinsToT != knbinToT) {
6dc9348d 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 }
340693af 647 meanTime[usefulBins]+=vTime[i];
648 meanTime2[usefulBins]+=vTime2[i];
649 meanToT[usefulBins]+=vToT[i];
650 meanToT2[usefulBins]+=vToT2[i];
6dc9348d 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++){
340693af 672 Int_t binN = (Int_t)((mToT[i]-minToT)/deltaToT)+1;
6dc9348d 673 histo->Fill(mToT[i],mTime[i]);
674 histo->SetBinError(binN,sigmaTime[i]);
675 }
676 return histo;
677}
678//_____________________________________________________________________________
679
340693af 680void AliTOFcalib::CorrectESDTime()
681{
682 //Calculate the corrected TOF time
6dc9348d 683 Int_t nelements = fESDsel->GetEntries();
684 for (Int_t i=0; i< nelements; i++) {
685 AliTOFcalibESD *element=(AliTOFcalibESD*)fESDsel->At(i);
d4ad0d6b 686 Int_t index = element->GetTOFCalChannel();
340693af 687 Float_t tToT = element->GetToT();
6dc9348d 688 //select the correspondent channel with its simulated ToT spectrum
689 //summing up everything, index = 0 for all channels:
6dc9348d 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];
85fc78e3 698 if (!fTOFCal) {
699 AliTOFGeometry *geom=new AliTOFGeometryV5();
700 fTOFCal = new AliTOFCal(geom);
701 fTOFCal->CreateArray();
702 delete geom;
703 }
340693af 704 AliTOFChannel * calChannel = fTOFCal->GetChannel(index);
6dc9348d 705 for (Int_t j = 0; j<6; j++){
340693af 706 par[j]=calChannel->GetSlewPar(j);
6dc9348d 707 }
340693af 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;
6dc9348d 710 }
711}
712//_____________________________________________________________________________
713
d4ad0d6b 714void AliTOFcalib::CorrectESDTime(AliESD *event){
340693af 715 //Calculate the corrected TOF time
d4ad0d6b 716
6dc9348d 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 }
340693af 728 UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
729 if(assignedTOFcluster==0){ // not matched
6dc9348d 730 continue;
731 }
d4ad0d6b 732 Int_t index = t->GetTOFCalChannel();
85fc78e3 733 if (!fTOFCal) {
734 AliTOFGeometry *geom=new AliTOFGeometryV5();
735 fTOFCal = new AliTOFCal(geom);
736 fTOFCal->CreateArray();
737 delete geom;
738 }
340693af 739 AliTOFChannel * calChannel = fTOFCal->GetChannel(index);
6dc9348d 740 Float_t par[6];
741 for (Int_t j = 0; j<6; j++){
340693af 742 par[j]=calChannel->GetSlewPar(j);
6dc9348d 743 }
340693af 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;
6dc9348d 747 }
748}
749//_____________________________________________________________________________
750
340693af 751void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun)
752{
753 //Write calibration parameters to the CDB
6dc9348d 754 AliCDBManager *man = AliCDBManager::Instance();
d4ad0d6b 755 if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
28dd10b6 756 Char_t *sel1 = "Par" ;
d4ad0d6b 757 Char_t out[100];
758 sprintf(out,"%s/%s",sel,sel1);
759 AliCDBId id(out,minrun,maxrun);
6dc9348d 760 AliCDBMetaData *md = new AliCDBMetaData();
d4ad0d6b 761 md->SetResponsible("Chiara Zampolli");
85fc78e3 762 if (!fTOFCal) {
763 AliTOFGeometry *geom=new AliTOFGeometryV5();
764 fTOFCal = new AliTOFCal(geom);
765 fTOFCal->CreateArray();
766 delete geom;
767 }
6dc9348d 768 man->Put(fTOFCal,id,md);
5eaae242 769 delete md;
6dc9348d 770}
771//_____________________________________________________________________________
772
d4ad0d6b 773void AliTOFcalib::WriteParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal){
340693af 774 //Write calibration parameters to the CDB
6dc9348d 775 AliCDBManager *man = AliCDBManager::Instance();
d4ad0d6b 776 if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
28dd10b6 777 Char_t *sel1 = "Par" ;
d4ad0d6b 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);
5eaae242 784 delete md;
d4ad0d6b 785}
786//_____________________________________________________________________________
787
1b20c168 788Bool_t AliTOFcalib::ReadParFromCDB(Char_t *sel, Int_t nrun)
340693af 789{
790 //Read calibration parameters from the CDB
d4ad0d6b 791 AliCDBManager *man = AliCDBManager::Instance();
792 if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
28dd10b6 793 Char_t *sel1 = "Par" ;
d4ad0d6b 794 Char_t out[100];
795 sprintf(out,"%s/%s",sel,sel1);
7a4f634e 796 if (!man->Get(out,nrun)) {
797 return kFALSE;
1b20c168 798 }
7a4f634e 799 AliCDBEntry *entry = man->Get(out,nrun);
800 if(!entry->GetObject()){
801 return kFALSE;
802 }
803
6dc9348d 804 AliTOFCal *cal =(AliTOFCal*)entry->GetObject();
805 fTOFCal = cal;
1b20c168 806 return kTRUE;
7a4f634e 807
6dc9348d 808}
809//_____________________________________________________________________________
340693af 810void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun)
811{
812 //Write Sim miscalibration parameters to the CDB
d4ad0d6b 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 }
85fc78e3 833 if (!fTOFSimCal) {
834 AliTOFGeometry *geom=new AliTOFGeometryV5();
835 fTOFSimCal = new AliTOFCal(geom);
836 fTOFSimCal->CreateArray();
837 delete geom;
838 }
d4ad0d6b 839 for(Int_t iTOFch=0; iTOFch<fTOFSimCal->NPads();iTOFch++){
340693af 840 AliTOFChannel * calChannel = fTOFSimCal->GetChannel(iTOFch);
841 calChannel->SetSlewPar(par);
d4ad0d6b 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");
28dd10b6 850 Char_t *sel1 = "SimPar" ;
d4ad0d6b 851 Char_t out[100];
852 sprintf(out,"%s/%s",sel,sel1);
853 AliCDBId id1(out,minrun,maxrun);
854 man->Put(fTOFSimCal,id1,md);
28dd10b6 855 Char_t *sel2 = "SimHisto" ;
d4ad0d6b 856 sprintf(out,"%s/%s",sel,sel2);
857 AliCDBId id2(out,minrun,maxrun);
858 man->Put(fTOFSimToT,id2,md);
5eaae242 859 delete md;
d4ad0d6b 860}
861
862//_____________________________________________________________________________
863void AliTOFcalib::WriteSimParOnCDB(Char_t *sel, Int_t minrun, Int_t maxrun, AliTOFCal *cal, TH1F * histo){
340693af 864 //Write Sim miscalibration parameters to the CDB
d4ad0d6b 865
d4ad0d6b 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");
28dd10b6 872 Char_t *sel1 = "SimPar" ;
d4ad0d6b 873 Char_t out[100];
874 sprintf(out,"%s/%s",sel,sel1);
875 AliCDBId id1(out,minrun,maxrun);
876 man->Put(fTOFSimCal,id1,md);
28dd10b6 877 Char_t *sel2 = "SimHisto" ;
d4ad0d6b 878 sprintf(out,"%s/%s",sel,sel2);
879 AliCDBId id2(out,minrun,maxrun);
880 man->Put(fTOFSimToT,id2,md);
5eaae242 881 delete md;
d4ad0d6b 882}
883//_____________________________________________________________________________
340693af 884void AliTOFcalib::ReadSimParFromCDB(Char_t *sel, Int_t nrun)
885{
886 //Read miscalibration parameters from the CDB
d4ad0d6b 887 AliCDBManager *man = AliCDBManager::Instance();
888 if(!man->IsDefaultStorageSet())man->SetDefaultStorage("local://$ALICE_ROOT");
7a4f634e 889
890 // The Slewing Pars
891
28dd10b6 892 Char_t *sel1 = "SimPar" ;
d4ad0d6b 893 Char_t out[100];
894 sprintf(out,"%s/%s",sel,sel1);
7a4f634e 895 if (!man->Get(out,nrun)) {
896 AliFatal("Exiting, no CDB object (SimPar) found!!!");
897 exit(0);
898 }
d4ad0d6b 899 AliCDBEntry *entry1 = man->Get(out,nrun);
7a4f634e 900 if(!entry1->GetObject()){
901 AliFatal("Exiting, no CDB object (SimPar) found!!!");
902 exit(0);
903 }
d4ad0d6b 904 AliTOFCal *cal =(AliTOFCal*)entry1->GetObject();
905 fTOFSimCal=cal;
7a4f634e 906
907 // The Tot Histo
908
28dd10b6 909 Char_t *sel2 = "SimHisto" ;
d4ad0d6b 910 sprintf(out,"%s/%s",sel,sel2);
7a4f634e 911 if (!man->Get(out,nrun)) {
912 AliFatal("Exiting, no CDB object (SimHisto) found!!!");
913 exit(0);
914 }
d4ad0d6b 915 AliCDBEntry *entry2 = man->Get(out,nrun);
7a4f634e 916 if(!entry2->GetObject()){
917 AliFatal("Exiting, no CDB object (SimHisto) found!!!");
918 exit(0);
919 }
d4ad0d6b 920 TH1F *histo =(TH1F*)entry2->GetObject();
921 fTOFSimToT=histo;
922}
923//_____________________________________________________________________________
6dc9348d 924
7037bd93 925//_____________________________________________________________________________
926void 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//_____________________________________________________________________________
941AliTOFRecoParam * 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);
7a4f634e 949 if (!man->Get(out,nrun)) {
950 AliFatal("Exiting, no CDB object (RecPar) found!!!");
951 exit(0);
952 }
7037bd93 953 AliCDBEntry *entry = man->Get(out,nrun);
7a4f634e 954 if(!entry->GetObject()){
955 AliFatal("Exiting, no CDB object (RecPar) found!!!");
956 exit(0);
957 }
958
7037bd93 959 AliTOFRecoParam *param=(AliTOFRecoParam*)entry->GetObject();
960 return param;
961}
962//_____________________________________________________________________________
963
340693af 964Int_t AliTOFcalib::GetIndex(Int_t *detId)
965{
966 //Retrieve calibration channel index
6dc9348d 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