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