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