1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 Revision 1.6 2007/10/05 10:44:11 zampolli
18 Oversight for debugging purpose fixed
20 Revision 1.5 2007/10/05 10:38:10 zampolli
23 Revision 1.4 2007/10/04 15:36:44 zampolli
24 Updates to new TOF offline calibration schema
26 Revision 1.3 2007/07/31 07:26:16 zampolli
27 Bug fixed in the $ field
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
52 #include <TBenchmark.h>
54 #include "AliTOFCalibTask.h"
55 #include "AliESDEvent.h"
56 #include "AliESDtrack.h"
58 #include "AliTOFArray.h"
60 //_________________________________________________________________________
61 AliTOFCalibTask::AliTOFCalibTask(const char *name) :
62 AliAnalysisTaskSE(name),
97 DefineOutput(1,TList::Class());
98 DefineOutput(2,TList::Class());
100 for (Int_t i=0;i<11;i++){
106 //______________________________________________________________________________
107 AliTOFCalibTask::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)
144 for (Int_t i=0;i<11;i++){
145 fassparticle[i]=calibtask.fassparticle[i];
149 //______________________________________________________________________________
150 AliTOFCalibTask:: ~AliTOFCalibTask()
154 AliInfo("TOF Calib Task: Deleting");
165 delete fhESDselected;
170 delete fhESDassTOFcl;
171 // delete fListOfHistos;
173 //______________________________________________________________________________
174 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
176 //assignment operator
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;
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;
209 for (Int_t i=0;i<11;i++){
210 fassparticle[i]=calibtask.fassparticle[i];
215 //--------------------------------------------------------------------------
216 void AliTOFCalibTask::BookHistos(){
220 AliInfo(Form("*** Booking Histograms %s", GetName())) ;
223 new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
225 new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
227 new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
229 new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
231 new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
233 new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
235 new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
238 // create the output list of histos
239 if (!fListOfHistos) fListOfHistos = new TList();
240 fListOfHistos->SetOwner();
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) ;
251 new TH1I("hESD","Number of analyzed ESDs",1,0,1);
253 new TH1I("hESDselected","Number of selected ESDs",1,0,1);
255 new TH1I("hESDkTOFout","Number of ESDs with kTOFout",1,0,1);
257 new TH1I("hESDkTIME","Number of ESDs with kTime",1,0,1);
259 new TH1I("hESDTRDcut","Number of ESDs with TRDcut",1,0,1);
261 new TH1I("hESDTIMEcut","Number of ESDs with TIMEcut",1,0,1);
263 new TH1I("hESDassTOFcl","Number of ESDs with assTOFcl",1,0,1);
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) ;
276 //----------------------------------------------------------------------------
277 void AliTOFCalibTask::DrawHistos(){
279 // drawing output histos
281 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
283 if (!gROOT->IsBatch()){
285 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
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);
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);
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");
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");
331 //canvasExpTime->Print("ExpTime.gif");
333 TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
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");
341 //canvasPID->Print("PID.gif");
343 TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
345 fhch->GetXaxis()->SetTitle("TOF ch");
347 Float_t meanTOFch = 0;
348 for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
349 meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
352 meanTOFch/=TOFCHANNELS;
353 AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
355 //canvasrndch->Print("rndch.gif");
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())) ;
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()));
377 //________________________________________________________________________
378 void AliTOFCalibTask::UserCreateOutputObjects()
380 // Create histograms and AliTOFArray
382 AliInfo("UserCreateObjects");
386 if (!fTOFArray) fTOFArray = new AliTOFArray(TOFCHANNELS);
387 for (Int_t i =0;i<TOFCHANNELS;i++){
388 fTOFArray->SetArray(i);
390 if (!fListArray) fListArray = new TList();
391 fListArray->SetOwner();
392 fListArray->Add(fTOFArray);
396 //______________________________________________________________________________
397 void AliTOFCalibTask::UserExec(Option_t * /*opt*/)
403 //******* The loop over events -----------------------------------------------
405 AliVEvent* fESD = fInputEvent ;
407 Error("UserExec","NO EVENT FOUND!");
412 fMinTime=22E3; //ns; not used
413 Int_t ntrk = fESD->GetNumberOfTracks() ;
420 AliESDtrack * t = (AliESDtrack*)fESD->GetTrack(ntrk) ;
421 //selecting only good quality tracks
422 if (!Select(t)) continue;
424 fhESDselected->Fill(0);
425 Int_t ich = Int_t(t->GetTOFCalChannel());
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[AliPID::kSPECIESC];
433 t->GetIntegratedTimes(expTime,AliPID::kSPECIESC);
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));
442 AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, currentSize/NIDX ));
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));
459 fnESDselected+=nselected;
460 PostData(1, fListOfHistos);
461 PostData(2, fListArray);
465 //_____________________________________________________________________________
466 void AliTOFCalibTask::Terminate(Option_t *)
468 // Processing when the event loop is ended
472 TH1::AddDirectory(0);
473 AliInfo("TOF Calib Task: End of events loop");
475 fListArray = dynamic_cast<TList*>(GetOutputData(2));
476 fTOFArray = dynamic_cast<AliTOFArray*>(fListArray->At(0));
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;
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));
492 AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
495 if (i%1000 == 0) AliInfo(Form("At channel %d",i));
496 if (!CombPID(fTOFArray->GetArray(i), size)) AliError("ERROR!!!!ERROR!!!");
498 bench.Stop("CombPID");
499 bench.Print("CombPID");
501 // saving data in a tree --> obsolete code; keeping for backup with new structure
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
512 tempArray->SetArray(i,NIDXSMALL*sizeChannel);
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));
523 AliInfo("Putting object for calibration in Reference data");
525 TFile *fileout = TFile::Open("outputTOFCalibration.root","RECREATE");
530 //_____________________________________________________________________________
532 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
534 // selecting good quality tracks
535 //track selection: reconstrution to TOF:
536 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
540 fhESDkTOFout->Fill(0);
541 //IsStartedTimeIntegral
542 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
547 if (t->GetStatus() & AliESDtrack::kTRDbackup) {
548 Float_t xout= t->GetOuterParam()->GetX();
549 if (xout<364.25 && xout > 300.) return 0;
552 fhESDTRDcut->Fill(0);
553 Double_t time=t->GetTOFsignal();
554 time*=1.E-3; // tof given in nanoseconds
555 if(time >= fMinTime){
559 fhESDTIMEcut->Fill(0);
561 Double_t mom=t->GetP();
562 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
563 // return 0; // skipping momentum cut
566 Int_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
567 if(assignedTOFcluster==-1){ // not matched
571 fhESDassTOFcl->Fill(0);
572 AliDebug(2,"selecting the track");
575 //_____________________________________________________________________________
577 Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
579 // to be re-implemented with new object
581 // discarding tracks with time-mintime < MINTIME
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;
590 AliDebug(1,Form("Mintime for channel %i = %f",ich, mintime));
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) {
598 AliDebug(1,Form("Deleting %i track from channel %i, time = %f",ndeleted, ich, time));
600 for (Int_t j=itr+1;j<ntracks;j++){
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];
624 //_____________________________________________________________________________
626 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
628 // track Combinatorial PID for calibration, only when more than 2 particles
634 AliDebug(1,"Number of tracks in channel smaller than 2, identifying every particle as if it was a pion!");
638 Int_t ntracksinchannel = size/NIDX;
640 Int_t nset = ntracksinchannel/6;
644 if (ntracksinchannel < 6) {
645 AliDebug(2,"Number of tracks in set smaller than 6, Combinatorial PID still applied to this set.");
649 AliDebug(2,Form("number of sets = %i", nset));
651 for (Int_t i=0; i< nset; i++) {
652 if (i<nset-1)ntrkinset=6;
654 if (ntracksinchannel<6){
655 ntrkinset=ntracksinchannel;
658 ntrkinset=6+Int_t((ntracksinchannel)%6);
661 AliDebug(2,Form("set = %i, number of tracks in set = %i", i,ntrkinset));
663 for (Int_t ii=0;ii<ntrkinset;ii++){
667 Float_t* timeofflight;
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];
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]));
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);
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);
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;
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]));
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;
725 else if (fassparticle[kk]==1){ //assigned to be a Ka
726 smallarray[idxextimePi]=smallarray[idxextimeKa];
727 smallarray[idxextimeKa]=-1;
728 smallarray[idxextimePr]=-1;
730 else if (fassparticle[kk]==2){ //assigned to be a Pr
731 smallarray[idxextimePi]=smallarray[idxextimePr];
732 smallarray[idxextimeKa]=-1;
733 smallarray[idxextimePr]=-1;
737 delete[] timeofflight;
744 //---------------------------------------------------------------------
745 Float_t AliTOFCalibTask::LoopCombPID(Int_t itrkinset, Int_t ntrkinset, Float_t **exptof, Float_t *texp, Float_t *timeofflight, Int_t *index, Float_t chisquarebest){
747 // performing combinatorial PID in recursive way
749 Int_t indextr = ntrkinset-itrkinset;
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);
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];
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++)
775 meantzero=meantzero/sumAllweights; // it is given in [ns]
777 // calculate chisquare
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++)
785 for (Int_t j=0;j<ntrkinset;j++){
786 if(index[j]==0)npion++;
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;
796 for (Int_t j=0;j<ntrkinset;j++){
797 fassparticle[j]=index[j];
800 chisquarebest=chisquare;
802 delete[] besttimezero;
803 delete[] bestchisquare;
804 delete[] bestweightedtimezero;
805 delete[] weightedtimezero;
809 return chisquarebest;