2e6d680bd8d7e5311ea2b2de3d5af28d919ed44f
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDResponse.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
17
18 //-----------------------------------------------------------------
19 //        Base class for handling the pid response               //
20 //        functions of all detectors                             //
21 //        and give access to the nsigmas                         //
22 //                                                               //
23 //   Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch //
24 //-----------------------------------------------------------------
25
26 #include <TList.h>
27 #include <TObjArray.h>
28 #include <TPRegexp.h>
29 #include <TF1.h>
30 #include <TSpline.h>
31 #include <TFile.h>
32
33 #include <AliVEvent.h>
34 #include <AliVTrack.h>
35 #include <AliLog.h>
36 #include <AliPID.h>
37 #include <AliOADBContainer.h>
38 #include <AliTRDPIDParams.h>
39 #include <AliTRDPIDReference.h>
40
41 #include "AliPIDResponse.h"
42
43 #include "AliCentrality.h"
44
45 ClassImp(AliPIDResponse);
46
47 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
48 TNamed("PIDResponse","PIDResponse"),
49 fITSResponse(isMC),
50 fTPCResponse(),
51 fTRDResponse(),
52 fTOFResponse(),
53 fEMCALResponse(),
54 fRange(5.),
55 fITSPIDmethod(kITSTruncMean),
56 fIsMC(isMC),
57 fOADBPath(),
58 fBeamType("PP"),
59 fLHCperiod(),
60 fMCperiodTPC(),
61 fMCperiodUser(),
62 fCurrentFile(),
63 fRecoPass(0),
64 fRecoPassUser(-1),
65 fRun(0),
66 fOldRun(0),
67 fArrPidResponseMaster(0x0),
68 fResolutionCorrection(0x0),
69 fTRDPIDParams(0x0),
70 fTRDPIDReference(0x0),
71 fTOFTimeZeroType(kBest_T0),
72 fTOFres(100.),
73 fTOFtail(1.1),
74 fEMCALPIDParams(0x0),
75 fCurrentEvent(0x0),
76 fCurrCentrality(0.0)
77 {
78   //
79   // default ctor
80   //
81   AliLog::SetClassDebugLevel("AliPIDResponse",10);
82   AliLog::SetClassDebugLevel("AliESDpid",10);
83   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
84
85   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
86 }
87
88 //______________________________________________________________________________
89 AliPIDResponse::~AliPIDResponse()
90 {
91   //
92   // dtor
93   //
94   delete fArrPidResponseMaster;
95   delete fTRDPIDParams;
96   delete fTRDPIDReference;
97 }
98
99 //______________________________________________________________________________
100 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
101 TNamed(other),
102 fITSResponse(other.fITSResponse),
103 fTPCResponse(other.fTPCResponse),
104 fTRDResponse(other.fTRDResponse),
105 fTOFResponse(other.fTOFResponse),
106 fEMCALResponse(other.fEMCALResponse),
107 fRange(other.fRange),
108 fITSPIDmethod(other.fITSPIDmethod),
109 fIsMC(other.fIsMC),
110 fOADBPath(other.fOADBPath),
111 fBeamType("PP"),
112 fLHCperiod(),
113 fMCperiodTPC(),
114 fMCperiodUser(other.fMCperiodUser),
115 fCurrentFile(),
116 fRecoPass(0),
117 fRecoPassUser(other.fRecoPassUser),
118 fRun(0),
119 fOldRun(0),
120 fArrPidResponseMaster(0x0),
121 fResolutionCorrection(0x0),
122 fTRDPIDParams(0x0),
123 fTRDPIDReference(0x0),
124 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
125 fTOFres(100.),
126 fTOFtail(1.1),
127 fEMCALPIDParams(0x0),
128 fCurrentEvent(0x0),
129 fCurrCentrality(0.0)
130 {
131   //
132   // copy ctor
133   //
134   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
135 }
136
137 //______________________________________________________________________________
138 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
139 {
140   //
141   // copy ctor
142   //
143   if(this!=&other) {
144     delete fArrPidResponseMaster;
145     TNamed::operator=(other);
146     fITSResponse=other.fITSResponse;
147     fTPCResponse=other.fTPCResponse;
148     fTRDResponse=other.fTRDResponse;
149     fTOFResponse=other.fTOFResponse;
150     fEMCALResponse=other.fEMCALResponse;
151     fRange=other.fRange;
152     fITSPIDmethod=other.fITSPIDmethod;
153     fOADBPath=other.fOADBPath;
154     fIsMC=other.fIsMC;
155     fBeamType="PP";
156     fLHCperiod="";
157     fMCperiodTPC="";
158     fMCperiodUser=other.fMCperiodUser;
159     fCurrentFile="";
160     fRecoPass=0;
161     fRecoPassUser=other.fRecoPassUser;
162     fRun=0;
163     fOldRun=0;
164     fArrPidResponseMaster=0x0;
165     fResolutionCorrection=0x0;
166     fTRDPIDParams=0x0;
167     fTRDPIDReference=0x0;
168     fEMCALPIDParams=0x0;
169     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
170     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
171     fTOFres=100.;
172     fTOFtail=1.1;
173     fCurrentEvent=other.fCurrentEvent;
174   }
175   return *this;
176 }
177
178 //______________________________________________________________________________
179 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
180 {
181   //
182   // NumberOfSigmas for 'detCode'
183   //
184
185   switch (detCode){
186     case kDetITS: return NumberOfSigmasITS(track, type); break;
187     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
188     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
189 //     case kDetTRD: return ComputeTRDProbability(track, type); break;
190 //     case kDetPHOS: return ComputePHOSProbability(track, type); break;
191 //     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
192 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
193     default: return -999.;
194   }
195
196 }
197
198 //______________________________________________________________________________
199 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
200
201   AliVCluster *matchedClus = NULL;
202
203   Double_t mom     = -1.; 
204   Double_t pt      = -1.; 
205   Double_t EovP    = -1.;
206   Double_t fClsE   = -1.;
207   
208   Int_t nMatchClus = -1;
209   Int_t charge     = 0;
210   
211   // Track matching
212   nMatchClus = track->GetEMCALcluster();
213   if(nMatchClus > -1){
214     
215     mom    = track->P();
216     pt     = track->Pt();
217     charge = track->Charge();
218     
219     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
220     
221     if(matchedClus){
222       
223       // matched cluster is EMCAL
224       if(matchedClus->IsEMCAL()){
225         
226         fClsE       = matchedClus->E();
227         EovP        = fClsE/mom;
228         
229         
230         // NSigma value really meaningful only for electrons!
231         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
232       }
233     }
234   }
235   
236   return -999;
237   
238 }
239
240 //______________________________________________________________________________
241 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
242
243   AliVCluster *matchedClus = NULL;
244
245   Double_t mom     = -1.; 
246   Double_t pt      = -1.; 
247   Double_t EovP    = -1.;
248   Double_t fClsE   = -1.;
249   
250   Int_t nMatchClus = -1;
251   Int_t charge     = 0;
252   
253   // Track matching
254   nMatchClus = track->GetEMCALcluster();
255   if(nMatchClus > -1){
256
257     mom    = track->P();
258     pt     = track->Pt();
259     charge = track->Charge();
260     
261     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
262     
263     if(matchedClus){
264       
265       // matched cluster is EMCAL
266       if(matchedClus->IsEMCAL()){
267         
268         fClsE       = matchedClus->E();
269         EovP        = fClsE/mom;
270         
271         // fill used EMCAL variables here
272         eop            = EovP; // E/p
273         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
274         showershape[1] = matchedClus->GetM02(); // long axis
275         showershape[2] = matchedClus->GetM20(); // short axis
276         showershape[3] = matchedClus->GetDispersion(); // dispersion
277         
278         // NSigma value really meaningful only for electrons!
279         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
280       }
281     }
282   }
283   return -999;
284 }
285
286
287 //______________________________________________________________________________
288 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
289 {
290   //
291   // Compute PID response of 'detCode'
292   //
293
294   switch (detCode){
295     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
296     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
297     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
298     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
299     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
300     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
301     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
302     default: return kDetNoSignal;
303   }
304 }
305
306 //______________________________________________________________________________
307 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
308 {
309   //
310   // Compute PID response for the ITS
311   //
312
313   // set flat distribution (no decision)
314   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
315
316   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
317     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
318   
319   Double_t mom=track->P();
320   Double_t dedx=track->GetITSsignal();
321   Bool_t isSA=kTRUE;
322   Double_t momITS=mom;
323   ULong_t trStatus=track->GetStatus();
324   if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
325   UChar_t clumap=track->GetITSClusterMap();
326   Int_t nPointsForPid=0;
327   for(Int_t i=2; i<6; i++){
328     if(clumap&(1<<i)) ++nPointsForPid;
329   }
330   
331   if(nPointsForPid<3) { // track not to be used for combined PID purposes
332     //       track->ResetStatus(AliVTrack::kITSpid);
333     return kDetNoSignal;
334   }
335
336   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
337   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
338     Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
339     Double_t bethe=fITSResponse.Bethe(momITS,mass);
340     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
341     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
342       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
343     } else {
344       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
345       mismatch=kFALSE;
346     }
347
348     // Check for particles heavier than (AliPID::kSPECIES - 1)
349     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
350
351   }
352
353   if (mismatch){
354     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
355     return kDetNoSignal;
356   }
357
358     
359   return kDetPidOk;
360 }
361 //______________________________________________________________________________
362 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
363 {
364   //
365   // Compute PID response for the TPC
366   //
367
368   // set flat distribution (no decision)
369   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
370
371   // check quality of the track
372   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
373
374   Double_t mom = track->GetTPCmomentum();
375
376   Double_t dedx=track->GetTPCsignal();
377   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
378
379   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
380     AliPID::EParticleType type=AliPID::EParticleType(j);
381     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
382     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
383     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
384       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
385     } else {
386       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
387       mismatch=kFALSE;
388     }
389
390     // TODO: Light nuclei, also in TPC pid response
391     
392     // Check for particles heavier than (AliPID::kSPECIES - 1)
393 //     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
394
395   }
396
397   if (mismatch){
398     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
399     return kDetNoSignal;
400   }
401
402   return kDetPidOk;
403 }
404 //______________________________________________________________________________
405 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
406 {
407   //
408   // Compute PID response for the
409   //
410
411   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
412
413   // set flat distribution (no decision)
414   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
415   
416   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
417   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
418   
419   Double_t time[AliPID::kSPECIESN];
420   track->GetIntegratedTimes(time);
421   
422   Double_t sigma[AliPID::kSPECIES];
423   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
424     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
425   }
426   
427   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
428   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
429     AliPID::EParticleType type=AliPID::EParticleType(j);
430     Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
431
432     Double_t sig = sigma[j];
433     if (TMath::Abs(nsigmas) > (fRange+2)) {
434       if(nsigmas < fTOFtail)
435         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
436       else
437         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
438     } else{
439       if(nsigmas < fTOFtail)
440         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
441       else
442         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
443     }
444
445     /* OLD Gaussian shape
446     if (TMath::Abs(nsigmas) > (fRange+2)) {
447       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
448     } else
449       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
450     */
451
452     if (TMath::Abs(nsigmas)<5.){
453       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
454       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
455     }
456   }
457
458   if (mismatch){
459     return kDetMismatch;    
460   }
461
462     // TODO: Light nuclei
463     
464   return kDetPidOk;
465 }
466 //______________________________________________________________________________
467 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
468 {
469   //
470   // Compute PID response for the
471   //
472
473   // set flat distribution (no decision)
474   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
475   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
476
477   Float_t mom[6];
478   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
479   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
480   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
481   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
482     mom[ilayer] = track->GetTRDmomentum(ilayer);
483     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
484       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
485     }
486   }
487   fTRDResponse.GetResponse(nslices, dedx, mom, p);
488   return kDetPidOk;
489 }
490 //______________________________________________________________________________
491 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
492 {
493   //
494   // Compute PID response for the EMCAL
495   //
496
497   AliVCluster *matchedClus = NULL;
498
499   Double_t mom     = -1.; 
500   Double_t pt      = -1.; 
501   Double_t EovP    = -1.;
502   Double_t fClsE   = -1.;
503   
504   Int_t nMatchClus = -1;
505   Int_t charge     = 0;
506   
507   // Track matching
508   nMatchClus = track->GetEMCALcluster();
509
510   if(nMatchClus > -1){
511
512     mom    = track->P();
513     pt     = track->Pt();
514     charge = track->Charge();
515     
516     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
517     
518     if(matchedClus){    
519
520     // matched cluster is EMCAL
521     if(matchedClus->IsEMCAL()){
522
523       fClsE       = matchedClus->E();
524       EovP        = fClsE/mom;
525       
526       
527       // compute the probabilities 
528       if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){     
529
530         // in case everything is OK
531         return kDetPidOk;
532         
533       }
534     }
535   }
536   }
537   
538   // in all other cases set flat distribution (no decision)
539   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
540   return kDetNoSignal;
541   
542 }
543 //______________________________________________________________________________
544 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
545 {
546   //
547   // Compute PID response for the PHOS
548   //
549
550   // set flat distribution (no decision)
551   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
552   return kDetNoSignal;
553 }
554 //______________________________________________________________________________
555 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
556 {
557   //
558   // Compute PID response for the HMPID
559   //
560
561   // set flat distribution (no decision)
562   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
563   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
564
565   track->GetHMPIDpid(p);
566   
567   return kDetPidOk;
568 }
569
570 //______________________________________________________________________________
571 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
572 {
573   //
574   // Apply settings for the current event
575   //
576   fRecoPass=pass;
577   
578   fCurrentEvent=0x0;
579   if (!event) return;
580   fCurrentEvent=event;
581   fRun=event->GetRunNumber();
582   
583   if (fRun!=fOldRun){
584     ExecNewRun();
585     fOldRun=fRun;
586   }
587   
588   //TPC resolution parametrisation PbPb
589   if ( fResolutionCorrection ){
590     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
591     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
592   }
593   
594   //TOF resolution
595   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
596
597   // Get and set centrality
598   AliCentrality *centrality = event->GetCentrality();
599   if(centrality){
600     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
601   }
602   else{
603     fCurrCentrality = -1;
604   }
605 }
606
607 //______________________________________________________________________________
608 void AliPIDResponse::ExecNewRun()
609 {
610   //
611   // Things to Execute upon a new run
612   //
613   SetRecoInfo();
614   
615   SetITSParametrisation();
616   
617   SetTPCPidResponseMaster();
618   SetTPCParametrisation();
619
620   SetTRDPidResponseMaster(); 
621   InitializeTRDResponse();
622
623   SetEMCALPidResponseMaster(); 
624   InitializeEMCALResponse();
625   
626   fTOFResponse.SetTimeResolution(fTOFres);
627 }
628
629 //_____________________________________________________
630 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
631 {
632   //
633   // Get TPC multiplicity in bins of 150
634   //
635   
636   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
637   Double_t tpcMulti=0.;
638   if(vertexTPC){
639     Double_t vertexContribTPC=vertexTPC->GetNContributors();
640     tpcMulti=vertexContribTPC/150.;
641     if (tpcMulti>20.) tpcMulti=20.;
642   }
643   
644   return tpcMulti;
645 }
646
647 //______________________________________________________________________________
648 void AliPIDResponse::SetRecoInfo()
649 {
650   //
651   // Set reconstruction information
652   //
653   
654   //reset information
655   fLHCperiod="";
656   fMCperiodTPC="";
657   
658   fBeamType="";
659     
660   fBeamType="PP";
661   
662   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
663   //find the period by run number (UGLY, but not stored in ESD and AOD... )
664   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
665   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
666   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
667   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
668   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
669   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
670   else if (fRun>=136851&&fRun<=139517) {
671     fLHCperiod="LHC10H";
672     fMCperiodTPC="LHC10H8";
673     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
674     fBeamType="PBPB";
675   }
676   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
677   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
678   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
679   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
680   else if (fRun>=166529) {
681     fLHCperiod="LHC11H";
682     fMCperiodTPC="LHC11A10";
683     fBeamType="PBPB";
684   }
685
686
687   //exception new pp MC productions from 2011
688   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
689 }
690
691 //______________________________________________________________________________
692 void AliPIDResponse::SetITSParametrisation()
693 {
694   //
695   // Set the ITS parametrisation
696   //
697 }
698
699 //______________________________________________________________________________
700 void AliPIDResponse::SetTPCPidResponseMaster()
701 {
702   //
703   // Load the TPC pid response functions from the OADB
704   //
705   //don't load twice for the moment
706    if (fArrPidResponseMaster) return;
707  
708
709   //reset the PID response functions
710   delete fArrPidResponseMaster;
711   fArrPidResponseMaster=0x0;
712   
713   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
714   
715   TFile *f=TFile::Open(fileName.Data());
716   if (f && f->IsOpen() && !f->IsZombie()){
717     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
718   }
719   delete f;
720   
721   if (!fArrPidResponseMaster){
722     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
723     return;
724   }
725   fArrPidResponseMaster->SetOwner();
726 }
727
728 //______________________________________________________________________________
729 void AliPIDResponse::SetTPCParametrisation()
730 {
731   //
732   // Change BB parametrisation for current run
733   //
734   
735   if (fLHCperiod.IsNull()) {
736     AliFatal("No period set, not changing parametrisation");
737     return;
738   }
739   
740   //
741   // Set default parametrisations for data and MC
742   //
743   
744   //data type
745   TString datatype="DATA";
746   //in case of mc fRecoPass is per default 1
747   if (fIsMC) {
748     datatype="MC";
749     fRecoPass=1;
750   }
751   
752   //
753   //reset old splines
754   //
755   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
756     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
757   }
758   
759   //
760   //set the new PID splines
761   //
762   TString period=fLHCperiod;
763   if (fArrPidResponseMaster){
764     TObject *grAll=0x0;
765     //for MC don't use period information
766 //     if (fIsMC) period="[A-Z0-9]*";
767     //for MC use MC period information
768     if (fIsMC) period=fMCperiodTPC;
769 //pattern for the default entry (valid for all particles)
770     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
771     
772     //loop over entries and filter them
773     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
774       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
775       if (responseFunction==0x0) continue;
776       TString responseName=responseFunction->GetName();
777       
778       if (!reg.MatchB(responseName)) continue;
779       
780       TObjArray *arr=reg.MatchS(responseName);
781       TString particleName=arr->At(1)->GetName();
782       delete arr;
783       if (particleName.IsNull()) continue;
784       if (particleName=="ALL") grAll=responseFunction;
785       else {
786         //find particle id
787         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
788           TString particle=AliPID::ParticleName(ispec);
789           particle.ToUpper();
790           if ( particle == particleName ){
791             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
792             fTPCResponse.SetUseDatabase(kTRUE);
793             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
794             break;
795           }
796         }
797       }
798     }
799     
800     //set default response function to all particles which don't have a specific one
801     if (grAll){
802       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
803         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
804           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
805           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
806         }
807       }
808     }
809   }
810   
811   //
812   // Setup resolution parametrisation
813   //
814   
815   //default
816   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
817   
818   if (fRun>=122195){
819     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
820   }
821   if (fArrPidResponseMaster)
822   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
823   
824   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
825 }
826
827 //______________________________________________________________________________
828 void AliPIDResponse::SetTRDPidResponseMaster()
829 {
830   //
831   // Load the TRD pid params and references from the OADB
832   //
833   if(fTRDPIDParams) return;
834   AliOADBContainer contParams("contParams"); 
835
836   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
837   if(statusPars){
838     AliError("Failed initializing PID Params from OADB");
839   } else {
840     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
841     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
842     if(!fTRDPIDParams){
843       AliError(Form("TRD Params not found in run %d", fRun));
844     }
845   }
846
847   AliOADBContainer contRefs("contRefs");
848   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
849   if(statusRefs){
850     AliInfo("Failed Loading References for TRD");
851   } else {
852     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
853     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
854     if(!fTRDPIDReference){
855       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
856     }
857   }
858 }
859
860 //______________________________________________________________________________
861 void AliPIDResponse::InitializeTRDResponse(){
862   //
863   // Set PID Params and references to the TRD PID response
864   // 
865   fTRDResponse.SetPIDParams(fTRDPIDParams);
866   fTRDResponse.Load(fTRDPIDReference);
867   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
868     fTRDslicesForPID[0] = 0;
869     fTRDslicesForPID[1] = 7;
870   }
871 }
872
873 //_________________________________________________________________________
874 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
875   //
876   // Check whether track is identified as electron under a given electron efficiency hypothesis
877   //
878   Double_t probs[AliPID::kSPECIES];
879   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
880
881   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
882   // Take mean of the TRD momenta in the given tracklets
883   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
884   Int_t nmomenta = 0;
885   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
886     if(vtrack->GetTRDmomentum(iPl) > 0.){
887       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
888     }
889   }
890   p = TMath::Mean(nmomenta, trdmomenta);
891
892   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
893 }
894
895 //______________________________________________________________________________
896 void AliPIDResponse::SetEMCALPidResponseMaster()
897 {
898   //
899   // Load the EMCAL pid response functions from the OADB
900   //
901   TObjArray* fEMCALPIDParamsRun      = NULL;
902   TObjArray* fEMCALPIDParamsPass     = NULL;
903
904   if(fEMCALPIDParams) return;
905   AliOADBContainer contParams("contParams"); 
906
907   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
908   if(statusPars){
909     AliError("Failed initializing PID Params from OADB");
910   } 
911   else {
912     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
913
914     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
915     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
916     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
917
918     if(!fEMCALPIDParams){
919       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
920       AliInfo("Will take the standard LHC11d instead ...");
921
922       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
923       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
924       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
925
926       if(!fEMCALPIDParams){
927         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
928       }
929     }
930   }
931 }
932
933 //______________________________________________________________________________
934 void AliPIDResponse::InitializeEMCALResponse(){
935   //
936   // Set PID Params to the EMCAL PID response
937   // 
938   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
939
940 }
941 //_________________________________________________________________________
942 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
943   //
944   // Set TOF response function
945   // Input option for event_time used
946   //
947   
948     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
949     if(t0spread < 10) t0spread = 80;
950
951     // T0 from TOF algorithm
952
953     Bool_t flagT0TOF=kFALSE;
954     Bool_t flagT0T0=kFALSE;
955     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
956     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
957     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
958
959     // T0-TOF arrays
960     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
961     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
962     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
963       estimatedT0event[i]=0.0;
964       estimatedT0resolution[i]=0.0;
965       startTimeMask[i] = 0;
966     }
967
968     Float_t resT0A=75,resT0C=65,resT0AC=55;
969     if(vevent->GetT0TOF()){ // check if T0 detector information is available
970         flagT0T0=kTRUE;
971     }
972
973
974     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
975
976     if (tofHeader) { // read global info and T0-TOF
977       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
978       t0spread = tofHeader->GetT0spread(); // read t0 sprad
979       if(t0spread < 10) t0spread = 80;
980
981       flagT0TOF=kTRUE;
982       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
983         startTime[i]=tofHeader->GetDefaultEventTimeVal();
984         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
985         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
986       }
987
988       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
989       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
990       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
991       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
992         Int_t icurrent = (Int_t)ibin->GetAt(j);
993         startTime[icurrent]=t0Bin->GetAt(j);
994         startTimeRes[icurrent]=t0ResBin->GetAt(j);
995         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
996       }
997     }
998
999     // for cut of 3 sigma on t0 spread
1000     Float_t t0cut = 3 * t0spread;
1001     if(t0cut < 500) t0cut = 500;
1002
1003     if(option == kFILL_T0){ // T0-FILL is used
1004         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1005           estimatedT0event[i]=0.0;
1006           estimatedT0resolution[i]=t0spread;
1007         }
1008         fTOFResponse.SetT0event(estimatedT0event);
1009         fTOFResponse.SetT0resolution(estimatedT0resolution);
1010     }
1011
1012     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1013         if(flagT0TOF){
1014             fTOFResponse.SetT0event(startTime);
1015             fTOFResponse.SetT0resolution(startTimeRes);
1016             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1017               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1018               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1019             }
1020         }
1021         else{
1022             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1023               estimatedT0event[i]=0.0;
1024               estimatedT0resolution[i]=t0spread;
1025               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1026             }
1027             fTOFResponse.SetT0event(estimatedT0event);
1028             fTOFResponse.SetT0resolution(estimatedT0resolution);
1029         }
1030     }
1031     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1032         Float_t t0AC=-10000;
1033         Float_t t0A=-10000;
1034         Float_t t0C=-10000;
1035         if(flagT0T0){
1036             t0AC= vevent->GetT0TOF()[0];
1037             t0A= vevent->GetT0TOF()[1];
1038             t0C= vevent->GetT0TOF()[2];
1039         }
1040
1041         Float_t t0t0Best = 0;
1042         Float_t t0t0BestRes = 9999;
1043         Int_t t0used=0;
1044         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1045             t0t0Best = t0AC;
1046             t0t0BestRes = resT0AC;
1047             t0used=6;
1048         }
1049         else if(TMath::Abs(t0C) < t0cut){
1050             t0t0Best = t0C;
1051             t0t0BestRes = resT0C;
1052             t0used=4;
1053         }
1054         else if(TMath::Abs(t0A) < t0cut){
1055             t0t0Best = t0A;
1056             t0t0BestRes = resT0A;
1057             t0used=2;
1058         }
1059
1060         if(flagT0TOF){ // if T0-TOF info is available
1061             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1062                 if(t0t0BestRes < 999){
1063                   if(startTimeRes[i] < t0spread){
1064                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1065                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1066                     estimatedT0event[i]=t0best / wtot;
1067                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1068                     startTimeMask[i] = t0used+1;
1069                   }
1070                   else {
1071                     estimatedT0event[i]=t0t0Best;
1072                     estimatedT0resolution[i]=t0t0BestRes;
1073                     startTimeMask[i] = t0used;
1074                   }
1075                 }
1076                 else{
1077                   estimatedT0event[i]=startTime[i];
1078                   estimatedT0resolution[i]=startTimeRes[i];
1079                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1080                 }
1081                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1082             }
1083             fTOFResponse.SetT0event(estimatedT0event);
1084             fTOFResponse.SetT0resolution(estimatedT0resolution);
1085         }
1086         else{ // if no T0-TOF info is available
1087             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1088               fTOFResponse.SetT0binMask(i,t0used);
1089               if(t0t0BestRes < 999){
1090                 estimatedT0event[i]=t0t0Best;
1091                 estimatedT0resolution[i]=t0t0BestRes;
1092               }
1093               else{
1094                 estimatedT0event[i]=0.0;
1095                 estimatedT0resolution[i]=t0spread;
1096               }
1097             }
1098             fTOFResponse.SetT0event(estimatedT0event);
1099             fTOFResponse.SetT0resolution(estimatedT0resolution);
1100         }
1101     }
1102
1103     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1104         Float_t t0AC=-10000;
1105         Float_t t0A=-10000;
1106         Float_t t0C=-10000;
1107         if(flagT0T0){
1108             t0AC= vevent->GetT0TOF()[0];
1109             t0A= vevent->GetT0TOF()[1];
1110             t0C= vevent->GetT0TOF()[2];
1111         }
1112
1113         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1114             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1115               estimatedT0event[i]=t0AC;
1116               estimatedT0resolution[i]=resT0AC;
1117               fTOFResponse.SetT0binMask(i,6);
1118             }
1119         }
1120         else if(TMath::Abs(t0C) < t0cut){
1121             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1122               estimatedT0event[i]=t0C;
1123               estimatedT0resolution[i]=resT0C;
1124               fTOFResponse.SetT0binMask(i,4);
1125             }
1126         }
1127         else if(TMath::Abs(t0A) < t0cut){
1128             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1129               estimatedT0event[i]=t0A;
1130               estimatedT0resolution[i]=resT0A;
1131               fTOFResponse.SetT0binMask(i,2);
1132             }
1133         }
1134         else{
1135             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1136               estimatedT0event[i]=0.0;
1137               estimatedT0resolution[i]=t0spread;
1138               fTOFResponse.SetT0binMask(i,0);
1139             }
1140         }
1141         fTOFResponse.SetT0event(estimatedT0event);
1142         fTOFResponse.SetT0resolution(estimatedT0resolution);
1143     }
1144     delete [] startTime;
1145     delete [] startTimeRes;
1146     delete [] startTimeMask;
1147     delete [] estimatedT0event;
1148     delete [] estimatedT0resolution;
1149 }