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