]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TOF/AliTOFCalibTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Memory leak corrected (Ch. Klein-Boesing)
[u/mrichter/AliRoot.git] / TOF / AliTOFCalibTask.cxx
CommitLineData
1ed77580 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/*
45af7330 16$Log$
1ed77580 17*/
18
19// task to perform TOF calibration
20// optimized for pp runs
21// expect a maximum of 100 entries per channel
22// different ways to calibrate are implemented
23// C. Zampolli
24
25#include <TChain.h>
26#include <TObject.h>
27#include <TCanvas.h>
28#include <TROOT.h>
29#include <TSystem.h>
30#include <TProfile.h>
31#include <TF1.h>
32#include <TTree.h>
33#include <TFile.h>
34#include <TH1F.h>
35#include <TH1I.h>
36#include <TH1D.h>
37#include <TH2F.h>
38
39#include "AliTOFCalibTask.h"
40#include "AliTOFChannelTask.h"
70ddd0ee 41#include "AliESDEvent.h"
1ed77580 42#include "AliESDtrack.h"
43#include "AliLog.h"
44//______________________________________________________________________________
45AliTOFCalibTask::AliTOFCalibTask(const char *name) :
46 AliAnalysisTask(name,""),
47 fdir(0),
48 fChain(0),
49 fESD(0),
50 fMinEntries(0),
51 fIch(0),
52 fToT(0),
53 fTime(0),
54 fExpTimePi(0),
55 fExpTimeKa(0),
56 fExpTimePr(0),
57 ftree(0x0),
58 fMinTime(0),
59 fbigarray(0x0),
60 findexarray(0x0),
61 fnESD(0),
62 fnESDselected(0),
63 fnESDkTOFout(0),
64 fnESDkTIME(0),
65 fnESDassTOFcl(0),
66 fnESDTIMEcut(0),
67 fnESDTRDcut(0),
68 fhToT(0),
69 fhTime(0),
70 fhExpTimePi(0),
71 fhExpTimeKa(0),
72 fhExpTimePr(0),
73 fhPID(0),
74 fhch(0),
75 fOutputContainer(0)
76{
77 // Constructor.
78 // Input slot #0 works with an Ntuple
79
80 DefineInput(0, TChain::Class());
81 // Output slot #0 writes into a TH1 container
82 DefineOutput(0, TObjArray::Class()) ;
83 fdir = gDirectory->GetPath();
70ddd0ee 84 fbigarray = new Float_t*[TOFCHANNELS];
85 findexarray = new Int_t[TOFCHANNELS];
86
87 for (Int_t i =0;i<TOFCHANNELS;i++){
88 fbigarray[i] = new Float_t[CHENTRIES];
89 findexarray[i]=0;
90 for (Int_t j =0;j<CHENTRIES;j++){
91 fbigarray[i][j]=-1;
92 }
93 }
1ed77580 94}
95
96//______________________________________________________________________________
97AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
98 AliAnalysisTask("AliTOFCalibTask",""),
99 fdir(0),
100 fChain(0),
101 fESD(0),
102 fMinEntries(0),
103 fIch(0),
104 fToT(0),
105 fTime(0),
106 fExpTimePi(0),
107 fExpTimeKa(0),
108 fExpTimePr(0),
109 ftree(0x0),
110 fMinTime(0),
111 fbigarray(0x0),
112 findexarray(0x0),
113 fnESD(0),
114 fnESDselected(0),
115 fnESDkTOFout(0),
116 fnESDkTIME(0),
117 fnESDassTOFcl(0),
118 fnESDTIMEcut(0),
119 fnESDTRDcut(0),
120 fhToT(0),
121 fhTime(0),
122 fhExpTimePi(0),
123 fhExpTimeKa(0),
124 fhExpTimePr(0),
125 fhPID(0),
126 fhch(0),
127 fOutputContainer(0)
128{
129 // Copy Constructor.
130 fdir=calibtask.fdir;
131 fChain=calibtask.fChain;
132 fESD=calibtask.fESD;
133 fMinEntries=calibtask.fMinEntries;
134 fIch=calibtask.fIch;
135 fToT=calibtask.fToT;
136 fTime=calibtask.fTime;
137 fExpTimePi=calibtask.fExpTimePi;
138 fExpTimeKa=calibtask.fExpTimeKa;
139 fExpTimePr=calibtask.fExpTimePr;
140 ftree=calibtask.ftree;
141 fMinTime=calibtask.fMinTime;
1ed77580 142 fnESD=calibtask.fnESD;
143 fnESDselected=calibtask.fnESDselected;
144 fnESDkTOFout=calibtask.fnESDkTOFout;
145 fnESDkTIME=calibtask.fnESDkTIME;
146 fnESDassTOFcl=calibtask.fnESDassTOFcl;
147 fnESDTIMEcut=calibtask.fnESDTIMEcut;
148 fnESDTRDcut=calibtask.fnESDTRDcut;
149 fhToT=calibtask.fhToT;
150 fhTime=calibtask.fhTime;
151 fhExpTimePi=calibtask.fhExpTimePi;
152 fhExpTimeKa=calibtask.fhExpTimeKa;
153 fhExpTimePr=calibtask.fhExpTimePr;
154 fhPID=calibtask.fhPID;
155 fhch=calibtask.fhch;
156 fOutputContainer=calibtask.fOutputContainer;
70ddd0ee 157
158 fbigarray = new Float_t*[TOFCHANNELS];
159 findexarray = new Int_t[TOFCHANNELS];
160
161 for (Int_t i =0;i<TOFCHANNELS;i++){
162 fbigarray[i] = new Float_t[CHENTRIES];
163 findexarray[i]=calibtask.findexarray[i];
164 for (Int_t j =0;j<CHENTRIES;j++){
165 fbigarray[i][j]=calibtask.fbigarray[i][j];
166 }
167 }
168
1ed77580 169}
170//______________________________________________________________________________
171AliTOFCalibTask:: ~AliTOFCalibTask()
172{
173 // destructor
174
175 AliInfo("TOF Calib Task: Deleting");
176 delete[] fbigarray;
177 delete findexarray;
178 delete fOutputContainer;
179 delete fhToT;
180 delete fhTime;
181 delete fhExpTimePi;
182 delete fhExpTimeKa;
183 delete fhExpTimePr;
184 delete fhPID;
185 delete fhch;
186}
187//______________________________________________________________________________
188AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
189{
190 //assignment operator
191 this->fdir=calibtask.fdir;
192 this->fChain=calibtask.fChain;
193 this->fESD=calibtask.fESD;
194 this->fMinEntries=calibtask.fMinEntries;
195 this->fIch=calibtask.fIch;
196 this->fToT=calibtask.fToT;
197 this->fTime=calibtask.fTime;
198 this->fExpTimePi=calibtask.fExpTimePi;
199 this->fExpTimeKa=calibtask.fExpTimeKa;
200 this->fExpTimePr=calibtask.fExpTimePr;
201 this->ftree=calibtask.ftree;
202 this->fMinTime=calibtask.fMinTime;
1ed77580 203 this->fnESD=calibtask.fnESD;
204 this->fnESDselected=calibtask.fnESDselected;
205 this->fnESDkTOFout=calibtask.fnESDkTOFout;
206 this->fnESDkTIME=calibtask.fnESDkTIME;
207 this->fnESDassTOFcl=calibtask.fnESDassTOFcl;
208 this->fnESDTIMEcut=calibtask.fnESDTIMEcut;
209 this->fnESDTRDcut=calibtask.fnESDTRDcut;
210 this->fhToT=calibtask.fhToT;
211 this->fhTime=calibtask.fhTime;
212 this->fhExpTimePi=calibtask.fhExpTimePi;
213 this->fhExpTimeKa=calibtask.fhExpTimeKa;
214 this->fhExpTimePr=calibtask.fhExpTimePr;
215 this->fOutputContainer=calibtask.fOutputContainer;
216 this->fhPID=calibtask.fhPID;
217 this->fhch=calibtask.fhch;
70ddd0ee 218
219 this->fbigarray = new Float_t*[TOFCHANNELS];
220 this->findexarray = new Int_t[TOFCHANNELS];
221 for (Int_t i =0;i<TOFCHANNELS;i++){
222 this->fbigarray[i] = new Float_t[CHENTRIES];
223 this->findexarray[i]=calibtask.findexarray[i];
224 for (Int_t j =0;j<CHENTRIES;j++){
225 this->fbigarray[i][j]=calibtask.fbigarray[i][j];
226 }
227 }
1ed77580 228 return *this;
229}
230//--------------------------------------------------------------------------
231void AliTOFCalibTask::BookHistos(){
232
233 //booking histos
234
235 AliInfo(Form("*** Booking Histograms %s", GetName())) ;
236
237 fhToT=
238 new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
239 fhTime=
240 new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
241 fhExpTimePi=
242 new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
243 fhExpTimeKa=
244 new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
245 fhExpTimePr=
246 new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
247 fhPID=
248 new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
249 fhch=
250 new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
251
252 //create the putput container
253 fOutputContainer = new TObjArray(7) ;
254 fOutputContainer->SetName(GetName()) ;
255
256 fOutputContainer->AddAt(fhToT, 0) ;
257 fOutputContainer->AddAt(fhTime, 1) ;
258 fOutputContainer->AddAt(fhExpTimePi, 2) ;
259 fOutputContainer->AddAt(fhExpTimeKa, 3) ;
260 fOutputContainer->AddAt(fhExpTimePr, 4) ;
261 fOutputContainer->AddAt(fhPID, 5) ;
262 fOutputContainer->AddAt(fhch, 6) ;
263
264}
265//----------------------------------------------------------------------------
266void AliTOFCalibTask::DrawHistos(){
267
268 // drawing output histos
269
270 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
271
272 TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
273 canvasToTTime->Divide(1,2);
274 canvasToTTime->cd(1);
275 fhToT->SetLineColor(4);
276 fhToT->GetXaxis()->SetTitle("ToT (ns)");
277 fhToT->Draw("hist");
278 canvasToTTime->cd(2);
279 fhTime->SetLineColor(4);
280 fhTime->GetXaxis()->SetTitle("Time (ns)");
281 fhTime->Draw("hist");
282 canvasToTTime->Update();
283 canvasToTTime->Print("ToTTime.gif");
284
285 TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
286 canvasExpTime->Divide(1,3);
287 canvasExpTime->cd(1);
288 fhExpTimePi->SetLineColor(4);
289 fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
290 fhExpTimePi->Draw("hist");
291 canvasExpTime->cd(2);
292 fhExpTimeKa->SetLineColor(4);
293 fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
294 fhExpTimeKa->Draw("hist");
295 canvasExpTime->cd(3);
296 fhExpTimePr->SetLineColor(4);
297 fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
298 fhExpTimePr->Draw("hist");
299
300 canvasExpTime->Print("ExpTime.gif");
301
302 TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
303 fhPID->GetXaxis()->SetTitle("Comb PID");
304 fhPID->GetXaxis()->SetBinLabel(1,"#pi");
305 fhPID->GetXaxis()->SetBinLabel(2,"K");
306 fhPID->GetXaxis()->SetBinLabel(3,"p");
307 fhPID->Draw("hist");
308
309 canvasPID->Print("PID.gif");
310
311 TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
312 fhch->GetXaxis()->SetTitle("TOF ch");
313 fhch->Draw("hist");
314 Float_t meanTOFch = 0;
315 for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
316 meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
317 }
318
319 meanTOFch/=TOFCHANNELS;
320 AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
321
322 canvasrndch->Print("rndch.gif");
323
324 char line[1024] ;
325 sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ;
326 gROOT->ProcessLine(line);
327 sprintf(line, ".!rm -fR *.gif");
328 gROOT->ProcessLine(line);
329 AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
330}
331
332//______________________________________________________________________________
333void AliTOFCalibTask::ConnectInputData(const Option_t*)
334{
335 // Initialization of branch container and histograms
336
337 // AliLog::SetClassDebugLevel("AliTOFCalibTask",1);
338 AliInfo(Form("*** Initialization of %s", GetName())) ;
339
340 // Get input data
341 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
342 if (!fChain) {
343 AliError(Form("Input 0 for %s not found\n", GetName()));
344 return ;
345 }
346
347 // One should first check if the branch address was taken by some other task
348 char ** address = (char **)GetBranchAddress(0, "ESD");
349 if (address) {
70ddd0ee 350 fESD = (AliESDEvent*)(*address);
1ed77580 351 } else {
70ddd0ee 352 fESD = new AliESDEvent();
353 AliInfo(" qui ok ");
354 // fESD = (AliESDEvent*)fChain->GetTree()->GetUserInfo()->FindObject("AliESDEvent");
1ed77580 355 fESD->ReadFromTree(fChain) ;
356 }
1ed77580 357
358 BookHistos();
359
360}
70ddd0ee 361//-----------------------------------------------------------------------
362Bool_t AliTOFCalibTask::Notify()
363{
364 // Initialisation of branch container and histograms
365
366 AliInfo(Form("*** We are in %s::Notify()", GetName())) ;
367
368 // Get input data
369 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
370 if (!fChain) {
371 AliError(Form("Input 0 for %s not found\n", GetName()));
372 return kFALSE;
373 }
374
375 // One should first check if the branch address was taken by some other task
376 char ** address = (char **)GetBranchAddress(0, "ESD") ;
377 if (address)
378 fESD = (AliESDEvent *)(*address) ;
379 else {
380 fESD = new AliESDEvent();
381 // fESD = (AliESDEvent*)fChain->GetTree()->GetUserInfo()->FindObject("AliESDEvent");
382 fESD->ReadFromTree(fChain) ;
383 //SetBranchAddress(0,"ESD",&fESD);
384 }
385
386 return kTRUE;
387}
388
1ed77580 389//________________________________________________________________________
390void AliTOFCalibTask::CreateOutputObjects()
391{
392// Create histograms
393}
394
395//______________________________________________________________________________
396void AliTOFCalibTask::Exec(Option_t * opt)
397{
398
399 // main
400
401 AliInfo(Form("*** Executing %s", GetName())) ;
402
403//******* The loop over events -----------------------------------------------
404
405 // Processing of one event
406 Long64_t entry = fChain->GetReadEntry() ;
407 if (!fESD) {
408 AliError("fESD is not connected to the input!") ;
409 return ;
410 }
411
412 if ( !((entry-1)%100) )
413 AliDebug(1,Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
414
415 // ************************ TOF *************************************
416
417 fMinTime=22E3; //ns; not used
418 Int_t ntrk = fESD->GetNumberOfTracks() ;
419 fnESD+=ntrk;
420 Int_t nselected = 0;
421 Int_t itr = -1;
422 while ( ntrk-- ) {
423 itr++;
424 AliESDtrack * t = fESD->GetTrack(ntrk) ;
425 //selecting only good quality tracks
426 if (!Select(t)) continue;
427 nselected++;
428 Int_t ich = Int_t(t->GetTOFCalChannel());
429 fhch->Fill(ich);
430 AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
431 Float_t tot = t->GetTOFsignalToT();
432 Float_t time = t->GetTOFsignalRaw();
70ddd0ee 433 AliDebug(2,Form(" track # %i in schannel %i, time = %f \n",ntrk,ich,time));
1ed77580 434 Double_t expTime[10];
435 t->GetIntegratedTimes(expTime);
436 Float_t expTimePi = expTime[2]*1.E-3;
437 Float_t expTimeKa = expTime[3]*1.E-3;
438 Float_t expTimePr = expTime[4]*1.E-3;
439 if (findexarray[ich]==(Int_t)(CHENTRIES/NIDX)) {
440 AliInfo(Form("too many tracks in channel %i, not storing any more...",ich));
441 continue;
442 }
443 findexarray[ich]++;
444 AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, findexarray[ich] ));
445 Int_t ientry=(findexarray[ich]-1)*NIDX;
446 fbigarray[ich][ientry+DELTAIDXTOT]=tot; //in ns
447 fbigarray[ich][ientry+DELTAIDXTIME]=time*1E-3; // in ns
448 fbigarray[ich][ientry+DELTAIDXEXTIMEPI]=expTimePi;
449 fbigarray[ich][ientry+DELTAIDXEXTIMEKA]=expTimeKa;
450 fbigarray[ich][ientry+DELTAIDXEXTIMEPR]=expTimePr;
451 fhToT->Fill(fbigarray[ich][ientry+DELTAIDXTOT]);
452 fhTime->Fill(fbigarray[ich][ientry+DELTAIDXTIME]);
453 fhExpTimePi->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPI]);
454 fhExpTimeKa->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEKA]);
455 fhExpTimePr->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPR]);
456 AliDebug(2, Form("track = %i, tot = %f, time = %f, and Exp time in TOF: pi = %f, K = %f, p = %f",itr, fbigarray[ich][ientry+DELTAIDXTOT], fbigarray[ich][ientry+DELTAIDXTIME], expTimePi,expTimeKa,expTimePr));
457 }
458 fnESDselected+=nselected;
459
460 PostData(0, fOutputContainer);
461}
462
463//_____________________________________________________________________________
464void AliTOFCalibTask::Terminate(Option_t *)
465{
466 // Processing when the event loop is ended
467
468 // some plots
469
470 TH1::AddDirectory(0);
471 AliInfo("TOF Calib Task: End of events loop");
472 AliInfo(Form(" Number of analyzed ESD tracks: %i\n",fnESD));
473 AliInfo(Form(" Number of selected ESD tracks: %i\n",fnESDselected));
474 AliInfo(Form(" Number of ESD tracks with kTOFout: %i\n",fnESDkTOFout));
475 AliInfo(Form(" Number of ESD tracks with kTIME: %i\n",fnESDkTIME));
476 AliInfo(Form(" Number of ESD tracks with TRDcut: %i\n",fnESDTRDcut));
477 AliInfo(Form(" Number of ESD tracks with TIMEcut: %i\n",fnESDTIMEcut));
478 AliInfo(Form(" Number of ESD tracks with assTOFcl: %i\n",fnESDassTOFcl));
479
480 for (Int_t i = 0;i<TOFCHANNELS;i++){
481 Int_t size=findexarray[i]*NIDX;
482 AliDebug(2, Form(" entries %i in channel %i ",findexarray[i],i));
483 if (findexarray[i]<6) {
484 AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
485 continue;
486 }
487 if (!CombPID(&fbigarray[i][0], size)) AliError("ERROR!!!!ERROR!!!");
488 }
70ddd0ee 489
1ed77580 490 DrawHistos();
491
492 // saving data in a tree
493
494 TDirectory *dir = gDirectory;
495 TFile *outFile = 0;
496 TTree * tree = 0x0;
497 Bool_t isThere=kFALSE;
498 const char *dirname = "./";
499 TString filename= "TOFCalib.root";
500
501 if((gSystem->FindFile(dirname,filename))!=NULL){
502 isThere=kTRUE;
503 }
504
505 if(!isThere){
506 AliInfo("File with tree for Calibration not there, creating it");
507 outFile = new TFile( "TOFCalib.root","RECREATE");
508 outFile->cd();
509 tree = new TTree("T", "Tree for TOF calibration");
510 Float_t p[CHENTRIESSMALL];
511 Int_t nentries;
512 tree->Branch("nentries",&nentries,"nentries/I");
513 tree->Branch("TOFentries",p,"TOFentries[nentries]/F");
514 for (Int_t i=0;i<TOFCHANNELS;i++){
515 nentries=findexarray[i]*(NIDXSMALL); // when filling small array,
516 // only first 3 floats taken
517 // into account
518 for (Int_t j=0; j<findexarray[i];j++){
519 for (Int_t k=0; k<NIDXSMALL;k++){
520 Int_t index1= j*NIDXSMALL+k; // index in small array
521 Int_t index2=j*NIDX+k; // index in big array
522 p[index1]=fbigarray[i][index2];
523 }
524 }
525 tree->Fill();
526 }
527 tree->Write("",TObject::kOverwrite);
528 outFile->Close();
529 //delete outFile;
530 //outFile=0x0;
531 //delete tree;
532 //tree=0x0;
533 }
534 else{
535 AliInfo("File with tree for Calibration already there, updating it");
536 outFile = new TFile( "TOFCalib.root","READ");
537 outFile->cd();
538 tree=(TTree*)outFile->Get("T");
539 Float_t p[CHENTRIESSMALL];
540 Float_t ptemp[MAXCHENTRIESSMALL];
541 Int_t nentries;
542 Int_t nentries1;
543 Int_t nentriestemp;
544 tree->SetBranchAddress("nentries",&nentries);
545 tree->SetBranchAddress("TOFentries",p);
546 TFile * outFile1 = new TFile( "TOFCalib.root","RECREATE");
547 outFile1->cd();
548 TTree *treeTemp = new TTree("T","Tree for TOF calibration");
549 treeTemp->Branch("nentries",&nentriestemp,"nentries/I");
550 treeTemp->Branch("TOFentries",ptemp,"TOFentries[nentries]/F");
551 for (Int_t i=0;i<TOFCHANNELS;i++){
552 tree->GetEntry(i);
553 nentries1=nentries;
554 if (nentries+findexarray[i]*(NIDXSMALL) > MAXCHENTRIESSMALL){
555 AliDebug(2, Form("findexarray[%i] = %i, nentries = %i",i,findexarray[i],nentries));
556 AliInfo(Form("Too many entries in channel %i, stopping at %i entries, no storing any more...",i,MAXCHENTRIESSMALL/NIDXSMALL));
557 findexarray[i]=(MAXCHENTRIESSMALL-nentries)/NIDXSMALL;
558 }
559 for (Int_t kk = 0;kk<nentries;kk++){
560 ptemp[kk]=p[kk];
561 }
562 nentries+=findexarray[i]*(NIDXSMALL); // when filling small array,
563 // only first 3 floats taken
564 // into account
565 nentriestemp=nentries;
566 for (Int_t j=0; j<findexarray[i];j++){
567 for (Int_t k=0; k<NIDXSMALL;k++){
568 Int_t index1= j*NIDXSMALL+k; // index in small array
569 Int_t index2=j*NIDX+k; // index in big array
570 ptemp[index1+nentries1]=fbigarray[i][index2];
571 }
572 }
573 treeTemp->Fill();
574 }
575 treeTemp->SetName("T");
576 treeTemp->SetTitle("Tree for TOF calibration");
577 outFile->Close();
578 treeTemp->Write("",TObject::kOverwrite);
579 outFile1->Close();
580 }
581 dir->cd();
582 Int_t calibrationStatus = 0;
70ddd0ee 583 // Int_t ch[2]={3,1003};
584 // calibrationStatus = Calibrate(2,ch,"save");
1ed77580 585 //calibrationStatus = Calibrate(3,4,"save");
586 //calibrationStatus = CalibrateFromProfile(3,"save");
70ddd0ee 587 // calibrationStatus = Calibrate(3,"save");
588 calibrationStatus = Calibrate("");
1ed77580 589 AliInfo(Form("Calibration Status = %i, bye.... \n",calibrationStatus));
590
591}
592//----------------------------------------------------------------------------
593Int_t AliTOFCalibTask::Calibrate(Int_t ichmin, Int_t ichmax, Option_t *optionSave, Option_t *optionFit){
594
595 // calibrating summing more than one channels
596 // computing calibration parameters
597 // Returning codes:
598 // 0 -> everything was ok
599 // 1 -> no tree for calibration found
600 // 2 -> not enough statistics to perform calibration
601 // 3 -> problems with arrays
602
603 TH1::AddDirectory(0);
604
605 AliInfo(Form("*** Calibrating Histograms %s, summing more channels, from channel %i, to channel %i, storing Calib Pars in channel %i", GetName(),ichmin,ichmax,ichmin)) ;
606 AliInfo(Form("Option for Saving histos = %s",optionSave )) ;
607 AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ;
608 const char *dirname = "./";
609 TString filename= "TOFCalib.root";
610
611 if((gSystem->FindFile(dirname,filename))==NULL){
612 AliInfo("No file with tree for calibration found! exiting....");
613 return 1;
614 }
615
616 TFile *file = new TFile("TOFCalib.root","READ");
617 TTree *tree = (TTree*)file->Get("T");
618 Float_t p[MAXCHENTRIESSMALL];
619 Int_t nentries;
620 tree->SetBranchAddress("nentries",&nentries);
621 tree->SetBranchAddress("TOFentries",p);
622
623 Float_t nentriesmean =0;
624 for (Int_t i=ichmin; i<ichmax; i++){
625 tree->GetEntry(i);
626 Int_t ntracks=nentries/3;
627 nentriesmean+=ntracks;
628 }
629
630 // nentriesmean/=(ichmax-ichmin);
631 if (nentriesmean < MEANENTRIES) {
632 AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
633 return 2;
634 }
635
636 //filling ToT and Time arrays
637
638 Int_t nbinToT = 100; // ToT bin width in Profile = 48.8 ps
639 Float_t minToT = 0; // ns
640 Float_t maxToT = 4.88; // ns
641 TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
642
643 TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
644 TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
645 Int_t ntracksTot = 0;
646 Int_t ntracks = 0;
647 Double_t binsProfile[101]; // sized larger than necessary, the correct
648 // dim being set in the booking of the profile
649 Int_t nusefulbins=0;
650 Float_t meantime=0;
651 for (Int_t i = ichmin;i<ichmax;i++){
652 tree->GetEntry(i);
653 ntracksTot+=nentries/3;
654 ntracks=nentries/3;
655 AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i,nentries, ntracks));
656 for (Int_t j=0;j<ntracks;j++){
657 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
658 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
659 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
660 Float_t tot = p[idxexToT];
661 hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
662 meantime+=p[idxexTime]-p[idxexExTime];
663 hToT->Fill(tot);
664 }
665 }
666 nusefulbins = FindBins(hToT,&binsProfile[0]);
667 meantime/=ntracksTot;
668 AliDebug(2, Form("meantime = %f",meantime));
669
670 for (Int_t j=1;j<=nusefulbins;j++) {
671 AliDebug(2,Form(" summing channels from %i to %i, nusefulbins = %i, bin %i = %f",ichmin,ichmax,nusefulbins,j,binsProfile[j]));
672 }
673
674 TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G"); // CHECK THE BUILD OPTION, PLEASE!!!!!!
675 TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
676 for (Int_t i=ichmin; i<ichmax; i++){
677 tree->GetEntry(i);
678 ntracks=nentries/3;
679 for (Int_t j=0;j<ntracks;j++){
680 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
681 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
682 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
683 Float_t tot = p[idxexToT];
684 Float_t time = p[idxexTime]-p[idxexExTime];
685 AliDebug (2, Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
686 hSlewingProf->Fill(tot,time); // if meantime is not used, the fill may be moved in the loop above
687 htimetot->Fill(tot,time-meantime); // if meantime is not used, the fill may be moved in the loop above
688 }
689 }
690 hSlewingProf->Fit("pol5",optionFit,"",0,4);
691 TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
692 Float_t par[6];
693 for(Int_t kk=0;kk<6;kk++){
694 par[kk]=calibfunc->GetParameter(kk);
695 AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
696 }
697
698 if(strstr(optionSave,"save")){
699 TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
700 fileProf->cd();
701 TString profName=Form("Profile%06i_%06i",ichmin,ichmax);
702 TString timeTotName=Form("TimeTot%06i_%06i",ichmin,ichmax);
703 TString totName=Form("Tot%06i_%06i",ichmin,ichmax);
704 TString deltaName=Form("Delta%06i_%06i",ichmin,ichmax);
705 hSlewingProf->Write(profName);
706 htimetot->Write(timeTotName);
707 hToT->Write(totName);
708 hdeltaTime->Write(deltaName);
709 fileProf->Close();
710 delete fileProf;
711 fileProf=0x0;
712 }
713
714 delete hToT;
715 hToT=0x0;
716 delete hSlewingProf;
717 hSlewingProf=0x0;
718 delete htimetot;
719 htimetot=0x0;
720 delete hdeltaTime;
721 hdeltaTime=0x0;
722 file->Close();
723 delete file;
724 file=0x0;
725
726 AliTOFChannelTask * calChannel = new AliTOFChannelTask();
727 calChannel->SetSlewPar(par);
728 // saving parameters in chmin
729 arrayCal->AddAt(calChannel,ichmin);
730 TFile *filecalib = new TFile("outCalArray.root","RECREATE");
731 arrayCal->Write("array",TObject::kSingleKey);
732 filecalib->Close();
733 delete filecalib;
734 filecalib = 0x0;
735 delete arrayCal;
736 arrayCal = 0x0;
737 delete calChannel;
738 calChannel =0x0;
739 return 0;
740}
741//----------------------------------------------------------------------------
742Int_t AliTOFCalibTask::Calibrate(Int_t i, Option_t *optionSave, Option_t *optionFit){
743
744 // computing calibration parameters for channel i
745 // Returning codes:
746 // 0 -> everything was ok
747 // 1 -> no tree for calibration found
748 // 2 -> not enough statistics to perform calibration
749 // 3 -> problems with arrays
750
751 TH1::AddDirectory(0);
752
753 AliInfo(Form("*** Calibrating Histograms (one channel) %s", GetName())) ;
754 AliInfo(Form("Option for Saving histos = %s",optionSave )) ;
755 AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ;
756 const char *dirname = "./";
757 TString filename= "TOFCalib.root";
758
759 if((gSystem->FindFile(dirname,filename))==NULL){
760 AliInfo("No file with tree for calibration found! exiting....");
761 return 1;
762 }
763
764 TFile *file = new TFile("TOFCalib.root","READ");
765 TTree *tree = (TTree*)file->Get("T");
766 Float_t p[MAXCHENTRIESSMALL];
767 Int_t nentries;
768 tree->SetBranchAddress("nentries",&nentries);
769 tree->SetBranchAddress("TOFentries",p);
770
771 tree->GetEntry(i);
772 Int_t nentriesmean=nentries/3;
773
774 if (nentriesmean < MEANENTRIES) {
775 AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
776 return 2;
777 }
778
779 //filling ToT and Time arrays
780
781 Int_t nbinToT = 100; // ToT bin width in Profile = 48.8 ps
782 Float_t minToT = 0; // ns
783 Float_t maxToT = 4.88; // ns
784 TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
785
786 TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
787 TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
788 Int_t ntracks = 0;
789 Double_t binsProfile[101]; // sized larger than necessary, the correct
790 // dim being set in the booking of the profile
791 Int_t nusefulbins=0;
792 Float_t meantime=0;
793 ntracks=nentries/3;
794 AliDebug(2,Form("channel %i, nentries = %i, ntracks = %i",i ,nentries, ntracks));
795 for (Int_t j=0;j<ntracks;j++){
796 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
797 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
798 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
799 Float_t tot = p[idxexToT];
800 meantime+=p[idxexTime]-p[idxexExTime];
801 hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
802 hToT->Fill(tot);
803 }
804
805 nusefulbins = FindBins(hToT,&binsProfile[0]);
806 meantime/=ntracks;
807 AliDebug(2,Form("meantime = %f",meantime));
808
809 for (Int_t j=1;j<=nusefulbins;j++) {
810 AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j]));
811 }
812
813 TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G"); // CHECK THE BUILD OPTION, PLEASE!!!!!!
814 TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
815 tree->GetEntry(i);
816 ntracks=nentries/3;
817 for (Int_t j=0;j<ntracks;j++){
818 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
819 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
820 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
821 Float_t tot = p[idxexToT];
822 Float_t time = p[idxexTime]-p[idxexExTime];
823 AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
824 hSlewingProf->Fill(tot,time); // if meantime is not used, the fill may be moved in the loop above
825 htimetot->Fill(tot,time-meantime); // if meantime is not used, the fill may be moved in the loop above
826 }
827 hSlewingProf->Fit("pol5",optionFit,"",0,4);
828 TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
829 Float_t par[6];
830 for(Int_t kk=0;kk<6;kk++){
831 par[kk]=calibfunc->GetParameter(kk);
832 AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
833 }
834
835 file->Close();
836 delete file;
837 file=0x0;
838
839 if(strstr(optionSave,"save")){
840 TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
841 fileProf->cd();
842 TString profName=Form("Profile%06i",i);
843 TString timeTotName=Form("TimeTot%06i",i);
844 TString totName=Form("Tot%06i",i);
845 TString deltaName=Form("Delta%06i",i);
846 hSlewingProf->Write(profName);
847 htimetot->Write(timeTotName);
848 hToT->Write(totName);
849 hdeltaTime->Write(deltaName);
850 fileProf->Close();
851 delete fileProf;
852 fileProf=0x0;
853 }
854
855 delete hToT;
856 hToT=0x0;
857 delete hSlewingProf;
858 hSlewingProf=0x0;
859 delete htimetot;
860 htimetot=0x0;
861 delete hdeltaTime;
862 hdeltaTime=0x0;
863
864 AliTOFChannelTask * calChannel = new AliTOFChannelTask();
865 calChannel->SetSlewPar(par);
866 arrayCal->AddAt(calChannel,i);
867 TFile *filecalib = new TFile("outCalArray.root","RECREATE");
868 arrayCal->Write("array",TObject::kSingleKey);
869 filecalib->Close();
870 delete filecalib;
871 filecalib=0x0;
872 delete arrayCal;
873 arrayCal = 0x0;
874 delete calChannel;
875 calChannel =0x0;
876
877 return 0;
878}
879//----------------------------------------------------------------------------
880Int_t AliTOFCalibTask::Calibrate(Int_t nch, Int_t *ch, Option_t *optionSave, Option_t *optionFit){
881
882 // calibrating an array of channels
883 // computing calibration parameters
884 // Returning codes:
885 // 0 -> everything was ok
886 // 1 -> no tree for calibration found
887 // 2 -> not enough statistics to perform calibration
888 // 3 -> problems with arrays
889
890 TH1::AddDirectory(0);
891
892 AliInfo(Form("*** Calibrating Histograms %s, number of channels = %i", GetName(),nch)) ;
893 AliInfo(Form("Option for Saving histos = %s",optionSave )) ;
894 AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ;
895 for (Int_t ich=0; ich<nch; ich++){
896 Int_t i = ch[ich];
897 AliInfo(Form("Calibrating channel = %i",i )) ;
898 }
899 const char *dirname = "./";
900 TString filename= "TOFCalib.root";
901
902 if((gSystem->FindFile(dirname,filename))==NULL){
903 AliInfo("No file with tree for calibration found! exiting....");
904 return 1;
905 }
906
907 TFile *file = new TFile("TOFCalib.root","READ");
908 TTree *tree = (TTree*)file->Get("T");
909 Float_t p[MAXCHENTRIESSMALL];
910 Int_t nentries;
911 tree->SetBranchAddress("nentries",&nentries);
912 tree->SetBranchAddress("TOFentries",p);
913
914 Float_t nentriesmean =0;
915 for (Int_t ich=0; ich<nch; ich++){
916 Int_t i = ch[ich];
917 tree->GetEntry(i);
918 Int_t ntracks=nentries/3;
919 nentriesmean+=ntracks;
920 }
921
922 nentriesmean/=nch;
923 if (nentriesmean < MEANENTRIES) {
924 AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
925 return 2;
926 }
927
928 //filling ToT and Time arrays
929
930 Int_t nbinToT = 100; // ToT bin width in Profile = 48.8 ps
931 Float_t minToT = 0; // ns
932 Float_t maxToT = 4.88; // ns
933 TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
934 arrayCal->SetOwner();
935 TFile * fileProf=0x0;
936 if(strstr(optionSave,"save")){
937 fileProf = new TFile("TOFCalibSave.root","recreate");
938 }
939 for (Int_t ich=0; ich<nch; ich++) {
940 TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
941 TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
942 Double_t binsProfile[101]; // sized larger than necessary, the correct
943 // dim being set in the booking of the profile
944 TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
945 Int_t ntracksTot = 0;
946 Int_t ntracks = 0;
947 Int_t nusefulbins=0;
948 Float_t meantime=0;
949 Int_t i = ch[ich];
950 AliDebug(2,Form("Calibrating channel %i",i));
951 tree->GetEntry(i);
952 ntracksTot+=nentries/3;
953 ntracks=nentries/3;
954 if (ntracks < MEANENTRIES) {
955 AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
956 continue;
957 }
958 for (Int_t j=0;j<ntracks;j++){
959 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
960 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
961 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
962 Float_t tot = p[idxexToT];
963 hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
964 meantime+=p[idxexTime]-p[idxexExTime];
965 hToT->Fill(tot);
966 }
967
968 nusefulbins = FindBins(hToT,&binsProfile[0]);
969 meantime/=ntracksTot;
970 for (Int_t j=1;j<=nusefulbins;j++) {
971 AliDebug(2,Form(" channel %i, nusefulbins = %i, bin %i = %f",i,nusefulbins,j,binsProfile[j]));
972 }
973
974 TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G"); // CHECK THE BUILD OPTION, PLEASE!!!!!!
975 tree->GetEntry(i);
976 ntracks=nentries/3;
977 for (Int_t j=0;j<ntracks;j++){
978 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
979 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
980 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
981 Float_t tot = p[idxexToT];
982 Float_t time = p[idxexTime]-p[idxexExTime];
983 AliDebug(2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
984 hSlewingProf->Fill(tot,time); // if meantime is not used, the fill may be moved in the loop above
985 htimetot->Fill(tot,time-meantime); // if meantime is not used, the fill may be moved in the loop above
986 }
987
988 hSlewingProf->Fit("pol5",optionFit,"",1,4);
989 TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
990 Float_t par[6];
991 for(Int_t kk=0;kk<6;kk++){
992 par[kk]=calibfunc->GetParameter(kk);
993 AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
994 }
995
996 if(strstr(optionSave,"save") && fileProf){
997 TString profName=Form("Profile%06i",i);
998 TString timeTotName=Form("TimeTot%06i",i);
999 TString totName=Form("Tot%06i",i);
1000 TString deltaName=Form("Delta%06i",i);
1001 fileProf->cd();
1002 hSlewingProf->Write(profName);
1003 htimetot->Write(timeTotName);
1004 hToT->Write(totName);
1005 hdeltaTime->Write(deltaName);
1006 }
1007
1008 AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1009 calChannel->SetSlewPar(par);
1010 arrayCal->AddAt(calChannel,i);
1011 delete hToT;
1012 hToT=0x0;
1013 delete hSlewingProf;
1014 hSlewingProf=0x0;
1015 delete htimetot;
1016 htimetot=0x0;
1017 delete hdeltaTime;
1018 hdeltaTime=0x0;
1019 }
1020
1021 if(strstr(optionSave,"save") && fileProf){
1022 fileProf->Close();
1023 delete fileProf;
1024 fileProf=0x0;
1025 }
1026 file->Close();
1027 delete file;
1028 file=0x0;
1029 TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1030 arrayCal->Write("array",TObject::kSingleKey);
1031 filecalib->Close();
1032 delete filecalib;
1033 filecalib=0x0;
1034 arrayCal->Clear();
1035 delete arrayCal;
1036 arrayCal=0x0;
1037
1038 return 0;
1039}
1040//----------------------------------------------------------------------------
1041Int_t AliTOFCalibTask::CalibrateFromProfile(Int_t i, Option_t *optionSave, Option_t *optionFit){
1042
1043 // computing calibration parameters using the old profiling algo
1044 // Returning codes:
1045 // 0 -> everything was ok
1046 // 1 -> no tree for calibration found
1047 // 2 -> not enough statistics to perform calibration
1048 // 3 -> problems with arrays
1049
1050 TH1::AddDirectory(0);
1051
1052 TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1053 AliInfo(Form("*** Calibrating Histograms From Profile %s", GetName())) ;
1054 AliInfo(Form("Option for Saving histos = %s",optionSave )) ;
1055 AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ;
1056 const char *dirname = "./";
1057 TString filename= "TOFCalib.root";
1058
1059 if((gSystem->FindFile(dirname,filename))==NULL){
1060 AliInfo("No file with tree for calibration found! exiting....");
1061 return 1;
1062 }
1063
1064 TFile *file = new TFile("TOFCalib.root","READ");
1065 TTree *tree = (TTree*)file->Get("T");
1066 Float_t p[MAXCHENTRIESSMALL];
1067 Int_t nentries;
1068 tree->SetBranchAddress("nentries",&nentries);
1069 tree->SetBranchAddress("TOFentries",p);
1070
1071 tree->GetEntry(i);
1072 Int_t ntracks=nentries/3;
1073
1074 if (ntracks < MEANENTRIES) {
1075 AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",ntracks));
1076 return 2;
1077 }
1078
1079 TH1F * hProf = new TH1F();
1080 hProf = Profile(i);
1081 hProf->Fit("pol5",optionFit,"",0,4);
1082 TF1 * calibfunc = (TF1*)hProf->GetFunction("pol5");
1083 Float_t par[6];
1084 for(Int_t kk=0;kk<6;kk++){
1085 par[kk]=calibfunc->GetParameter(kk);
1086 AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1087 }
1088
1089 if(strstr(optionSave,"save")){
1090 TFile * fileProf = new TFile("TOFCalibSave.root","recreate");
1091 fileProf->cd();
1092 TString profName=Form("Profile%06i",i);
1093 hProf->Write(profName);
1094 fileProf->Close();
1095 delete fileProf;
1096 fileProf=0x0;
1097 }
1098
1099 delete hProf;
1100 hProf=0x0;
1101 file->Close();
1102 delete file;
1103 file=0x0;
1104 AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1105 calChannel->SetSlewPar(par);
1106 arrayCal->AddAt(calChannel,i);
1107 TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1108 arrayCal->Write("array",TObject::kSingleKey);
1109 filecalib->Close();
1110 delete filecalib;
1111 filecalib = 0x0;
1112 delete calChannel;
1113 delete arrayCal;
1114 return 0;
1115}
1116//----------------------------------------------------------------------------
1117Int_t AliTOFCalibTask::Calibrate(Option_t *optionSave, Option_t *optionFit){
1118
1119 // calibrating the whole TOF
1120 // computing calibration parameters
1121 // Returning codes:
1122 // 0 -> everything was ok
1123 // 1 -> no tree for calibration found
1124 // 2 -> not enough statistics to perform calibration
1125 // 3 -> problems with arrays
1126
1127 TH1::AddDirectory(0);
1128
1129 AliInfo(Form("*** Calibrating Histograms %s, all channels", GetName())) ;
1130 AliInfo(Form("Option for Saving histos = %s",optionSave )) ;
1131 AliInfo(Form("Option for Fitting Profile histos = %s",optionFit )) ;
1132 const char *dirname = "./";
1133 TString filename= "TOFCalib.root";
1134
1135 if((gSystem->FindFile(dirname,filename))==NULL){
1136 AliInfo("No file with tree for calibration found! exiting....");
1137 return 1;
1138 }
1139
1140 TFile * fileProf=0x0;
1141 if(strstr(optionSave,"save")){
1142 fileProf = new TFile("TOFCalibSave.root","recreate");
1143 }
1144 TFile *file = new TFile("TOFCalib.root","READ");
1145 TTree *tree = (TTree*)file->Get("T");
1146 Float_t p[MAXCHENTRIESSMALL];
1147 Int_t nentries;
1148 tree->SetBranchAddress("nentries",&nentries);
1149 tree->SetBranchAddress("TOFentries",p);
1150
1151 Float_t nentriesmean =0;
1152 for (Int_t i=0; i<TOFCHANNELS; i++){
1153 tree->GetEntry(i);
1154 Int_t ntracks=nentries/3;
1155 nentriesmean+=ntracks;
1156 }
1157
1158 nentriesmean/=TOFCHANNELS;
1159 if (nentriesmean < MEANENTRIES) {
1160 AliInfo(Form(" Too small mean number of entires per channel (mean number = %f) not calibrating and exiting.....",nentriesmean));
1161 return 2;
1162 }
1163
1164 //filling ToT and Time arrays
1165
1166 Int_t nbinToT = 100; // ToT bin width in Profile = 50.0 ps
1167 Float_t minToT = 0; // ns
1168 Float_t maxToT = 4.88;// ns
1169 TObjArray * arrayCal = new TObjArray(TOFCHANNELS);
1170 arrayCal->SetOwner();
1171 for (Int_t ii=0; ii<TOFCHANNELS; ii++) {
1172 TH1F *hToT = new TH1F("htot","htot",nbinToT, minToT, maxToT);
1173 TH1F *hdeltaTime = new TH1F("hdeltaTime","hdeltaTime",200,2,4);
1174 TH2F * htimetot = new TH2F("htimetot","htimetot",nbinToT, minToT, maxToT,600,-5,10);
1175 if (ii%1000 == 0) AliDebug(1,Form("Calibrating channel %i ",ii));
1176 Int_t i = 3;
1177 Int_t nusefulbins=0;
1178 Double_t binsProfile[101]; // sized larger than necessary, the correct
1179 // dim being set in the booking of the profile
1180 tree->GetEntry(i);
1181 Int_t ntracks=nentries/3;
1182 if (ntracks < MEANENTRIES) {
1183 AliInfo(Form(" Too small mean number of entires in channel %i (number of tracks = %f), not calibrating channel and continuing.....",i,ntracks));
1184 continue;
1185 }
1186 Float_t meantime=0;
1187 for (Int_t j=0;j<ntracks;j++){
1188 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
1189 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
1190 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
1191 Float_t tot = p[idxexToT];
1192 hdeltaTime->Fill(p[idxexTime]-p[idxexExTime]);
1193 meantime+=p[idxexTime]-p[idxexExTime];
1194 hToT->Fill(tot);
1195 }
1196 nusefulbins = FindBins(hToT,&binsProfile[0]);
1197 meantime/=ntracks;
1198 for (Int_t j=0;j<nusefulbins;j++) {
1199 AliDebug(2,Form(" channel %i, usefulbin = %i, low edge = %f",i,j,binsProfile[j]));
1200 }
1201 TProfile* hSlewingProf = new TProfile("hSlewingProf", "hSlewingProf",nusefulbins, binsProfile, "G"); // CHECK THE BUILD OPTION, PLEASE!!!!!!
1202 for (Int_t j=0;j<ntracks;j++){
1203 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
1204 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
1205 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
1206 Float_t tot = p[idxexToT];
1207 Float_t time = p[idxexTime]-p[idxexExTime];
1208 AliDebug (2,Form("track = %i, time = %f, tot = %f, time-meantime = %f",j,time, tot, time-meantime));
70ddd0ee 1209 hSlewingProf->Fill(tot,time); // if meantime is not used, the fill may be moved in the loop above
1ed77580 1210 htimetot->Fill(tot,time-meantime); // if meantime is not used, the fill may be moved in the loop above
1211 }
1212
1213 // hSlewingProf->Fit("pol5",optionFit,"",1,4);
1214 //TF1 * calibfunc = (TF1*)hSlewingProf->GetFunction("pol5");
1215 Float_t par[6];
1216 for(Int_t kk=0;kk<6;kk++){
1217 // par[kk]=calibfunc->GetParameter(kk);
1218 par[kk]=kk;
1219 AliDebug(2,Form("parameter %i = %f",kk,par[kk]));
1220 }
1221
1222 if(strstr(optionSave,"save") && fileProf){
1223 TString profName=Form("Profile%06i",ii);
1224 TString timeTotName=Form("TimeTot%06i",ii);
1225 TString totName=Form("Tot%06i",ii);
1226 TString deltaName=Form("Delta%06i",ii);
1227 fileProf->cd();
1228 hSlewingProf->Write(profName);
1229 htimetot->Write(timeTotName);
1230 hToT->Write(totName);
1231 hdeltaTime->Write(deltaName);
1232 }
1233 AliTOFChannelTask * calChannel = new AliTOFChannelTask();
1234 calChannel->SetSlewPar(par);
1235 arrayCal->AddAt(calChannel,ii);
1236
1237 delete hToT;
1238 hToT=0x0;
1239 delete hSlewingProf;
1240 hSlewingProf=0x0;
1241 delete htimetot;
1242 htimetot=0x0;
1243 delete hdeltaTime;
1244 hdeltaTime=0x0;
1245 }
1246
1247 if(strstr(optionSave,"save")){
1248 fileProf->Close();
1249 delete fileProf;
1250 fileProf=0x0;
1251 }
1252
1253 file->Close();
1254 delete file;
1255 file=0x0;
1256 TFile *filecalib = new TFile("outCalArray.root","RECREATE");
1257 arrayCal->Write("array", TObject::kSingleKey);
1258 filecalib->Close();
1259 delete filecalib;
1260 filecalib=0x0;
1261 arrayCal->Clear();
1262 delete arrayCal;
1263 arrayCal=0x0;
1264
1265 return 0;
1266}
1267
1268//_____________________________________________________________________________
1269
1270Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
1271
1272 // selecting good quality tracks
1273 //track selection: reconstrution to TOF:
1274 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
1275 return 0;
1276 }
1277 fnESDkTOFout++;
1278 //IsStartedTimeIntegral
1279 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
1280 return 0;
1281 }
1282 fnESDkTIME++;
1283 if (t->GetStatus() & AliESDtrack::kTRDbackup) {
1284 Float_t xout= t->GetOuterParam()->GetX();
1285 if (xout<364.25 && xout > 300.) return 0;
1286 }
1287 fnESDTRDcut++;
1288 Double_t time=t->GetTOFsignal();
1289 time*=1.E-3; // tof given in nanoseconds
1290 if(time >= fMinTime){
1291 return 0;
1292 }
1293 fnESDTIMEcut++;
1294
1295 Double_t mom=t->GetP();
1296 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
1297 // return 0; // skipping momentum cut
1298 }
1299
1300 UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
1301 if(assignedTOFcluster==0){ // not matched
1302 return 0;
1303 }
1304 fnESDassTOFcl++;
1305 return 1;
1306}
1307//_____________________________________________________________________________
1308
1309Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
1310
1311 //track Combinatorial PID for calibration
1312
1313 //cout << " In comb PID " << endl;
1314 Float_t t0offset=0;
1315 const Int_t kntracksinset=6;
1316 if (kntracksinset != 6) {
1317 AliDebug(1,"Number of tracks in set smaller than expected one, identifying every particle as if it was a pion!");
1318 return 0;
1319 }
1320 Float_t exptof[kntracksinset][3];
1321 Int_t assparticle[kntracksinset];
1322 Float_t timeofflight[kntracksinset];
1323 Float_t texp[kntracksinset];
1324 Float_t timezero[kntracksinset];
1325 Float_t weightedtimezero[kntracksinset];
1326 Float_t besttimezero[kntracksinset];
1327 Float_t bestchisquare[kntracksinset];
1328 Float_t bestweightedtimezero[kntracksinset];
1329
1330 for (Int_t i=0;i<kntracksinset;i++){
1331 assparticle[i]=3;
1332 timeofflight[i]=0;
1333 texp[i]=0;
1334 timezero[i]=0;
1335 weightedtimezero[i]=0;
1336 besttimezero[i]=0;
1337 bestchisquare[i]=0;
1338 bestweightedtimezero[i]=0;
1339 }
1340
1341 Int_t nset= (Int_t)(size/kntracksinset/NIDX);
1342 for (Int_t i=0; i< nset; i++) {
1343 for (Int_t j=0; j<kntracksinset; j++) {
1344 Int_t idxtime = ((kntracksinset*i+j)*NIDX)+DELTAIDXTIME;
1345 Int_t idxextimePi = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPI;
1346 Int_t idxextimeKa = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEKA;
1347 Int_t idxextimePr = ((kntracksinset*i+j)*NIDX)+DELTAIDXEXTIMEPR;
1348 Double_t time=smallarray[idxtime]; // TOF time in ns
1349 timeofflight[j]=time+t0offset;
1350 exptof[j][0]=smallarray[idxextimePi];
1351 exptof[j][1]=smallarray[idxextimeKa];
1352 exptof[j][2]=smallarray[idxextimePr];
1353 AliDebug(2,Form("j = %i, Time = %f, and Exp time in PID: pi = %f, K = %f, p = %f",j,timeofflight[j],exptof[j][0],exptof[j][1],exptof[j][2]));
1354 }
1355 Float_t t0best=999.;
1356 Float_t et0best=999.;
1357 Float_t chisquarebest=999.;
1358 for (Int_t i1=0; i1<3;i1++) {
1359 texp[0]=exptof[0][i1];
1360 for (Int_t i2=0; i2<3;i2++) {
1361 texp[1]=exptof[1][i2];
1362 for (Int_t i3=0; i3<3;i3++) {
1363 texp[2]=exptof[2][i3];
1364 for (Int_t i4=0; i4<3;i4++) {
1365 texp[3]=exptof[3][i4];
1366 for (Int_t i5=0; i5<3;i5++) {
1367 texp[4]=exptof[4][i5];
1368 for (Int_t i6=0; i6<3;i6++) {
1369 texp[5]=exptof[5][i6];
1370
1371 Float_t sumAllweights=0.;
1372 Float_t meantzero=0.;
1373 Float_t emeantzero=0.;
1374
1375 for (Int_t itz=0; itz<kntracksinset;itz++) {
1376 timezero[itz]=texp[itz]-timeofflight[itz];
1377 weightedtimezero[itz]=timezero[itz]/TRACKERROR;
1378 sumAllweights+=1./TRACKERROR;
1379 meantzero+=weightedtimezero[itz];
1380
1381 } // end loop for (Int_t itz=0; itz<15;itz++)
1382
1383 meantzero=meantzero/sumAllweights; // it is given in [ns]
1384 emeantzero=sqrt(1./sumAllweights); // it is given in [ns]
1385
1386 // calculate chisquare
1387
1388 Float_t chisquare=0.;
1389 for (Int_t icsq=0; icsq<kntracksinset;icsq++) {
1390 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
1391 } // end loop for (Int_t icsq=0; icsq<15;icsq++)
1392
1393 Int_t npion=0;
1394 if(i1==0)npion++;
1395 if(i2==0)npion++;
1396 if(i3==0)npion++;
1397 if(i4==0)npion++;
1398 if(i5==0)npion++;
1399 if(i6==0)npion++;
1400
1401 if(chisquare<=chisquarebest && ((Float_t) npion/ ((Float_t) kntracksinset)>0.3)){
1402 for(Int_t iqsq = 0; iqsq<kntracksinset; iqsq++) {
1403 besttimezero[iqsq]=timezero[iqsq];
1404 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
1405 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR;
1406 }
1407
1408 assparticle[0]=i1;
1409 assparticle[1]=i2;
1410 assparticle[2]=i3;
1411 assparticle[3]=i4;
1412 assparticle[4]=i5;
1413 assparticle[5]=i6;
1414
1415 chisquarebest=chisquare;
1416 t0best=meantzero;
1417 et0best=emeantzero;
1418 } // close if(dummychisquare<=chisquare)
1419 } // end loop on i6
1420 } // end loop on i5
1421 } // end loop on i4
1422 } // end loop on i3
1423 } // end loop on i2
1424 } // end loop on i1
1425
1426 Float_t confLevel=999;
1427 if(chisquarebest<999.){
1428 Double_t dblechisquare=(Double_t)chisquarebest;
1429 confLevel=(Float_t)TMath::Prob(dblechisquare,kntracksinset-1);
1430 }
1431 // assume they are all pions for fake sets
1432 if(confLevel<0.01 || confLevel==999. ){
1433 for (Int_t itrk=0; itrk<kntracksinset; itrk++)assparticle[itrk]=0;
1434 }
1435
1436 AliDebug(2,Form(" Best Assignment, selection %i %i %i %i %i %i ",assparticle[0],assparticle[1],assparticle[2],assparticle[3],assparticle[4],assparticle[5]));
1437
1438 for (Int_t kk=0;kk<kntracksinset;kk++){
1439 Int_t idxextimePi = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPI;
1440 Int_t idxextimeKa = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEKA;
1441 Int_t idxextimePr = ((kntracksinset*i+kk)*NIDX)+DELTAIDXEXTIMEPR;
1442 // storing in third slot of smallarray the assigned expected time
1443 fhPID->Fill(assparticle[kk]);
70ddd0ee 1444 // assparticle[kk]=0; //assuming all particles are pions
1ed77580 1445 if (assparticle[kk]==0){ //assigned to be a Pi
1446 smallarray[idxextimeKa]=-1;
1447 smallarray[idxextimePr]=-1;
1448 }
1449 else if (assparticle[kk]==1){ //assigned to be a Ka
1450 smallarray[idxextimePi]=smallarray[idxextimeKa];
1451 smallarray[idxextimeKa]=-1;
1452 smallarray[idxextimePr]=-1;
1453 }
1454 else if (assparticle[kk]==2){ //assigned to be a Pr
1455 smallarray[idxextimePi]=smallarray[idxextimePr];
1456 smallarray[idxextimeKa]=-1;
1457 smallarray[idxextimePr]=-1;
1458 }
1459 }
1460 }
1461 return 1;
1462}
1463//-----------------------------------------------------------------------
1464TH1F* AliTOFCalibTask::Profile(Int_t ich)
1465{
1466 // profiling algo
1467
1468 TFile *file = new TFile("TOFCalib.root","READ");
1469 TTree *tree = (TTree*)file->Get("T");
1470 Float_t p[MAXCHENTRIESSMALL];
1471 Int_t nentries;
1472 tree->SetBranchAddress("nentries",&nentries);
1473 tree->SetBranchAddress("TOFentries",p);
1474
1475 //Prepare histograms for Slewing Correction
1476 const Int_t knbinToT = 100;
1477 Int_t nbinTime = 200;
1478 Float_t minTime = -5.5; //ns
1479 Float_t maxTime = 5.5; //ns
1480 Float_t minToT = 0; //ns
1481 Float_t maxToT = 5.; //ns
1482 Float_t deltaToT = (maxToT-minToT)/knbinToT;
1483 Double_t mTime[knbinToT+1],mToT[knbinToT+1],meanTime[knbinToT+1], meanTime2[knbinToT+1],vToT[knbinToT+1], vToT2[knbinToT+1],meanToT[knbinToT+1],meanToT2[knbinToT+1],vTime[knbinToT+1],vTime2[knbinToT+1],xlow[knbinToT+1],sigmaTime[knbinToT+1];
1484 Int_t n[knbinToT+1], nentrx[knbinToT+1];
1485 Double_t sigmaToT[knbinToT+1];
1486 for (Int_t i = 0; i < knbinToT+1 ; i++){
1487 mTime[i]=0;
1488 mToT[i]=0;
1489 n[i]=0;
1490 meanTime[i]=0;
1491 meanTime2[i]=0;
1492 vToT[i]=0;
1493 vToT2[i]=0;
1494 meanToT[i]=0;
1495 meanToT2[i]=0;
1496 vTime[i]=0;
1497 vTime2[i]=0;
1498 xlow[i]=0;
1499 sigmaTime[i]=0;
1500 sigmaToT[i]=0;
1501 n[i]=0;
1502 nentrx[i]=0;
1503 }
1504 TH2F* hSlewing = new TH2F("hSlewing", "hSlewing", knbinToT, minToT, maxToT, nbinTime, minTime, maxTime);
1505 Int_t ntracks = 0;
1506 TH1F *histo = new TH1F("histo", "1D Time vs ToT", knbinToT, minToT, maxToT);
1507 tree->GetEntry(ich);
1508 ntracks=nentries/3;
1509 for (Int_t j=0;j<ntracks;j++){
1510 Int_t idxexToT = (j* NIDXSMALL)+DELTAIDXTOT;
1511 Int_t idxexTime = (j* NIDXSMALL)+DELTAIDXTIME;
1512 Int_t idxexExTime = (j* NIDXSMALL)+DELTAIDXPID;
1513 Float_t tot = p[idxexToT];
1514 Float_t time = p[idxexTime]-p[idxexExTime];
1515 Int_t nx = (Int_t)((tot-minToT)/deltaToT)+1;
1516 if ((tot != 0) && ( time!= 0)){
1517 vTime[nx]+=time;
1518 vTime2[nx]+=time*time;
1519 vToT[nx]+=tot;
1520 vToT2[nx]+=tot*tot;
1521 nentrx[nx]++;
1522 hSlewing->Fill(tot,time);
1523 }
1524 }
1525 Int_t nbinsToT=hSlewing->GetNbinsX();
1526 if (nbinsToT != knbinToT) {
1527 AliError("Profile :: incompatible numbers of bins");
1528 return 0x0;
1529 }
1530
1531 Int_t usefulBins=0;
1532 for (Int_t i=1;i<=nbinsToT;i++){
1533 if (nentrx[i]!=0){
1534 n[usefulBins]+=nentrx[i];
1535 if (n[usefulBins]==0 && i == nbinsToT) {
1536 break;
1537 }
1538 meanTime[usefulBins]+=vTime[i];
1539 meanTime2[usefulBins]+=vTime2[i];
1540 meanToT[usefulBins]+=vToT[i];
1541 meanToT2[usefulBins]+=vToT2[i];
1542 if (n[usefulBins]<10 && i!=nbinsToT) continue;
1543 mTime[usefulBins]=meanTime[usefulBins]/n[usefulBins];
1544 mToT[usefulBins]=meanToT[usefulBins]/n[usefulBins];
1545 sigmaTime[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1546 *(meanTime2[usefulBins]-meanTime[usefulBins]
1547 *meanTime[usefulBins]/n[usefulBins]));
1548 if ((1./n[usefulBins]/n[usefulBins]
1549 *(meanToT2[usefulBins]-meanToT[usefulBins]
1550 *meanToT[usefulBins]/n[usefulBins]))< 0) {
1551 AliError(" too small radical" );
1552 sigmaToT[usefulBins]=0;
1553 }
1554 else{
1555 sigmaToT[usefulBins]=TMath::Sqrt(1./n[usefulBins]/n[usefulBins]
1556 *(meanToT2[usefulBins]-meanToT[usefulBins]
1557 *meanToT[usefulBins]/n[usefulBins]));
1558 }
1559 usefulBins++;
1560 }
1561 }
1562 for (Int_t i=0;i<usefulBins;i++){
1563 Int_t binN = (Int_t)((mToT[i]-minToT)/deltaToT)+1;
1564 histo->Fill(mToT[i],mTime[i]);
1565 histo->SetBinError(binN,sigmaTime[i]);
1566 }
1567 delete file;
1568 file=0x0;
1569 delete hSlewing;
1570 hSlewing=0x0;
1571
1572 return histo;
1573}
1574//----------------------------------------------------------------------------
1575Int_t AliTOFCalibTask::FindBins(TH1F* h, Double_t *binsProfile) const{
1576
1577 // to determine the bins for ToT histo
1578
1579 Int_t cont = 0;
1580 Int_t startBin = 1;
1581 Int_t nbin = h->GetNbinsX();
70ddd0ee 1582 Int_t nentries = (Int_t)h->GetEntries();
1ed77580 1583 Float_t max = h->GetBinLowEdge(nbin);
1584 Int_t nusefulbins=0;
1585 Int_t maxcont=0;
1586 // setting maxvalue of entries per bin
70ddd0ee 1587 if (nentries <= 60) maxcont = 2;
1588 else if (nentries <= 100) maxcont = 5;
1589 else if (nentries <= 500) maxcont = 10;
1ed77580 1590 else maxcont = 20;
1591 for (Int_t j=1;j<=nbin;j++) {
1592 cont += (Int_t)h->GetBinContent(j);
1593 if (j<nbin){
1594 if (cont>=maxcont){
1595 nusefulbins++;
1596 binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1597 cont=0;
1598 startBin=j+1;
1599 continue;
1600 }
1601 }
1602 else{
1603 if (cont>=maxcont){
1604 nusefulbins++;
1605 binsProfile[nusefulbins-1]=h->GetBinLowEdge(startBin);
1606 binsProfile[nusefulbins]=max;
1607 }
1608 else {
1609 binsProfile[nusefulbins]=h->GetBinLowEdge(startBin);
1610 }
1611 }
1612 }
1613 return nusefulbins;
1614}