]> 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
just Fixing Log info
[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
16/*$Log$
17author: 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
55extern TROOT *gROOT;
56extern TStyle *gStyle;
57
58ClassImp(AliTOFcalib)
59
60const Int_t AliTOFcalib::fgkchannel = 5000;
61const 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//_______________________________________________________________________
63AliTOFcalib::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
84AliTOFcalib::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
123AliTOFcalib::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
146AliTOFcalib::~AliTOFcalib()
147{
148 delete fArrayToT;
149 delete fArrayTime;
150 delete flistFunc;
151 delete[] fhToT;
152 delete fTOFCal;
153 delete fESDsel;
154}
155//____________________________________________________________________________
156
157void 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
174void 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
209void 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
231TF1* 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
288TClonesArray* 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
361void 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
444void 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
606void 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
678TH1F* 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
771void 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
804void 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
864void 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
873void 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
888Int_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