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