EMCAL Response in OADB
[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 ClassImp(AliPIDResponse);
44
45 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
46 TNamed("PIDResponse","PIDResponse"),
47 fITSResponse(isMC),
48 fTPCResponse(),
49 fTRDResponse(),
50 fTOFResponse(),
51 fEMCALResponse(),
52 fRange(5.),
53 fITSPIDmethod(kITSTruncMean),
54 fIsMC(isMC),
55 fOADBPath(),
56 fBeamType("PP"),
57 fLHCperiod(),
58 fMCperiodTPC(),
59 fMCperiodUser(),
60 fCurrentFile(),
61 fRecoPass(0),
62 fRecoPassUser(-1),
63 fRun(0),
64 fOldRun(0),
65 fArrPidResponseMaster(0x0),
66 fResolutionCorrection(0x0),
67 fTRDPIDParams(0x0),
68 fTRDPIDReference(0x0),
69 fTOFTimeZeroType(kBest_T0),
70 fTOFres(100.),
71 fEMCALPIDParams(0x0),
72 fCurrentEvent(0x0)
73 {
74   //
75   // default ctor
76   //
77   AliLog::SetClassDebugLevel("AliPIDResponse",10);
78   AliLog::SetClassDebugLevel("AliESDpid",10);
79   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
80
81   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
82 }
83
84 //______________________________________________________________________________
85 AliPIDResponse::~AliPIDResponse()
86 {
87   //
88   // dtor
89   //
90   delete fArrPidResponseMaster;
91   delete fTRDPIDParams;
92   delete fTRDPIDReference;
93 }
94
95 //______________________________________________________________________________
96 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
97 TNamed(other),
98 fITSResponse(other.fITSResponse),
99 fTPCResponse(other.fTPCResponse),
100 fTRDResponse(other.fTRDResponse),
101 fTOFResponse(other.fTOFResponse),
102 fEMCALResponse(other.fEMCALResponse),
103 fRange(other.fRange),
104 fITSPIDmethod(other.fITSPIDmethod),
105 fIsMC(other.fIsMC),
106 fOADBPath(other.fOADBPath),
107 fBeamType("PP"),
108 fLHCperiod(),
109 fMCperiodTPC(),
110 fMCperiodUser(other.fMCperiodUser),
111 fCurrentFile(),
112 fRecoPass(0),
113 fRecoPassUser(other.fRecoPassUser),
114 fRun(0),
115 fOldRun(0),
116 fArrPidResponseMaster(0x0),
117 fResolutionCorrection(0x0),
118 fTRDPIDParams(0x0),
119 fTRDPIDReference(0x0),
120 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
121 fTOFres(100.),
122 fEMCALPIDParams(0x0),
123 fCurrentEvent(0x0)
124 {
125   //
126   // copy ctor
127   //
128   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
129 }
130
131 //______________________________________________________________________________
132 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
133 {
134   //
135   // copy ctor
136   //
137   if(this!=&other) {
138     delete fArrPidResponseMaster;
139     TNamed::operator=(other);
140     fITSResponse=other.fITSResponse;
141     fTPCResponse=other.fTPCResponse;
142     fTRDResponse=other.fTRDResponse;
143     fTOFResponse=other.fTOFResponse;
144     fEMCALResponse=other.fEMCALResponse;
145     fRange=other.fRange;
146     fITSPIDmethod=other.fITSPIDmethod;
147     fOADBPath=other.fOADBPath;
148     fIsMC=other.fIsMC;
149     fBeamType="PP";
150     fLHCperiod="";
151     fMCperiodTPC="";
152     fMCperiodUser=other.fMCperiodUser;
153     fCurrentFile="";
154     fRecoPass=0;
155     fRecoPassUser=other.fRecoPassUser;
156     fRun=0;
157     fOldRun=0;
158     fArrPidResponseMaster=0x0;
159     fResolutionCorrection=0x0;
160     fTRDPIDParams=0x0;
161     fTRDPIDReference=0x0;
162     fEMCALPIDParams=0x0;
163     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
164     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
165     fTOFres=100.;
166     fCurrentEvent=other.fCurrentEvent;
167   }
168   return *this;
169 }
170
171 //______________________________________________________________________________
172 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
173 {
174   //
175   // NumberOfSigmas for 'detCode'
176   //
177
178   switch (detCode){
179     case kDetITS: return NumberOfSigmasITS(track, type); break;
180     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
181     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
182 //     case kDetTRD: return ComputeTRDProbability(track, type); break;
183 //     case kDetPHOS: return ComputePHOSProbability(track, type); break;
184 //     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
185 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
186     default: return -999.;
187   }
188
189 }
190
191 //______________________________________________________________________________
192 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
193
194   AliVCluster *matchedClus = NULL;
195
196   Double_t mom     = -1.; 
197   Double_t pt      = -1.; 
198   Double_t EovP    = -1.;
199   Double_t fClsE   = -1.;
200   
201   Int_t nMatchClus = -1;
202   Int_t charge     = 0;
203   
204   // Track matching
205   nMatchClus = track->GetEMCALcluster();
206   if(nMatchClus > -1){
207
208     mom    = track->P();
209     pt     = track->Pt();
210     charge = track->Charge();
211     
212     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
213     
214     if(matchedClus){
215       
216     // matched cluster is EMCAL
217     if(matchedClus->IsEMCAL()){
218       
219       fClsE       = matchedClus->E();
220       EovP        = fClsE/mom;
221       
222       
223       // NSigma value really meaningful only for electrons!
224       return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
225     }
226   }
227   }
228
229   return -999;
230   
231 }
232
233 //______________________________________________________________________________
234 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
235 {
236   //
237   // Compute PID response of 'detCode'
238   //
239
240   switch (detCode){
241     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
242     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
243     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
244     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
245     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
246     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
247     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
248     default: return kDetNoSignal;
249   }
250 }
251
252 //______________________________________________________________________________
253 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
254 {
255   //
256   // Compute PID response for the ITS
257   //
258
259   // set flat distribution (no decision)
260   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
261
262   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
263     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
264   
265   Double_t mom=track->P();
266   Double_t dedx=track->GetITSsignal();
267   Bool_t isSA=kTRUE;
268   Double_t momITS=mom;
269   ULong_t trStatus=track->GetStatus();
270   if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
271   UChar_t clumap=track->GetITSClusterMap();
272   Int_t nPointsForPid=0;
273   for(Int_t i=2; i<6; i++){
274     if(clumap&(1<<i)) ++nPointsForPid;
275   }
276   
277   if(nPointsForPid<3) { // track not to be used for combined PID purposes
278     //       track->ResetStatus(AliVTrack::kITSpid);
279     return kDetNoSignal;
280   }
281
282   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
283   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
284     Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
285     Double_t bethe=fITSResponse.Bethe(momITS,mass);
286     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
287     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
288       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
289     } else {
290       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
291       mismatch=kFALSE;
292     }
293
294     // Check for particles heavier than (AliPID::kSPECIES - 1)
295     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
296
297   }
298
299   if (mismatch){
300     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
301     return kDetNoSignal;
302   }
303
304     
305   return kDetPidOk;
306 }
307 //______________________________________________________________________________
308 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
309 {
310   //
311   // Compute PID response for the TPC
312   //
313
314   // set flat distribution (no decision)
315   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
316
317   // check quality of the track
318   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
319
320   Double_t mom = track->GetTPCmomentum();
321
322   Double_t dedx=track->GetTPCsignal();
323   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
324
325   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
326     AliPID::EParticleType type=AliPID::EParticleType(j);
327     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
328     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
329     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
330       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
331     } else {
332       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
333       mismatch=kFALSE;
334     }
335
336     // TODO: Light nuclei, also in TPC pid response
337     
338     // Check for particles heavier than (AliPID::kSPECIES - 1)
339 //     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
340
341   }
342
343   if (mismatch){
344     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
345     return kDetNoSignal;
346   }
347
348   return kDetPidOk;
349 }
350 //______________________________________________________________________________
351 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
352 {
353   //
354   // Compute PID response for the
355   //
356
357   // set flat distribution (no decision)
358   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
359   
360   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
361   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
362   
363   Double_t time[AliPID::kSPECIESN];
364   track->GetIntegratedTimes(time);
365   
366   Double_t sigma[AliPID::kSPECIES];
367   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
368     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
369   }
370   
371   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
372   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
373     AliPID::EParticleType type=AliPID::EParticleType(j);
374     Double_t nsigmas=NumberOfSigmasTOF(track,type);
375     
376     Double_t sig = sigma[j];
377     if (TMath::Abs(nsigmas) > (fRange+2)) {
378       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
379     } else
380       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
381
382     if (TMath::Abs(nsigmas)<5.){
383       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
384       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
385     }
386   }
387
388   if (mismatch){
389     return kDetMismatch;    
390   }
391
392     // TODO: Light nuclei
393     
394   return kDetPidOk;
395 }
396 //______________________________________________________________________________
397 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
398 {
399   //
400   // Compute PID response for the
401   //
402
403   // set flat distribution (no decision)
404   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
405   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
406
407   Float_t mom[6];
408   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
409   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
410   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
411   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
412     mom[ilayer] = track->GetTRDmomentum(ilayer);
413     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
414       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
415     }
416   }
417   fTRDResponse.GetResponse(nslices, dedx, mom, p);
418   return kDetPidOk;
419 }
420 //______________________________________________________________________________
421 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
422 {
423   //
424   // Compute PID response for the EMCAL
425   //
426
427   AliVCluster *matchedClus = NULL;
428
429   Double_t mom     = -1.; 
430   Double_t pt      = -1.; 
431   Double_t EovP    = -1.;
432   Double_t fClsE   = -1.;
433   
434   Int_t nMatchClus = -1;
435   Int_t charge     = 0;
436   
437   // Track matching
438   nMatchClus = track->GetEMCALcluster();
439
440   if(nMatchClus > -1){
441
442     mom    = track->P();
443     pt     = track->Pt();
444     charge = track->Charge();
445     
446     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
447     
448     if(matchedClus){    
449
450     // matched cluster is EMCAL
451     if(matchedClus->IsEMCAL()){
452
453       fClsE       = matchedClus->E();
454       EovP        = fClsE/mom;
455       
456       
457       // compute the probabilities 
458       if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){     
459
460         // in case everything is OK
461         return kDetPidOk;
462         
463       }
464     }
465   }
466   }
467   
468   // in all other cases set flat distribution (no decision)
469   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
470   return kDetNoSignal;
471   
472 }
473 //______________________________________________________________________________
474 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
475 {
476   //
477   // Compute PID response for the PHOS
478   //
479
480   // set flat distribution (no decision)
481   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
482   return kDetNoSignal;
483 }
484 //______________________________________________________________________________
485 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
486 {
487   //
488   // Compute PID response for the HMPID
489   //
490
491   // set flat distribution (no decision)
492   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
493   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
494
495   track->GetHMPIDpid(p);
496   
497   return kDetPidOk;
498 }
499
500 //______________________________________________________________________________
501 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
502 {
503   //
504   // Apply settings for the current event
505   //
506   fRecoPass=pass;
507   
508   fCurrentEvent=0x0;
509   if (!event) return;
510   fCurrentEvent=event;
511   fRun=event->GetRunNumber();
512   
513   if (fRun!=fOldRun){
514     ExecNewRun();
515     fOldRun=fRun;
516   }
517   
518   //TPC resolution parametrisation PbPb
519   if ( fResolutionCorrection ){
520     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
521     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
522   }
523   
524   //TOF resolution
525   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
526   
527 }
528
529 //______________________________________________________________________________
530 void AliPIDResponse::ExecNewRun()
531 {
532   //
533   // Things to Execute upon a new run
534   //
535   SetRecoInfo();
536   
537   SetITSParametrisation();
538   
539   SetTPCPidResponseMaster();
540   SetTPCParametrisation();
541
542   SetTRDPidResponseMaster(); 
543   InitializeTRDResponse();
544
545   SetEMCALPidResponseMaster(); 
546   InitializeEMCALResponse();
547   
548   fTOFResponse.SetTimeResolution(fTOFres);
549 }
550
551 //_____________________________________________________
552 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
553 {
554   //
555   // Get TPC multiplicity in bins of 150
556   //
557   
558   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
559   Double_t tpcMulti=0.;
560   if(vertexTPC){
561     Double_t vertexContribTPC=vertexTPC->GetNContributors();
562     tpcMulti=vertexContribTPC/150.;
563     if (tpcMulti>20.) tpcMulti=20.;
564   }
565   
566   return tpcMulti;
567 }
568
569 //______________________________________________________________________________
570 void AliPIDResponse::SetRecoInfo()
571 {
572   //
573   // Set reconstruction information
574   //
575   
576   //reset information
577   fLHCperiod="";
578   fMCperiodTPC="";
579   
580   fBeamType="";
581     
582   fBeamType="PP";
583   
584   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
585   //find the period by run number (UGLY, but not stored in ESD and AOD... )
586   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
587   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
588   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
589   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
590   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
591   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
592   else if (fRun>=136851&&fRun<=139517) {
593     fLHCperiod="LHC10H";
594     fMCperiodTPC="LHC10H8";
595     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
596     fBeamType="PBPB";
597   }
598   else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
599
600   //exception new pp MC productions from 2011
601   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
602 }
603
604 //______________________________________________________________________________
605 void AliPIDResponse::SetITSParametrisation()
606 {
607   //
608   // Set the ITS parametrisation
609   //
610 }
611
612 //______________________________________________________________________________
613 void AliPIDResponse::SetTPCPidResponseMaster()
614 {
615   //
616   // Load the TPC pid response functions from the OADB
617   //
618   //don't load twice for the moment
619    if (fArrPidResponseMaster) return;
620  
621
622   //reset the PID response functions
623   delete fArrPidResponseMaster;
624   fArrPidResponseMaster=0x0;
625   
626   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
627   
628   TFile *f=TFile::Open(fileName.Data());
629   if (f && f->IsOpen() && !f->IsZombie()){
630     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
631   }
632   delete f;
633   
634   if (!fArrPidResponseMaster){
635     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
636     return;
637   }
638   fArrPidResponseMaster->SetOwner();
639 }
640
641 //______________________________________________________________________________
642 void AliPIDResponse::SetTPCParametrisation()
643 {
644   //
645   // Change BB parametrisation for current run
646   //
647   
648   if (fLHCperiod.IsNull()) {
649     AliFatal("No period set, not changing parametrisation");
650     return;
651   }
652   
653   //
654   // Set default parametrisations for data and MC
655   //
656   
657   //data type
658   TString datatype="DATA";
659   //in case of mc fRecoPass is per default 1
660   if (fIsMC) {
661     datatype="MC";
662     fRecoPass=1;
663   }
664   
665   //
666   //reset old splines
667   //
668   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
669     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
670   }
671   
672   //
673   //set the new PID splines
674   //
675   TString period=fLHCperiod;
676   if (fArrPidResponseMaster){
677     TObject *grAll=0x0;
678     //for MC don't use period information
679 //     if (fIsMC) period="[A-Z0-9]*";
680     //for MC use MC period information
681     if (fIsMC) period=fMCperiodTPC;
682 //pattern for the default entry (valid for all particles)
683     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
684     
685     //loop over entries and filter them
686     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
687       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
688       if (responseFunction==0x0) continue;
689       TString responseName=responseFunction->GetName();
690       
691       if (!reg.MatchB(responseName)) continue;
692       
693       TObjArray *arr=reg.MatchS(responseName);
694       TString particleName=arr->At(1)->GetName();
695       delete arr;
696       if (particleName.IsNull()) continue;
697       if (particleName=="ALL") grAll=responseFunction;
698       else {
699         //find particle id
700         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
701           TString particle=AliPID::ParticleName(ispec);
702           particle.ToUpper();
703           if ( particle == particleName ){
704             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
705             fTPCResponse.SetUseDatabase(kTRUE);
706             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
707             break;
708           }
709         }
710       }
711     }
712     
713     //set default response function to all particles which don't have a specific one
714     if (grAll){
715       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
716         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
717           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
718           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
719         }
720       }
721     }
722   }
723   
724   //
725   // Setup resolution parametrisation
726   //
727   
728   //default
729   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
730   
731   if (fRun>=122195){
732     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
733   }
734   if (fArrPidResponseMaster)
735   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
736   
737   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
738 }
739
740 //______________________________________________________________________________
741 void AliPIDResponse::SetTRDPidResponseMaster()
742 {
743   //
744   // Load the TRD pid params and references from the OADB
745   //
746   if(fTRDPIDParams) return;
747   AliOADBContainer contParams("contParams"); 
748
749   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
750   if(statusPars){
751     AliError("Failed initializing PID Params from OADB");
752   } else {
753     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
754     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
755     if(!fTRDPIDParams){
756       AliError(Form("TRD Params not found in run %d", fRun));
757     }
758   }
759
760   AliOADBContainer contRefs("contRefs");
761   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
762   if(statusRefs){
763     AliInfo("Failed Loading References for TRD");
764   } else {
765     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
766     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
767     if(!fTRDPIDReference){
768       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
769     }
770   }
771 }
772
773 //______________________________________________________________________________
774 void AliPIDResponse::InitializeTRDResponse(){
775   //
776   // Set PID Params and references to the TRD PID response
777   // 
778   fTRDResponse.SetPIDParams(fTRDPIDParams);
779   fTRDResponse.Load(fTRDPIDReference);
780   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
781     fTRDslicesForPID[0] = 0;
782     fTRDslicesForPID[1] = 7;
783   }
784 }
785
786 //_________________________________________________________________________
787 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
788   //
789   // Check whether track is identified as electron under a given electron efficiency hypothesis
790   //
791   Double_t probs[AliPID::kSPECIES];
792   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
793
794   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
795   // Take mean of the TRD momenta in the given tracklets
796   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
797   Int_t nmomenta = 0;
798   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
799     if(vtrack->GetTRDmomentum(iPl) > 0.){
800       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
801     }
802   }
803   p = TMath::Mean(nmomenta, trdmomenta);
804
805   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
806 }
807
808 //______________________________________________________________________________
809 void AliPIDResponse::SetEMCALPidResponseMaster()
810 {
811   //
812   // Load the EMCAL pid response functions from the OADB
813   //
814   TObjArray* fEMCALPIDParamsRun      = NULL;
815   TObjArray* fEMCALPIDParamsPass     = NULL;
816
817   if(fEMCALPIDParams) return;
818   AliOADBContainer contParams("contParams"); 
819
820   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
821   if(statusPars){
822     AliError("Failed initializing PID Params from OADB");
823   } 
824   else {
825     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
826
827     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
828     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
829     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
830
831     if(!fEMCALPIDParams){
832       AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
833       AliInfo("Will take the standard LHC11a pass2 instead ...");
834
835       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(144871));
836       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",2)));
837       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
838
839       if(!fEMCALPIDParams){
840         AliError(Form("DEFAULT EMCAL Params (LHC11a pass2) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));       
841       }
842     }
843   }
844 }
845
846 //______________________________________________________________________________
847 void AliPIDResponse::InitializeEMCALResponse(){
848   //
849   // Set PID Params to the EMCAL PID response
850   // 
851   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
852
853 }