]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TOF/AliTOFCalibTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Coverity RESOURCE_LEAK fixes
[u/mrichter/AliRoot.git] / TOF / AliTOFCalibTask.cxx
... / ...
CommitLineData
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$
17Revision 1.6 2007/10/05 10:44:11 zampolli
18Oversight for debugging purpose fixed
19
20Revision 1.5 2007/10/05 10:38:10 zampolli
21Oversight fixed
22
23Revision 1.4 2007/10/04 15:36:44 zampolli
24Updates to new TOF offline calibration schema
25
26Revision 1.3 2007/07/31 07:26:16 zampolli
27Bug fixed in the $ field
28
29*/
30
31// task to perform TOF calibration
32// optimized for pp runs
33// expect a maximum of 100 entries per channel
34// different ways to calibrate are implemented
35// C. Zampolli
36
37#include <TObject.h>
38#include <TCanvas.h>
39#include <TROOT.h>
40#include <TSystem.h>
41#include <TProfile.h>
42#include <TF1.h>
43#include <TTree.h>
44#include <TFile.h>
45#include <TH1F.h>
46#include <TH1I.h>
47#include <TH1D.h>
48#include <TH2F.h>
49#include <TGrid.h>
50#include <TList.h>
51#include <TArrayF.h>
52#include <TBenchmark.h>
53
54#include "AliTOFCalibTask.h"
55#include "AliESDEvent.h"
56#include "AliESDtrack.h"
57#include "AliLog.h"
58#include "AliTOFArray.h"
59
60//_________________________________________________________________________
61AliTOFCalibTask::AliTOFCalibTask(const char *name) :
62 AliAnalysisTaskSE(name),
63 fESD(0),
64 fToT(0),
65 fTime(0),
66 fExpTimePi(0),
67 fExpTimeKa(0),
68 fExpTimePr(0),
69 fMinTime(0),
70 fTOFArray(0x0),
71 fnESD(0),
72 fnESDselected(0),
73 fnESDkTOFout(0),
74 fnESDkTIME(0),
75 fnESDassTOFcl(0),
76 fnESDTIMEcut(0),
77 fnESDTRDcut(0),
78 fhToT(0),
79 fhTime(0),
80 fhExpTimePi(0),
81 fhExpTimeKa(0),
82 fhExpTimePr(0),
83 fhPID(0),
84 fhch(0),
85 fhESD(0),
86 fhESDselected(0),
87 fhESDkTOFout(0),
88 fhESDkTIME(0),
89 fhESDassTOFcl(0),
90 fhESDTIMEcut(0),
91 fhESDTRDcut(0),
92 fListOfHistos(0x0),
93 fListArray(0x0)
94{
95 // Constructor.
96
97 DefineOutput(1,TList::Class());
98 DefineOutput(2,TList::Class());
99
100 for (Int_t i=0;i<11;i++){
101 fassparticle[i]=-1;
102 }
103
104}
105
106//______________________________________________________________________________
107AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
108 AliAnalysisTaskSE("AliTOFCalibTask"),
109 fESD(calibtask.fESD),
110 fToT(calibtask.fToT),
111 fTime(calibtask.fTime),
112 fExpTimePi(calibtask.fExpTimePi),
113 fExpTimeKa(calibtask.fExpTimeKa),
114 fExpTimePr(calibtask.fExpTimePr),
115 fMinTime(calibtask.fMinTime),
116 fTOFArray(calibtask.fTOFArray),
117 fnESD(calibtask.fnESD),
118 fnESDselected(calibtask.fnESDselected),
119 fnESDkTOFout(calibtask.fnESDkTOFout),
120 fnESDkTIME(calibtask.fnESDkTIME),
121 fnESDassTOFcl(calibtask.fnESDassTOFcl),
122 fnESDTIMEcutcalibtask.fnESDTIMEcut(),
123 fnESDTRDcutcalibtask.fnESDTRDcut(),
124 fhToT(calibtask.fhToT),
125 fhTime(calibtask.fhTime),
126 fhExpTimePi(calibtask.fhExpTimePi),
127 fhExpTimeKa(calibtask.fhExpTimeKa),
128 fhExpTimePr(calibtask.fhExpTimePr),
129 fhPID(calibtask.fhPID),
130 fhch(calibtask.fhch),
131 fhESD(calibtask.fhESD),
132 fhESDselected(calibtask.fhESDselected),
133 fhESDkTOFout(calibtask.fhESDkTOFout),
134 fhESDkTIME(calibtask.fhESDkTIME),
135 fhESDassTOFcl(calibtask.fhESDassTOFcl),
136 fhESDTIMEcut(calibtask.fhESDTIMEcut),
137 fhESDTRDcut(calibtask.fhESDTRDcut),
138 fListOfHistos(calibtask.fListOfHistos),
139 fListArray(calibtask.fListArray)
140
141{
142 // Copy Constructor.
143
144 for (Int_t i=0;i<11;i++){
145 fassparticle[i]=calibtask.fassparticle[i];
146 }
147
148}
149//______________________________________________________________________________
150AliTOFCalibTask:: ~AliTOFCalibTask()
151{
152 // destructor
153
154 AliInfo("TOF Calib Task: Deleting");
155
156 delete fTOFArray;
157 delete fhToT;
158 delete fhTime;
159 delete fhExpTimePi;
160 delete fhExpTimeKa;
161 delete fhExpTimePr;
162 delete fhPID;
163 delete fhch;
164 delete fhESD;
165 delete fhESDselected;
166 delete fhESDkTOFout;
167 delete fhESDkTIME;
168 delete fhESDTRDcut;
169 delete fhESDTIMEcut;
170 delete fhESDassTOFcl;
171 // delete fListOfHistos;
172}
173//______________________________________________________________________________
174AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
175{
176 //assignment operator
177 fESD=calibtask.fESD;
178 fToT=calibtask.fToT;
179 fTime=calibtask.fTime;
180 fExpTimePi=calibtask.fExpTimePi;
181 fExpTimeKa=calibtask.fExpTimeKa;
182 fExpTimePr=calibtask.fExpTimePr;
183 fMinTime=calibtask.fMinTime;
184 fTOFArray=calibtask.fTOFArray;
185 fnESD=calibtask.fnESD;
186 fnESDselected=calibtask.fnESDselected;
187 fnESDkTOFout=calibtask.fnESDkTOFout;
188 fnESDkTIME=calibtask.fnESDkTIME;
189 fnESDassTOFcl=calibtask.fnESDassTOFcl;
190 fnESDTIMEcut=calibtask.fnESDTIMEcut;
191 fnESDTRDcut=calibtask.fnESDTRDcut;
192 fhToT=calibtask.fhToT;
193 fhTime=calibtask.fhTime;
194 fhExpTimePi=calibtask.fhExpTimePi;
195 fhExpTimeKa=calibtask.fhExpTimeKa;
196 fhExpTimePr=calibtask.fhExpTimePr;
197 fhPID=calibtask.fhPID;
198 fhch=calibtask.fhch;
199 fhESD=calibtask.fhESD;
200 fhESDselected=calibtask.fhESDselected;
201 fhESDkTOFout=calibtask.fhESDkTOFout;
202 fhESDkTIME=calibtask.fhESDkTIME;
203 fhESDassTOFcl=calibtask.fhESDassTOFcl;
204 fhESDTIMEcut=calibtask.fhESDTIMEcut;
205 fhESDTRDcut=calibtask.fhESDTRDcut;
206 fListOfHistos=calibtask.fListOfHistos;
207 fListArray=calibtask.fListArray;
208
209 for (Int_t i=0;i<11;i++){
210 fassparticle[i]=calibtask.fassparticle[i];
211 }
212
213 return *this;
214}
215//--------------------------------------------------------------------------
216void AliTOFCalibTask::BookHistos(){
217
218 //booking histos
219
220 AliInfo(Form("*** Booking Histograms %s", GetName())) ;
221
222 fhToT=
223 new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
224 fhTime=
225 new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
226 fhExpTimePi=
227 new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
228 fhExpTimeKa=
229 new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
230 fhExpTimePr=
231 new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
232 fhPID=
233 new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
234 fhch=
235 new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
236
237
238 // create the output list of histos
239 if (!fListOfHistos) fListOfHistos = new TList();
240 fListOfHistos->SetOwner();
241
242 fListOfHistos->AddAt(fhToT, 0) ;
243 fListOfHistos->AddAt(fhTime, 1) ;
244 fListOfHistos->AddAt(fhExpTimePi, 2) ;
245 fListOfHistos->AddAt(fhExpTimeKa, 3) ;
246 fListOfHistos->AddAt(fhExpTimePr, 4) ;
247 fListOfHistos->AddAt(fhPID, 5) ;
248 fListOfHistos->AddAt(fhch, 6) ;
249
250 fhESD=
251 new TH1I("hESD","Number of analyzed ESDs",1,0,1);
252 fhESDselected=
253 new TH1I("hESDselected","Number of selected ESDs",1,0,1);
254 fhESDkTOFout=
255 new TH1I("hESDkTOFout","Number of ESDs with kTOFout",1,0,1);
256 fhESDkTIME=
257 new TH1I("hESDkTIME","Number of ESDs with kTime",1,0,1);
258 fhESDTRDcut=
259 new TH1I("hESDTRDcut","Number of ESDs with TRDcut",1,0,1);
260 fhESDTIMEcut=
261 new TH1I("hESDTIMEcut","Number of ESDs with TIMEcut",1,0,1);
262 fhESDassTOFcl=
263 new TH1I("hESDassTOFcl","Number of ESDs with assTOFcl",1,0,1);
264
265 fListOfHistos->AddAt(fhESD, 7) ;
266 fListOfHistos->AddAt(fhESDselected, 8) ;
267 fListOfHistos->AddAt(fhESDkTOFout, 9) ;
268 fListOfHistos->AddAt(fhESDkTIME, 10) ;
269 fListOfHistos->AddAt(fhESDTRDcut, 11) ;
270 fListOfHistos->AddAt(fhESDTIMEcut, 12) ;
271 fListOfHistos->AddAt(fhESDassTOFcl, 13) ;
272
273 return;
274
275}
276//----------------------------------------------------------------------------
277void AliTOFCalibTask::DrawHistos(){
278
279 // drawing output histos
280
281 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
282
283 if (!gROOT->IsBatch()){
284
285 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
286
287 fhToT = (TH1F*)fListOfHistos->At(0);
288 fhTime = (TH1F*)fListOfHistos->At(1);
289 fhExpTimePi = (TH1F*)fListOfHistos->At(2);
290 fhExpTimeKa = (TH1F*)fListOfHistos->At(3);
291 fhExpTimePr = (TH1F*)fListOfHistos->At(4);
292 fhPID = (TH1I*)fListOfHistos->At(5);
293 fhch = (TH1D*)fListOfHistos->At(6);
294
295 fhESD = (TH1I*)fListOfHistos->At(7);
296 fhESDselected = (TH1I*)fListOfHistos->At(8);
297 fhESDkTOFout = (TH1I*)fListOfHistos->At(9);
298 fhESDkTIME = (TH1I*)fListOfHistos->At(10);
299 fhESDassTOFcl = (TH1I*)fListOfHistos->At(11);
300 fhESDTIMEcut = (TH1I*)fListOfHistos->At(12);
301 fhESDTRDcut = (TH1I*)fListOfHistos->At(13);
302
303 TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
304 canvasToTTime->Divide(1,2);
305 canvasToTTime->cd(1);
306 fhToT->SetLineColor(4);
307 fhToT->GetXaxis()->SetTitle("ToT (ns)");
308 fhToT->DrawCopy("hist");
309 canvasToTTime->cd(2);
310 fhTime->SetLineColor(4);
311 fhTime->GetXaxis()->SetTitle("Time (ns)");
312 fhTime->DrawCopy("hist");
313 canvasToTTime->Update();
314 //canvasToTTime->Print("ToTTime.gif");
315
316 TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
317 canvasExpTime->Divide(1,3);
318 canvasExpTime->cd(1);
319 fhExpTimePi->SetLineColor(4);
320 fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
321 fhExpTimePi->DrawCopy("hist");
322 canvasExpTime->cd(2);
323 fhExpTimeKa->SetLineColor(4);
324 fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
325 fhExpTimeKa->DrawCopy("hist");
326 canvasExpTime->cd(3);
327 fhExpTimePr->SetLineColor(4);
328 fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
329 fhExpTimePr->DrawCopy("hist");
330
331 //canvasExpTime->Print("ExpTime.gif");
332
333 TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
334 canvasPID->cd();
335 fhPID->GetXaxis()->SetTitle("Comb PID");
336 fhPID->GetXaxis()->SetBinLabel(1,"#pi");
337 fhPID->GetXaxis()->SetBinLabel(2,"K");
338 fhPID->GetXaxis()->SetBinLabel(3,"p");
339 fhPID->DrawCopy("hist");
340
341 //canvasPID->Print("PID.gif");
342
343 TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
344 canvasrndch->cd();
345 fhch->GetXaxis()->SetTitle("TOF ch");
346 fhch->Draw("hist");
347 Float_t meanTOFch = 0;
348 for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
349 meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
350 }
351
352 meanTOFch/=TOFCHANNELS;
353 AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
354
355 //canvasrndch->Print("rndch.gif");
356
357 /*
358 char line[1024] ;
359 sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ;
360 gROOT->ProcessLine(line);
361 sprintf(line, ".!rm -fR *.gif");
362 gROOT->ProcessLine(line);
363 AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
364 */
365
366 AliInfo(Form(" Number of analyzed ESD tracks: %i",(Int_t)fhESD->GetEntries()));
367 AliInfo(Form(" Number of selected ESD tracks: %i",(Int_t)fhESDselected->GetEntries()));
368 AliInfo(Form(" Number of ESD tracks with kTOFout: %i",(Int_t)fhESDkTOFout->GetEntries()));
369 AliInfo(Form(" Number of ESD tracks with kTIME: %i",(Int_t)fhESDkTIME->GetEntries()));
370 AliInfo(Form(" Number of ESD tracks with TRDcut: %i",(Int_t)fhESDTRDcut->GetEntries()));
371 AliInfo(Form(" Number of ESD tracks with TIMEcut: %i",(Int_t)fhESDTIMEcut->GetEntries()));
372 AliInfo(Form(" Number of ESD tracks with assTOFcl: %i",(Int_t)fhESDassTOFcl->GetEntries()));
373 }
374 return;
375}
376
377//________________________________________________________________________
378void AliTOFCalibTask::UserCreateOutputObjects()
379{
380// Create histograms and AliTOFArray
381
382 AliInfo("UserCreateObjects");
383 OpenFile(1);
384 BookHistos();
385
386 if (!fTOFArray) fTOFArray = new AliTOFArray(TOFCHANNELS);
387 for (Int_t i =0;i<TOFCHANNELS;i++){
388 fTOFArray->SetArray(i);
389 }
390 if (!fListArray) fListArray = new TList();
391 fListArray->SetOwner();
392 fListArray->Add(fTOFArray);
393 return;
394}
395
396//______________________________________________________________________________
397void AliTOFCalibTask::UserExec(Option_t * /*opt*/)
398{
399
400 // main
401
402
403//******* The loop over events -----------------------------------------------
404
405 AliVEvent* fESD = fInputEvent ;
406 if (!fESD) {
407 Error("UserExec","NO EVENT FOUND!");
408 return;
409 }
410
411
412 fMinTime=22E3; //ns; not used
413 Int_t ntrk = fESD->GetNumberOfTracks() ;
414 fnESD+=ntrk;
415 Int_t nselected = 0;
416 Int_t itr = -1;
417 while ( ntrk-- ) {
418 fhESD->Fill(0);
419 itr++;
420 AliESDtrack * t = (AliESDtrack*)fESD->GetTrack(ntrk) ;
421 //selecting only good quality tracks
422 if (!Select(t)) continue;
423 nselected++;
424 fhESDselected->Fill(0);
425 Int_t ich = Int_t(t->GetTOFCalChannel());
426 fhch->Fill(ich);
427 // ich=3; //only for debug purpose
428 AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
429 Float_t tot = t->GetTOFsignalToT();
430 Float_t time = t->GetTOFsignalRaw();
431 AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time));
432 Double_t expTime[10];
433 t->GetIntegratedTimes(expTime);
434 Float_t expTimePi = expTime[2]*1.E-3;
435 Float_t expTimeKa = expTime[3]*1.E-3;
436 Float_t expTimePr = expTime[4]*1.E-3;
437 Int_t currentSize = fTOFArray->GetArraySize(ich);
438 if (currentSize==CHENTRIES){
439 AliDebug(2,Form("too many tracks in channel %i, not storing any more...",ich));
440 continue;
441 }
442 AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, currentSize/NIDX ));
443
444 fTOFArray->ReSetArraySize(ich,currentSize+5);
445 fTOFArray->SetAt(ich,currentSize+DELTAIDXTOT,tot);
446 fTOFArray->SetAt(ich,currentSize+DELTAIDXTIME,time*1E-3);
447 fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEPI,expTimePi);
448 fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEKA,expTimeKa);
449 fTOFArray->SetAt(ich,currentSize+DELTAIDXEXTIMEPR,expTimePr);
450 fhToT->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTOT));
451 fhTime->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTIME));
452 fhExpTimePi->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEPI));
453 fhExpTimeKa->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEKA));
454 fhExpTimePr->Fill(fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXEXTIMEPR));
455 AliDebug(2,Form("track = %i, tot = %f, time = %f, and Exp time in TOF: pi = %f, K = %f, p = %f",itr, fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTOT), fTOFArray->GetArrayAt(ich,currentSize+DELTAIDXTIME), expTimePi,expTimeKa,expTimePr));
456
457
458 }
459 fnESDselected+=nselected;
460 PostData(1, fListOfHistos);
461 PostData(2, fListArray);
462
463}
464
465//_____________________________________________________________________________
466void AliTOFCalibTask::Terminate(Option_t *)
467{
468 // Processing when the event loop is ended
469
470 // some plots
471
472 TH1::AddDirectory(0);
473 AliInfo("TOF Calib Task: End of events loop");
474 DrawHistos();
475 fListArray = dynamic_cast<TList*>(GetOutputData(2));
476 fTOFArray = dynamic_cast<AliTOFArray*>(fListArray->At(0));
477
478 for (Int_t ich = 0;ich<TOFCHANNELS;ich++){
479 Int_t entriesChannel = fTOFArray->GetArraySize(ich)/NIDX;
480 if (entriesChannel>0){
481 Int_t ncuttime = SelectOnTime(fTOFArray->GetArray(ich),entriesChannel,ich);
482 fnESDselected-=ncuttime;
483 }
484 }
485
486 TBenchmark bench;
487 bench.Start("CombPID");
488 for (Int_t i = 0;i<TOFCHANNELS;i++){
489 Int_t size=fTOFArray->GetArraySize(i);
490 AliDebug(2, Form(" entries %i in channel %i ",size/NIDX,i));
491 if (size/NIDX<=2) {
492 AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
493 continue;
494 }
495 if (i%1000 == 0) AliInfo(Form("At channel %d",i));
496 if (!CombPID(fTOFArray->GetArray(i), size)) AliError("ERROR!!!!ERROR!!!");
497 }
498 bench.Stop("CombPID");
499 bench.Print("CombPID");
500
501 // saving data in a tree --> obsolete code; keeping for backup with new structure
502 // using AliTOFArray
503 /*
504 AliInfo("Building tree for Calibration");
505 TTree * tree = new TTree("T", "Tree for TOF calibration");
506 AliTOFArray* tempArray = new AliTOFArray(TOFCHANNELS);
507 tree->Branch("TOFentries","AliTOFArray",&tempArray,32000,0);
508 for (Int_t i = 0;i<TOFCHANNELS;i++){
509 Int_t sizeChannel = fTOFArray->GetArraySize(i)/NIDX;
510 if (i==3) AliDebug(2,Form("Entries in channel %d = %d ",i,sizeChannel)); // just for debug
511 if (sizeChannel!=0){
512 tempArray->SetArray(i,NIDXSMALL*sizeChannel);
513 }
514 for (Int_t j =0; j<sizeChannel;j++){
515 for (Int_t k=0; k<NIDXSMALL;k++){
516 tempArray->SetAt(i,j*NIDXSMALL+k,fTOFArray->GetArrayAt(i,j*NIDX+k));
517 }
518 }
519 }
520 tree->Fill();
521 */
522
523 AliInfo("Putting object for calibration in Reference data");
524
525 TFile *fileout = TFile::Open("outputTOFCalibration.root","RECREATE");
526 tempArray->Write();
527 fileout->Close();
528 delete tempArray;
529}
530//_____________________________________________________________________________
531
532Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
533
534 // selecting good quality tracks
535 //track selection: reconstrution to TOF:
536 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
537 return 0;
538 }
539 fnESDkTOFout++;
540 fhESDkTOFout->Fill(0);
541 //IsStartedTimeIntegral
542 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
543 return 0;
544 }
545 fnESDkTIME++;
546 fhESDkTIME->Fill(0);
547 if (t->GetStatus() & AliESDtrack::kTRDbackup) {
548 Float_t xout= t->GetOuterParam()->GetX();
549 if (xout<364.25 && xout > 300.) return 0;
550 }
551 fnESDTRDcut++;
552 fhESDTRDcut->Fill(0);
553 Double_t time=t->GetTOFsignal();
554 time*=1.E-3; // tof given in nanoseconds
555 if(time >= fMinTime){
556 return 0;
557 }
558 fnESDTIMEcut++;
559 fhESDTIMEcut->Fill(0);
560
561 Double_t mom=t->GetP();
562 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
563 // return 0; // skipping momentum cut
564 }
565
566 Int_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
567 if(assignedTOFcluster==-1){ // not matched
568 return 0;
569 }
570 fnESDassTOFcl++;
571 fhESDassTOFcl->Fill(0);
572 AliDebug(2,"selecting the track");
573 return 1;
574}
575//_____________________________________________________________________________
576
577Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
578
579 // to be re-implemented with new object
580
581 // discarding tracks with time-mintime < MINTIME
582
583 Int_t ndeleted=0;
584 Float_t mintime = 1E6;
585 for (Int_t itr = 0;itr<ntracks;itr++){
586 Int_t ientry=itr*NIDX;
587 Float_t time = charray[ientry+DELTAIDXTIME];// in ns
588 if (time<mintime) mintime = time;
589 }
590 AliDebug(1,Form("Mintime for channel %i = %f",ich, mintime));
591 /*
592 Int_t nTracksInChannel = fTOFArray->GetArray(ich)->GetSize();
593 for (Int_t itr = 0;itr<ntracks;itr++){
594 Int_t ientry=itr*NIDX;
595 Float_t time = charray[ientry+DELTAIDXTIME];// in ns
596 if ((time-mintime)>MINTIME) {
597 ndeleted++;
598 AliDebug(1,Form("Deleting %i track from channel %i, time = %f",ndeleted, ich, time));
599 nTracksInChannel;
600 for (Int_t j=itr+1;j<ntracks;j++){
601 Int_t ientry=j*NIDX;
602 Int_t idxtot = ientry+DELTAIDXTOT;
603 Int_t idxtime = ientry+DELTAIDXTIME;
604 Int_t idxextimePi = ientry+DELTAIDXEXTIMEPI;
605 Int_t idxextimeKa = ientry+DELTAIDXEXTIMEKA;
606 Int_t idxextimePr = ientry+DELTAIDXEXTIMEPR;
607 Int_t ientrydel=(j-1)*NIDX;
608 Int_t idxtotdel = ientrydel+DELTAIDXTOT;
609 Int_t idxtimedel = ientrydel+DELTAIDXTIME;
610 Int_t idxextimePidel = ientrydel+DELTAIDXEXTIMEPI;
611 Int_t idxextimeKadel = ientrydel+DELTAIDXEXTIMEKA;
612 Int_t idxextimePrdel = ientrydel+DELTAIDXEXTIMEPR;
613 charray[idxtotdel]=charray[idxtot];
614 charray[idxtimedel]=charray[idxtime];
615 charray[idxextimePidel]=charray[idxextimePi];
616 charray[idxextimeKadel]=charray[idxextimeKa];
617 charray[idxextimePrdel]=charray[idxextimePr];
618 }
619 }
620 }
621 */
622 return ndeleted;
623}
624//_____________________________________________________________________________
625
626Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
627
628 // track Combinatorial PID for calibration, only when more than 2 particles
629 // fall in channel
630
631 Float_t t0offset=0;
632
633 if (size/NIDX<=2){
634 AliDebug(1,"Number of tracks in channel smaller than 2, identifying every particle as if it was a pion!");
635 return 0;
636 }
637
638 Int_t ntracksinchannel = size/NIDX;
639 Int_t ntrkinset=-1;
640 Int_t nset = ntracksinchannel/6;
641
642 if (nset ==0) {
643 nset=1;
644 if (ntracksinchannel < 6) {
645 AliDebug(2,"Number of tracks in set smaller than 6, Combinatorial PID still applied to this set.");
646 }
647 }
648
649 AliDebug(2,Form("number of sets = %i", nset));
650 // loop on sets
651 for (Int_t i=0; i< nset; i++) {
652 if (i<nset-1)ntrkinset=6;
653 else {
654 if (ntracksinchannel<6){
655 ntrkinset=ntracksinchannel;
656 }
657 else{
658 ntrkinset=6+Int_t((ntracksinchannel)%6);
659 }
660 }
661 AliDebug(2,Form("set = %i, number of tracks in set = %i", i,ntrkinset));
662
663 for (Int_t ii=0;ii<ntrkinset;ii++){
664 fassparticle[ii]=-1;
665 }
666 Float_t **exptof;
667 Float_t* timeofflight;
668 Float_t* texp;
669 Int_t* index;
670
671 exptof = new Float_t*[ntrkinset];
672 timeofflight = new Float_t[ntrkinset];
673 texp = new Float_t[ntrkinset];
674 index = new Int_t[ntrkinset];
675 for (Int_t ii=0;ii<ntrkinset;ii++){
676 exptof[ii] = new Float_t[3];
677 timeofflight[ii]=0;
678 texp[ii]=0;
679 index[ii]=-1;
680 }
681
682 for (Int_t j=0; j<ntrkinset; j++) {
683 Int_t idxtime = ((6*i+j)*NIDX)+DELTAIDXTIME;
684 Int_t idxextimePi = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPI;
685 Int_t idxextimeKa = ((6*i+j)*NIDX)+DELTAIDXEXTIMEKA;
686 Int_t idxextimePr = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPR;
687 AliDebug(2,Form("idxtime = %i, idxextimePi = %i, idxextimeKa = %i, idxextimePr = %i", idxtime, idxextimePi, idxextimeKa, idxextimePr));
688 Double_t time=smallarray[idxtime]; // TOF time in ns
689 timeofflight[j]=time+t0offset;
690 exptof[j][0]=smallarray[idxextimePi];
691 exptof[j][1]=smallarray[idxextimeKa];
692 exptof[j][2]=smallarray[idxextimePr];
693 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]));
694 }
695
696
697 Float_t chisquarebest=999.;
698 AliDebug(2,Form(" Set = %i with %i tracks ", i,ntrkinset));
699 chisquarebest = LoopCombPID(ntrkinset, ntrkinset,exptof,&texp[0],&timeofflight[0], &index[0],chisquarebest);
700
701 Float_t confLevel=999;
702 if(chisquarebest<999.){
703 Double_t dblechisquare=(Double_t)chisquarebest;
704 confLevel=(Float_t)TMath::Prob(dblechisquare,ntrkinset-1);
705 }
706
707 // assume they are all pions for fake sets
708 if(confLevel<0.01 || confLevel==999. ){
709 for (Int_t itrk=0; itrk<ntrkinset; itrk++)fassparticle[itrk]=0;
710 }
711
712 AliDebug(2,Form(" Best Assignment, selection %i %i %i %i %i %i %i %i %i %i",fassparticle[0],fassparticle[1],fassparticle[2],fassparticle[3],fassparticle[4],fassparticle[5],fassparticle[6],fassparticle[7],fassparticle[8],fassparticle[9]));
713
714 for (Int_t kk=0;kk<ntrkinset;kk++){
715 Int_t idxextimePi = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPI;
716 Int_t idxextimeKa = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEKA;
717 Int_t idxextimePr = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPR;
718 // storing in third slot of smallarray the assigned expected time
719 fhPID->Fill(fassparticle[kk]);
720 // fassparticle[kk]=0; //assuming all particles are pions
721 if (fassparticle[kk]==0){ //assigned to be a Pi
722 smallarray[idxextimeKa]=-1;
723 smallarray[idxextimePr]=-1;
724 }
725 else if (fassparticle[kk]==1){ //assigned to be a Ka
726 smallarray[idxextimePi]=smallarray[idxextimeKa];
727 smallarray[idxextimeKa]=-1;
728 smallarray[idxextimePr]=-1;
729 }
730 else if (fassparticle[kk]==2){ //assigned to be a Pr
731 smallarray[idxextimePi]=smallarray[idxextimePr];
732 smallarray[idxextimeKa]=-1;
733 smallarray[idxextimePr]=-1;
734 }
735 }
736 delete[] exptof;
737 delete[] timeofflight;
738 delete[] texp;
739 delete[] index;
740
741 }
742 return 1;
743}
744//---------------------------------------------------------------------
745Float_t AliTOFCalibTask::LoopCombPID(Int_t itrkinset, Int_t ntrkinset, Float_t **exptof, Float_t *texp, Float_t *timeofflight, Int_t *index, Float_t chisquarebest){
746
747 // performing combinatorial PID in recursive way
748
749 Int_t indextr = ntrkinset-itrkinset;
750
751 for (index[indextr]=0;index[indextr]<3;index[indextr]++){
752 Int_t ii = index[indextr];
753 texp[indextr]=exptof[indextr][ii];
754 if (indextr<ntrkinset-1){
755 chisquarebest = LoopCombPID(ntrkinset-indextr-1,ntrkinset,exptof,&texp[0],&timeofflight[0],&index[0], chisquarebest);
756 }
757
758 else {
759 Float_t *besttimezero = new Float_t[ntrkinset];
760 Float_t *bestchisquare = new Float_t[ntrkinset];
761 Float_t *bestweightedtimezero = new Float_t[ntrkinset];
762 Float_t sumAllweights=0.;
763 Float_t meantzero=0.;
764 Float_t *weightedtimezero = new Float_t[ntrkinset];
765 Float_t *timezero = new Float_t[ntrkinset];
766
767 AliDebug(2,Form("Configuration = %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, so far chisquarebest = %f ",index[0],index[1],index[2],index[3],index[4],index[5],index[6],index[7],index[8],index[9],chisquarebest));
768 for (Int_t itz=0; itz<ntrkinset;itz++) {
769 timezero[itz]=texp[itz]-timeofflight[itz];
770 weightedtimezero[itz]=timezero[itz]/TRACKERROR;
771 sumAllweights+=1./TRACKERROR;
772 meantzero+=weightedtimezero[itz];
773 } // end loop for (Int_t itz=0; itz<ntrkinset;itz++)
774
775 meantzero=meantzero/sumAllweights; // it is given in [ns]
776
777 // calculate chisquare
778
779 Float_t chisquare=0.;
780 for (Int_t icsq=0; icsq<ntrkinset;icsq++) {
781 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
782 } // end loop for (Int_t icsq=0; icsq<ntrkinset;icsq++)
783
784 Int_t npion=0;
785 for (Int_t j=0;j<ntrkinset;j++){
786 if(index[j]==0)npion++;
787 }
788
789 if(chisquare<=chisquarebest && ((Float_t) npion/ ((Float_t) ntrkinset)>0.3)){
790 for(Int_t iqsq = 0; iqsq<ntrkinset; iqsq++) {
791 besttimezero[iqsq]=timezero[iqsq];
792 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
793 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR;
794 }
795
796 for (Int_t j=0;j<ntrkinset;j++){
797 fassparticle[j]=index[j];
798 }
799
800 chisquarebest=chisquare;
801 }
802 delete[] besttimezero;
803 delete[] bestchisquare;
804 delete[] bestweightedtimezero;
805 delete[] weightedtimezero;
806 delete[] timezero;
807 }
808 }
809 return chisquarebest;
810}