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
50 #include <TDirectory.h>
53 #include "AliTOFCalibTask.h"
54 #include "AliESDEvent.h"
55 #include "AliESDtrack.h"
57 #include "AliCDBManager.h"
58 #include "AliCDBMetaData.h"
60 #include "AliCDBStorage.h"
62 //______________________________________________________________________________
63 AliTOFCalibTask::AliTOFCalibTask(const char *name) :
64 AliAnalysisTask(name,""),
94 // Input slot #0 works with an Ntuple
96 DefineInput(0, TChain::Class());
97 // Output slot #0 writes into a TH1 container
98 DefineOutput(0, TObjArray::Class()) ;
99 fdir = gDirectory->GetPath();
100 fbigarray = new Float_t*[TOFCHANNELS];
101 findexarray = new Int_t[TOFCHANNELS];
103 for (Int_t i=0;i<11;i++){
107 for (Int_t i =0;i<TOFCHANNELS;i++){
108 fbigarray[i] = new Float_t[CHENTRIES];
110 for (Int_t j =0;j<CHENTRIES;j++){
116 //______________________________________________________________________________
117 AliTOFCalibTask::AliTOFCalibTask(const AliTOFCalibTask &calibtask) :
118 AliAnalysisTask("AliTOFCalibTask",""),
149 fChain=calibtask.fChain;
152 fTime=calibtask.fTime;
153 fExpTimePi=calibtask.fExpTimePi;
154 fExpTimeKa=calibtask.fExpTimeKa;
155 fExpTimePr=calibtask.fExpTimePr;
156 fMinTime=calibtask.fMinTime;
157 fnESD=calibtask.fnESD;
158 fnESDselected=calibtask.fnESDselected;
159 fnESDkTOFout=calibtask.fnESDkTOFout;
160 fnESDkTIME=calibtask.fnESDkTIME;
161 fnESDassTOFcl=calibtask.fnESDassTOFcl;
162 fnESDTIMEcut=calibtask.fnESDTIMEcut;
163 fnESDTRDcut=calibtask.fnESDTRDcut;
164 fhToT=calibtask.fhToT;
165 fhTime=calibtask.fhTime;
166 fhExpTimePi=calibtask.fhExpTimePi;
167 fhExpTimeKa=calibtask.fhExpTimeKa;
168 fhExpTimePr=calibtask.fhExpTimePr;
169 fhPID=calibtask.fhPID;
172 fOutputContainer=calibtask.fOutputContainer;
174 fbigarray = new Float_t*[TOFCHANNELS];
175 findexarray = new Int_t[TOFCHANNELS];
177 for (Int_t i=0;i<11;i++){
178 fassparticle[i]=calibtask.fassparticle[i];
181 for (Int_t i =0;i<TOFCHANNELS;i++){
182 fbigarray[i] = new Float_t[CHENTRIES];
183 findexarray[i]=calibtask.findexarray[i];
184 for (Int_t j =0;j<CHENTRIES;j++){
185 fbigarray[i][j]=calibtask.fbigarray[i][j];
190 //______________________________________________________________________________
191 AliTOFCalibTask:: ~AliTOFCalibTask()
195 AliInfo("TOF Calib Task: Deleting");
197 delete[] findexarray;
198 delete fOutputContainer;
207 //______________________________________________________________________________
208 AliTOFCalibTask& AliTOFCalibTask::operator=(const AliTOFCalibTask &calibtask)
210 //assignment operator
211 this->fdir=calibtask.fdir;
212 this->fChain=calibtask.fChain;
213 this->fESD=calibtask.fESD;
214 this->fToT=calibtask.fToT;
215 this->fTime=calibtask.fTime;
216 this->fExpTimePi=calibtask.fExpTimePi;
217 this->fExpTimeKa=calibtask.fExpTimeKa;
218 this->fExpTimePr=calibtask.fExpTimePr;
219 this->fMinTime=calibtask.fMinTime;
220 this->fnESD=calibtask.fnESD;
221 this->fnESDselected=calibtask.fnESDselected;
222 this->fnESDkTOFout=calibtask.fnESDkTOFout;
223 this->fnESDkTIME=calibtask.fnESDkTIME;
224 this->fnESDassTOFcl=calibtask.fnESDassTOFcl;
225 this->fnESDTIMEcut=calibtask.fnESDTIMEcut;
226 this->fnESDTRDcut=calibtask.fnESDTRDcut;
227 this->fhToT=calibtask.fhToT;
228 this->fhTime=calibtask.fhTime;
229 this->fhExpTimePi=calibtask.fhExpTimePi;
230 this->fhExpTimeKa=calibtask.fhExpTimeKa;
231 this->fhExpTimePr=calibtask.fhExpTimePr;
232 this->fOutputContainer=calibtask.fOutputContainer;
233 this->fhPID=calibtask.fhPID;
234 this->fhch=calibtask.fhch;
236 this->fbigarray = new Float_t*[TOFCHANNELS];
237 this->findexarray = new Int_t[TOFCHANNELS];
239 for (Int_t i=0;i<11;i++){
240 this->fassparticle[i]=calibtask.fassparticle[i];
242 for (Int_t i =0;i<TOFCHANNELS;i++){
243 this->fbigarray[i] = new Float_t[CHENTRIES];
244 this->findexarray[i]=calibtask.findexarray[i];
245 for (Int_t j =0;j<CHENTRIES;j++){
246 this->fbigarray[i][j]=calibtask.fbigarray[i][j];
251 //--------------------------------------------------------------------------
252 void AliTOFCalibTask::BookHistos(){
256 AliInfo(Form("*** Booking Histograms %s", GetName())) ;
259 new TH1F("hToT", " ToT distribution (ns) ", 400, 0, 40);
261 new TH1F("hTime", " Time distribution (ns) ", 400, 0, 40);
263 new TH1F("hExpTimePi", " Exp Time distribution, Pi (ns) ", 400, 0, 40);
265 new TH1F("hExpTimeKa", " Exp Time distribution, Ka (ns) ", 400, 0, 40);
267 new TH1F("hExpTimePr", " Exp Time distribution, Pr (ns) ", 400, 0, 40);
269 new TH1I("hPID", " Combinatorial PID identity ", 3, 0, 3);
271 new TH1D("hch", " TOF channel ", TOFCHANNELS, 0, TOFCHANNELS);
273 //create the putput container
274 fOutputContainer = new TObjArray(7) ;
275 fOutputContainer->SetName(GetName()) ;
277 fOutputContainer->AddAt(fhToT, 0) ;
278 fOutputContainer->AddAt(fhTime, 1) ;
279 fOutputContainer->AddAt(fhExpTimePi, 2) ;
280 fOutputContainer->AddAt(fhExpTimeKa, 3) ;
281 fOutputContainer->AddAt(fhExpTimePr, 4) ;
282 fOutputContainer->AddAt(fhPID, 5) ;
283 fOutputContainer->AddAt(fhch, 6) ;
286 //----------------------------------------------------------------------------
287 void AliTOFCalibTask::DrawHistos(){
289 // drawing output histos
291 AliInfo(Form("*** Drawing Histograms %s", GetName())) ;
293 TCanvas * canvasToTTime = new TCanvas("canvasToTTime", " ToT and Time ",400, 30, 550, 630) ;
294 canvasToTTime->Divide(1,2);
295 canvasToTTime->cd(1);
296 fhToT->SetLineColor(4);
297 fhToT->GetXaxis()->SetTitle("ToT (ns)");
299 canvasToTTime->cd(2);
300 fhTime->SetLineColor(4);
301 fhTime->GetXaxis()->SetTitle("Time (ns)");
302 fhTime->Draw("hist");
303 canvasToTTime->Update();
304 canvasToTTime->Print("ToTTime.gif");
306 TCanvas * canvasExpTime = new TCanvas("canvasExpTime", " Expected Times ",400, 30, 550, 630) ;
307 canvasExpTime->Divide(1,3);
308 canvasExpTime->cd(1);
309 fhExpTimePi->SetLineColor(4);
310 fhExpTimePi->GetXaxis()->SetTitle("Exp Time (ns), #pi");
311 fhExpTimePi->Draw("hist");
312 canvasExpTime->cd(2);
313 fhExpTimeKa->SetLineColor(4);
314 fhExpTimeKa->GetXaxis()->SetTitle("Exp Time (ns), K");
315 fhExpTimeKa->Draw("hist");
316 canvasExpTime->cd(3);
317 fhExpTimePr->SetLineColor(4);
318 fhExpTimePr->GetXaxis()->SetTitle("Exp Time (ns), p");
319 fhExpTimePr->Draw("hist");
321 canvasExpTime->Print("ExpTime.gif");
323 TCanvas * canvasPID = new TCanvas("canvasPID", " Combinatorial PID ",400, 30, 550, 400);
324 fhPID->GetXaxis()->SetTitle("Comb PID");
325 fhPID->GetXaxis()->SetBinLabel(1,"#pi");
326 fhPID->GetXaxis()->SetBinLabel(2,"K");
327 fhPID->GetXaxis()->SetBinLabel(3,"p");
330 canvasPID->Print("PID.gif");
332 TCanvas * canvasrndch = new TCanvas("canvasrndch", " TOF channel ",400, 30, 550, 400);
333 fhch->GetXaxis()->SetTitle("TOF ch");
335 Float_t meanTOFch = 0;
336 for (Int_t ibin=0;ibin<TOFCHANNELS;ibin++){
337 meanTOFch+=(Float_t)fhch->GetBinContent(ibin+1);
340 meanTOFch/=TOFCHANNELS;
341 AliDebug(1,Form(" Mean number of tracks/channel = %f ",meanTOFch));
343 canvasrndch->Print("rndch.gif");
346 sprintf(line, ".!tar -zcvf %s.tar.gz *.gif", GetName()) ;
347 gROOT->ProcessLine(line);
348 sprintf(line, ".!rm -fR *.gif");
349 gROOT->ProcessLine(line);
350 AliInfo(Form("*** TOF Calib Task: plots saved in %s.tar.gz...\n", GetName())) ;
353 //______________________________________________________________________________
354 void AliTOFCalibTask::ConnectInputData(const Option_t*)
356 // Initialization of branch container and histograms
358 // AliLog::SetClassDebugLevel("AliTOFCalibTask",1);
359 AliInfo(Form("*** Initialization of %s", GetName())) ;
362 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
364 AliError(Form("Input 0 for %s not found\n", GetName()));
368 // One should first check if the branch address was taken by some other task
369 char ** address = (char **)GetBranchAddress(0, "ESD");
371 fESD = (AliESDEvent*)(*address);
373 fESD = new AliESDEvent();
375 fESD->ReadFromTree(fChain) ;
380 //-----------------------------------------------------------------------
381 Bool_t AliTOFCalibTask::Notify()
383 // Initialisation of branch container and histograms
385 AliInfo(Form("*** We are in %s::Notify()", GetName())) ;
388 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
390 AliError(Form("Input 0 for %s not found\n", GetName()));
394 char ** address = (char **)GetBranchAddress(0, "ESD");
396 fESD = (AliESDEvent*)(*address);
398 fESD = new AliESDEvent();
400 fESD->ReadFromTree(fChain) ;
405 //________________________________________________________________________
406 void AliTOFCalibTask::CreateOutputObjects()
411 //______________________________________________________________________________
412 void AliTOFCalibTask::Exec(Option_t * opt)
417 AliInfo(Form("*** Executing %s", GetName())) ;
419 //******* The loop over events -----------------------------------------------
421 // Processing of one event
422 Long64_t entry = fChain->GetReadEntry() ;
424 AliError("fESD is not connected to the input!") ;
428 if ( !((entry-1)%100) )
429 AliDebug(1,Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
431 // ************************ TOF *************************************
433 fMinTime=22E3; //ns; not used
434 Int_t ntrk = fESD->GetNumberOfTracks() ;
440 AliESDtrack * t = fESD->GetTrack(ntrk) ;
441 //selecting only good quality tracks
442 if (!Select(t)) continue;
444 Int_t ich = Int_t(t->GetTOFCalChannel());
446 // ich=3; //only for debug purpose
447 AliDebug(2,Form(" ESD in channel %i, filling array %i",t->GetTOFCalChannel(),ich));
448 Float_t tot = t->GetTOFsignalToT();
449 Float_t time = t->GetTOFsignalRaw();
450 AliDebug(2,Form(" track # %i in channel %i, time = %f \n",ntrk,ich,time));
451 Double_t expTime[10];
452 t->GetIntegratedTimes(expTime);
453 Float_t expTimePi = expTime[2]*1.E-3;
454 Float_t expTimeKa = expTime[3]*1.E-3;
455 Float_t expTimePr = expTime[4]*1.E-3;
456 if (findexarray[ich]==(Int_t)(CHENTRIES/NIDX)) {
457 AliInfo(Form("too many tracks in channel %i, not storing any more...",ich));
461 AliDebug(2,Form("tracks in channel %i = %i, storing... ",ich, findexarray[ich] ));
462 Int_t ientry=(findexarray[ich]-1)*NIDX;
463 fbigarray[ich][ientry+DELTAIDXTOT]=tot; //in ns
464 fbigarray[ich][ientry+DELTAIDXTIME]=time*1E-3; // in ns
465 fbigarray[ich][ientry+DELTAIDXEXTIMEPI]=expTimePi;
466 fbigarray[ich][ientry+DELTAIDXEXTIMEKA]=expTimeKa;
467 fbigarray[ich][ientry+DELTAIDXEXTIMEPR]=expTimePr;
468 fhToT->Fill(fbigarray[ich][ientry+DELTAIDXTOT]);
469 fhTime->Fill(fbigarray[ich][ientry+DELTAIDXTIME]);
470 fhExpTimePi->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPI]);
471 fhExpTimeKa->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEKA]);
472 fhExpTimePr->Fill(fbigarray[ich][ientry+DELTAIDXEXTIMEPR]);
473 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));
475 fnESDselected+=nselected;
477 PostData(0, fOutputContainer);
480 //_____________________________________________________________________________
481 void AliTOFCalibTask::Terminate(Option_t *)
483 // Processing when the event loop is ended
487 TH1::AddDirectory(0);
488 TDirectory *dir = gDirectory;
489 AliInfo("TOF Calib Task: End of events loop");
490 for (Int_t ich = 0;ich<TOFCHANNELS;ich++){
491 if (findexarray[ich]>0){
492 Int_t ncuttime = SelectOnTime(&fbigarray[ich][0],findexarray[ich],ich);
493 fnESDselected-=ncuttime;
497 AliInfo(Form(" Number of analyzed ESD tracks: %i\n",fnESD));
498 AliInfo(Form(" Number of selected ESD tracks: %i\n",fnESDselected));
499 AliInfo(Form(" Number of ESD tracks with kTOFout: %i\n",fnESDkTOFout));
500 AliInfo(Form(" Number of ESD tracks with kTIME: %i\n",fnESDkTIME));
501 AliInfo(Form(" Number of ESD tracks with TRDcut: %i\n",fnESDTRDcut));
502 AliInfo(Form(" Number of ESD tracks with TIMEcut: %i\n",fnESDTIMEcut));
503 AliInfo(Form(" Number of ESD tracks with assTOFcl: %i\n",fnESDassTOFcl));
505 for (Int_t i = 0;i<TOFCHANNELS;i++){
506 Int_t size=findexarray[i]*NIDX;
507 AliDebug(2, Form(" entries %i in channel %i ",findexarray[i],i));
508 if (findexarray[i]<=2) {
509 AliDebug(1, Form(" not enough statistics for combined PID for channel %i, putting all the tracks as if they were pions",i));
512 if (!CombPID(&fbigarray[i][0], size)) AliError("ERROR!!!!ERROR!!!");
518 // saving data in a tree
520 AliInfo("Building tree for Calibration");
521 TTree * tree = new TTree("T", "Tree for TOF calibration");
522 Float_t p[CHENTRIESSMALL];
524 tree->Branch("nentries",&nentries,"nentries/I");
525 tree->Branch("TOFentries",p,"TOFentries[nentries]/F");
526 for (Int_t i=0;i<TOFCHANNELS;i++){
527 nentries=findexarray[i]*(NIDXSMALL); // when filling small array,
528 // only first 3 floats taken
530 for (Int_t j=0; j<findexarray[i];j++){
531 for (Int_t k=0; k<NIDXSMALL;k++){
532 Int_t index1= j*NIDXSMALL+k; // index in small array
533 Int_t index2=j*NIDX+k; // index in big array
534 p[index1]=fbigarray[i][index2];
540 AliInfo("Putting tree for calibration in Reference data");
543 Char_t filename[100];
544 sprintf(filename,"alien:///alice/cern.ch/user/c/czampolli/TOFCalibReference_%i.root",frun);
545 TGrid::Connect("alien://");
546 TFile *filegrid = TFile::Open(filename,"CREATE");
551 //_____________________________________________________________________________
553 Bool_t AliTOFCalibTask::Select(AliESDtrack *t){
555 // selecting good quality tracks
556 //track selection: reconstrution to TOF:
557 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) {
561 //IsStartedTimeIntegral
562 if ((t->GetStatus()&AliESDtrack::kTIME)==0) {
566 if (t->GetStatus() & AliESDtrack::kTRDbackup) {
567 Float_t xout= t->GetOuterParam()->GetX();
568 if (xout<364.25 && xout > 300.) return 0;
571 Double_t time=t->GetTOFsignal();
572 time*=1.E-3; // tof given in nanoseconds
573 if(time >= fMinTime){
578 Double_t mom=t->GetP();
579 if (!(mom<=UPPERMOMBOUND && mom>=LOWERMOMBOUND)){
580 // return 0; // skipping momentum cut
583 UInt_t assignedTOFcluster=t->GetTOFcluster();//index of the assigned TOF cluster, >0 ?
584 if(assignedTOFcluster==0){ // not matched
590 //_____________________________________________________________________________
592 Int_t AliTOFCalibTask::SelectOnTime(Float_t *charray, Int_t ntracks, Int_t ich){
594 // discarding tracks with time-mintime < MINTIME
597 Float_t mintime = 1E6;
598 for (Int_t itr = 0;itr<ntracks;itr++){
599 Int_t ientry=itr*NIDX;
600 Float_t time = charray[ientry+DELTAIDXTIME];// in ns
601 if (time<mintime) mintime = time;
603 AliDebug(1,Form("Mintime for channel %i = %f",ich, mintime));
604 for (Int_t itr = 0;itr<ntracks;itr++){
605 Int_t ientry=itr*NIDX;
606 Float_t time = charray[ientry+DELTAIDXTIME];// in ns
607 if ((time-mintime)>MINTIME) {
609 AliDebug(1,Form("Deleting %i track from channel %i, time = %f",ndeleted, ich, time));
611 for (Int_t j=itr+1;j<ntracks;j++){
613 Int_t idxtot = ientry+DELTAIDXTOT;
614 Int_t idxtime = ientry+DELTAIDXTIME;
615 Int_t idxextimePi = ientry+DELTAIDXEXTIMEPI;
616 Int_t idxextimeKa = ientry+DELTAIDXEXTIMEKA;
617 Int_t idxextimePr = ientry+DELTAIDXEXTIMEPR;
618 Int_t ientrydel=(j-1)*NIDX;
619 Int_t idxtotdel = ientrydel+DELTAIDXTOT;
620 Int_t idxtimedel = ientrydel+DELTAIDXTIME;
621 Int_t idxextimePidel = ientrydel+DELTAIDXEXTIMEPI;
622 Int_t idxextimeKadel = ientrydel+DELTAIDXEXTIMEKA;
623 Int_t idxextimePrdel = ientrydel+DELTAIDXEXTIMEPR;
624 charray[idxtotdel]=charray[idxtot];
625 charray[idxtimedel]=charray[idxtime];
626 charray[idxextimePidel]=charray[idxextimePi];
627 charray[idxextimeKadel]=charray[idxextimeKa];
628 charray[idxextimePrdel]=charray[idxextimePr];
634 //_____________________________________________________________________________
636 Bool_t AliTOFCalibTask::CombPID(Float_t *smallarray, Int_t size){
638 // track Combinatorial PID for calibration, only when more than 2 particles
644 AliDebug(1,"Number of tracks in channel smaller than 2, identifying every particle as if it was a pion!");
648 Int_t ntracksinchannel = size/NIDX;
650 Int_t nset = ntracksinchannel/6;
654 if (ntracksinchannel < 6) {
655 AliDebug(2,"Number of tracks in set smaller than 6, Combinatorial PID still applied to this set.");
659 AliInfo(Form("number of sets = %i", nset));
661 for (Int_t i=0; i< nset; i++) {
662 if (i<nset-1)ntrkinset=6;
664 if (ntracksinchannel<6){
665 ntrkinset=ntracksinchannel;
668 ntrkinset=6+Int_t((ntracksinchannel)%6);
671 AliInfo(Form("set = %i, number of tracks in set = %i", i,ntrkinset));
673 for (Int_t ii=0;ii<ntrkinset;ii++){
677 Float_t* timeofflight;
681 exptof = new Float_t*[ntrkinset];
682 timeofflight = new Float_t[ntrkinset];
683 texp = new Float_t[ntrkinset];
684 index = new Int_t[ntrkinset];
685 for (Int_t ii=0;ii<ntrkinset;ii++){
686 exptof[ii] = new Float_t[3];
692 for (Int_t j=0; j<ntrkinset; j++) {
693 Int_t idxtime = ((6*i+j)*NIDX)+DELTAIDXTIME;
694 Int_t idxextimePi = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPI;
695 Int_t idxextimeKa = ((6*i+j)*NIDX)+DELTAIDXEXTIMEKA;
696 Int_t idxextimePr = ((6*i+j)*NIDX)+DELTAIDXEXTIMEPR;
697 AliDebug(2,Form("idxtime = %i, idxextimePi = %i, idxextimeKa = %i, idxextimePr = %i", idxtime, idxextimePi, idxextimeKa, idxextimePr));
698 Double_t time=smallarray[idxtime]; // TOF time in ns
699 timeofflight[j]=time+t0offset;
700 exptof[j][0]=smallarray[idxextimePi];
701 exptof[j][1]=smallarray[idxextimeKa];
702 exptof[j][2]=smallarray[idxextimePr];
703 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]));
707 Float_t chisquarebest=999.;
708 AliInfo(Form(" Set = %i with %i tracks ", i,ntrkinset));
709 chisquarebest = LoopCombPID(ntrkinset, ntrkinset,exptof,&texp[0],&timeofflight[0], &index[0],chisquarebest);
711 Float_t confLevel=999;
712 if(chisquarebest<999.){
713 Double_t dblechisquare=(Double_t)chisquarebest;
714 confLevel=(Float_t)TMath::Prob(dblechisquare,ntrkinset-1);
717 // assume they are all pions for fake sets
718 if(confLevel<0.01 || confLevel==999. ){
719 for (Int_t itrk=0; itrk<ntrkinset; itrk++)fassparticle[itrk]=0;
722 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]));
724 for (Int_t kk=0;kk<ntrkinset;kk++){
725 Int_t idxextimePi = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPI;
726 Int_t idxextimeKa = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEKA;
727 Int_t idxextimePr = ((6*i+kk)*NIDX)+DELTAIDXEXTIMEPR;
728 // storing in third slot of smallarray the assigned expected time
729 fhPID->Fill(fassparticle[kk]);
730 // fassparticle[kk]=0; //assuming all particles are pions
731 if (fassparticle[kk]==0){ //assigned to be a Pi
732 smallarray[idxextimeKa]=-1;
733 smallarray[idxextimePr]=-1;
735 else if (fassparticle[kk]==1){ //assigned to be a Ka
736 smallarray[idxextimePi]=smallarray[idxextimeKa];
737 smallarray[idxextimeKa]=-1;
738 smallarray[idxextimePr]=-1;
740 else if (fassparticle[kk]==2){ //assigned to be a Pr
741 smallarray[idxextimePi]=smallarray[idxextimePr];
742 smallarray[idxextimeKa]=-1;
743 smallarray[idxextimePr]=-1;
747 delete[] timeofflight;
754 //---------------------------------------------------------------------
755 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){
757 Int_t indextr = ntrkinset-itrkinset;
759 for (index[indextr]=0;index[indextr]<3;index[indextr]++){
760 Int_t ii = index[indextr];
761 texp[indextr]=exptof[indextr][ii];
762 if (indextr<ntrkinset-1){
763 chisquarebest = LoopCombPID(ntrkinset-indextr-1,ntrkinset,exptof,&texp[0],&timeofflight[0],&index[0], chisquarebest);
767 Float_t *besttimezero = new Float_t[ntrkinset];
768 Float_t *bestchisquare = new Float_t[ntrkinset];
769 Float_t *bestweightedtimezero = new Float_t[ntrkinset];
770 Float_t sumAllweights=0.;
771 Float_t meantzero=0.;
772 Float_t *weightedtimezero = new Float_t[ntrkinset];
773 Float_t *timezero = new Float_t[ntrkinset];
775 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));
776 for (Int_t itz=0; itz<ntrkinset;itz++) {
777 timezero[itz]=texp[itz]-timeofflight[itz];
778 weightedtimezero[itz]=timezero[itz]/TRACKERROR;
779 sumAllweights+=1./TRACKERROR;
780 meantzero+=weightedtimezero[itz];
781 } // end loop for (Int_t itz=0; itz<ntrkinset;itz++)
783 meantzero=meantzero/sumAllweights; // it is given in [ns]
785 // calculate chisquare
787 Float_t chisquare=0.;
788 for (Int_t icsq=0; icsq<ntrkinset;icsq++) {
789 chisquare+=(timezero[icsq]-meantzero)*(timezero[icsq]-meantzero)/TRACKERROR;
790 } // end loop for (Int_t icsq=0; icsq<ntrkinset;icsq++)
793 for (Int_t j=0;j<ntrkinset;j++){
794 if(index[j]==0)npion++;
797 if(chisquare<=chisquarebest && ((Float_t) npion/ ((Float_t) ntrkinset)>0.3)){
798 for(Int_t iqsq = 0; iqsq<ntrkinset; iqsq++) {
799 besttimezero[iqsq]=timezero[iqsq];
800 bestweightedtimezero[iqsq]=weightedtimezero[iqsq];
801 bestchisquare[iqsq]=(timezero[iqsq]-meantzero)*(timezero[iqsq]-meantzero)/TRACKERROR;
804 for (Int_t j=0;j<ntrkinset;j++){
805 fassparticle[j]=index[j];
808 chisquarebest=chisquare;
810 delete[] besttimezero;
811 delete[] bestchisquare;
812 delete[] bestweightedtimezero;
813 delete[] weightedtimezero;
817 return chisquarebest;