]> 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
changes to include AliPerformanceMach component (Michael Knichel)
[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$
ea23eb9f 17Revision 1.6 2007/10/05 10:44:11 zampolli
18Oversight for debugging purpose fixed
19
b8047c99 20Revision 1.5 2007/10/05 10:38:10 zampolli
21Oversight fixed
22
1e76f1ec 23Revision 1.4 2007/10/04 15:36:44 zampolli
24Updates to new TOF offline calibration schema
25
e26353e8 26Revision 1.3 2007/07/31 07:26:16 zampolli
27Bug fixed in the $ field
28
1ed77580 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
1ed77580 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>
e26353e8 49#include <TGrid.h>
4db98a6a 50#include <TList.h>
51#include <TArrayF.h>
52#include <TBenchmark.h>
1ed77580 53
54#include "AliTOFCalibTask.h"
70ddd0ee 55#include "AliESDEvent.h"
1ed77580 56#include "AliESDtrack.h"
57#include "AliLog.h"
4db98a6a 58#include "AliTOFArray.h"
e26353e8 59
4db98a6a 60//_________________________________________________________________________
1ed77580 61AliTOFCalibTask::AliTOFCalibTask(const char *name) :
4db98a6a 62 AliAnalysisTaskSE(name),
1ed77580 63 fESD(0),
1ed77580 64 fToT(0),
65 fTime(0),
66 fExpTimePi(0),
67 fExpTimeKa(0),
68 fExpTimePr(0),
1ed77580 69 fMinTime(0),
4db98a6a 70 fTOFArray(0x0),
1ed77580 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),
4db98a6a 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)
1ed77580 94{
95 // Constructor.
1ed77580 96
4db98a6a 97 DefineOutput(1,TList::Class());
98 DefineOutput(2,TList::Class());
70ddd0ee 99
e26353e8 100 for (Int_t i=0;i<11;i++){
101 fassparticle[i]=-1;
102 }
103
1ed77580 104}
105
106//______________________________________________________________________________
107AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
4db98a6a 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
1ed77580 141{
142 // Copy Constructor.
70ddd0ee 143
e26353e8 144 for (Int_t i=0;i<11;i++){
145 fassparticle[i]=calibtask.fassparticle[i];
146 }
147
1ed77580 148}
149//______________________________________________________________________________
150AliTOFCalibTask:: ~AliTOFCalibTask()
151{
152 // destructor
153
154 AliInfo("TOF Calib Task: Deleting");
4db98a6a 155
156 delete fTOFArray;
1ed77580 157 delete fhToT;
158 delete fhTime;
159 delete fhExpTimePi;
160 delete fhExpTimeKa;
161 delete fhExpTimePr;
162 delete fhPID;
163 delete fhch;
4db98a6a 164 delete fhESD;
165 delete fhESDselected;
166 delete fhESDkTOFout;
167 delete fhESDkTIME;
168 delete fhESDTRDcut;
169 delete fhESDTIMEcut;
170 delete fhESDassTOFcl;
171 // delete fListOfHistos;
1ed77580 172}
173//______________________________________________________________________________
174AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
175{
176 //assignment operator
4db98a6a 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;
e26353e8 208
209 for (Int_t i=0;i<11;i++){
4db98a6a 210 fassparticle[i]=calibtask.fassparticle[i];
e26353e8 211 }
4db98a6a 212
1ed77580 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
1ed77580 237
4db98a6a 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
1ed77580 275}
276//----------------------------------------------------------------------------
277void AliTOFCalibTask::DrawHistos(){
278
279 // drawing output histos
280
281 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
282
4db98a6a 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()));
1ed77580 373 }
4db98a6a 374 return;
70ddd0ee 375}
376
1ed77580 377//________________________________________________________________________
4db98a6a 378void AliTOFCalibTask::UserCreateOutputObjects()
1ed77580 379{
4db98a6a 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;
1ed77580 394}
395
396//______________________________________________________________________________
4db98a6a 397void AliTOFCalibTask::UserExec(Option_t * /*opt*/)
1ed77580 398{
399
400 // main
401
1ed77580 402
403//******* The loop over events -----------------------------------------------
404
4db98a6a 405 AliVEvent* fESD = fInputEvent ;
1ed77580 406 if (!fESD) {
4db98a6a 407 Error("UserExec","NO EVENT FOUND!");
408 return;
1ed77580 409 }
410
1ed77580 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-- ) {
4db98a6a 418 fhESD->Fill(0);
1ed77580 419 itr++;
4db98a6a 420 AliESDtrack * t = (AliESDtrack*)fESD->GetTrack(ntrk) ;
1ed77580 421 //selecting only good quality tracks
422 if (!Select(t)) continue;
423 nselected++;
4db98a6a 424 fhESDselected->Fill(0);
1ed77580 425 Int_t ich = Int_t(t->GetTOFCalChannel());
426 fhch->Fill(ich);
b8047c99 427 // ich=3; //only for debug purpose
1ed77580 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();
e26353e8 431 AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time));
1ed77580 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;
4db98a6a 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));
1ed77580 440 continue;
441 }
4db98a6a 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
1ed77580 458 }
459 fnESDselected+=nselected;
4db98a6a 460 PostData(1, fListOfHistos);
461 PostData(2, fListArray);
1ed77580 462
1ed77580 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");
4db98a6a 474 DrawHistos();
475 fListArray = dynamic_cast<TList*>(GetOutputData(2));
476 fTOFArray = dynamic_cast<AliTOFArray*>(fListArray->At(0));
477
e26353e8 478 for (Int_t ich = 0;ich<TOFCHANNELS;ich++){
4db98a6a 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 }
e26353e8 484 }
485
4db98a6a 486 TBenchmark bench;
487 bench.Start("CombPID");
1ed77580 488 for (Int_t i = 0;i<TOFCHANNELS;i++){
4db98a6a 489 Int_t size=fTOFArray->GetArraySize(i);
490 AliDebug(2, Form(" entries %i in channel %i ",size/NIDX,i));
491 if (size/NIDX<=2) {
1ed77580 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 }
4db98a6a 495 if (i%1000 == 0) AliInfo(Form("At channel %d",i));
496 if (!CombPID(fTOFArray->GetArray(i), size)) AliError("ERROR!!!!ERROR!!!");
1ed77580 497 }
4db98a6a 498 bench.Stop("CombPID");
499 bench.Print("CombPID");
1ed77580 500
4db98a6a 501 // saving data in a tree --> obsolete code; keeping for backup with new structure
502 // using AliTOFArray
503 /*
e26353e8 504 AliInfo("Building tree for Calibration");
505 TTree * tree = new TTree("T", "Tree for TOF calibration");
4db98a6a 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 }
1ed77580 519 }
4db98a6a 520 tree->Fill();
521 */
1ed77580 522
4db98a6a 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;
1ed77580 529}
1ed77580 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++;
4db98a6a 540 fhESDkTOFout->Fill(0);
1ed77580 541 //IsStartedTimeIntegral
542 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
543 return 0;
544 }
545 fnESDkTIME++;
4db98a6a 546 fhESDkTIME->Fill(0);
1ed77580 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++;
4db98a6a 552 fhESDTRDcut->Fill(0);
1ed77580 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++;
4db98a6a 559 fhESDTIMEcut->Fill(0);
1ed77580 560
561 Double_t mom=t->GetP();
562 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
563 // return 0; // skipping momentum cut
564 }
565
21a8ed8d 566 Int_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
567 if(assignedTOFcluster==-1){ // not matched
1ed77580 568 return 0;
569 }
570 fnESDassTOFcl++;
4db98a6a 571 fhESDassTOFcl->Fill(0);
572 AliDebug(2,"selecting the track");
1ed77580 573 return 1;
574}
575//_____________________________________________________________________________
576
e26353e8 577Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
578
4db98a6a 579 // to be re-implemented with new object
580
e26353e8 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));
4db98a6a 591 /*
592 Int_t nTracksInChannel = fTOFArray->GetArray(ich)->GetSize();
e26353e8 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));
4db98a6a 599 nTracksInChannel;
e26353e8 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 }
4db98a6a 621 */
e26353e8 622 return ndeleted;
623}
624//_____________________________________________________________________________
625
1ed77580 626Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
627
e26353e8 628 // track Combinatorial PID for calibration, only when more than 2 particles
629 // fall in channel
1ed77580 630
1ed77580 631 Float_t t0offset=0;
e26353e8 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!");
1ed77580 635 return 0;
636 }
e26353e8 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 }
1ed77580 647 }
648
4db98a6a 649 AliDebug(2,Form("number of sets = %i", nset));
e26353e8 650 // loop on sets
1ed77580 651 for (Int_t i=0; i< nset; i++) {
e26353e8 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 }
4db98a6a 661 AliDebug(2,Form("set = %i, number of tracks in set = %i", i,ntrkinset));
e26353e8 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));
1ed77580 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 }
e26353e8 695
696
1ed77580 697 Float_t chisquarebest=999.;
4db98a6a 698 AliDebug(2,Form(" Set = %i with %i tracks ", i,ntrkinset));
e26353e8 699 chisquarebest = LoopCombPID(ntrkinset, ntrkinset,exptof,&texp[0],&timeofflight[0], &index[0],chisquarebest);
700
1ed77580 701 Float_t confLevel=999;
702 if(chisquarebest<999.){
703 Double_t dblechisquare=(Double_t)chisquarebest;
e26353e8 704 confLevel=(Float_t)TMath::Prob(dblechisquare,ntrkinset-1);
1ed77580 705 }
e26353e8 706
1ed77580 707 // assume they are all pions for fake sets
708 if(confLevel<0.01 || confLevel==999. ){
e26353e8 709 for (Int_t itrk=0; itrk<ntrkinset; itrk++)fassparticle[itrk]=0;
1ed77580 710 }
e26353e8 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]));
1ed77580 713
e26353e8 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;
1ed77580 718 // storing in third slot of smallarray the assigned expected time
e26353e8 719 fhPID->Fill(fassparticle[kk]);
720 // fassparticle[kk]=0; //assuming all particles are pions
721 if (fassparticle[kk]==0){ //assigned to be a Pi
1ed77580 722 smallarray[idxextimeKa]=-1;
723 smallarray[idxextimePr]=-1;
724 }
e26353e8 725 else if (fassparticle[kk]==1){ //assigned to be a Ka
1ed77580 726 smallarray[idxextimePi]=smallarray[idxextimeKa];
727 smallarray[idxextimeKa]=-1;
728 smallarray[idxextimePr]=-1;
729 }
e26353e8 730 else if (fassparticle[kk]==2){ //assigned to be a Pr
1ed77580 731 smallarray[idxextimePi]=smallarray[idxextimePr];
732 smallarray[idxextimeKa]=-1;
733 smallarray[idxextimePr]=-1;
734 }
735 }
ea23eb9f 736 delete[] exptof;
737 delete[] timeofflight;
738 delete[] texp;
739 delete[] index;
e26353e8 740
1ed77580 741 }
742 return 1;
743}
e26353e8 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
4db98a6a 747 // performing combinatorial PID in recursive way
748
e26353e8 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);
1ed77580 756 }
e26353e8 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++;
1ed77580 787 }
e26353e8 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;
1ed77580 801 }
ea23eb9f 802 delete[] besttimezero;
803 delete[] bestchisquare;
804 delete[] bestweightedtimezero;
805 delete[] weightedtimezero;
806 delete[] timezero;
1ed77580 807 }
808 }
e26353e8 809 return chisquarebest;
1ed77580 810}