]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/STEERBase/AliPIDResponse.cxx
Fix
[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&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
599   else if (fRun>=166529) {
600     fLHCperiod="LHC11H";
601     fMCperiodTPC="LHC11A10";
602     fBeamType="PBPB";
603   }
604
605
606   //exception new pp MC productions from 2011
607   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
608 }
609
610 //______________________________________________________________________________
611 void AliPIDResponse::SetITSParametrisation()
612 {
613   //
614   // Set the ITS parametrisation
615   //
616 }
617
618 //______________________________________________________________________________
619 void AliPIDResponse::SetTPCPidResponseMaster()
620 {
621   //
622   // Load the TPC pid response functions from the OADB
623   //
624   //don't load twice for the moment
625    if (fArrPidResponseMaster) return;
626  
627
628   //reset the PID response functions
629   delete fArrPidResponseMaster;
630   fArrPidResponseMaster=0x0;
631   
632   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
633   
634   TFile *f=TFile::Open(fileName.Data());
635   if (f && f->IsOpen() && !f->IsZombie()){
636     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
637   }
638   delete f;
639   
640   if (!fArrPidResponseMaster){
641     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
642     return;
643   }
644   fArrPidResponseMaster->SetOwner();
645 }
646
647 //______________________________________________________________________________
648 void AliPIDResponse::SetTPCParametrisation()
649 {
650   //
651   // Change BB parametrisation for current run
652   //
653   
654   if (fLHCperiod.IsNull()) {
655     AliFatal("No period set, not changing parametrisation");
656     return;
657   }
658   
659   //
660   // Set default parametrisations for data and MC
661   //
662   
663   //data type
664   TString datatype="DATA";
665   //in case of mc fRecoPass is per default 1
666   if (fIsMC) {
667     datatype="MC";
668     fRecoPass=1;
669   }
670   
671   //
672   //reset old splines
673   //
674   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
675     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
676   }
677   
678   //
679   //set the new PID splines
680   //
681   TString period=fLHCperiod;
682   if (fArrPidResponseMaster){
683     TObject *grAll=0x0;
684     //for MC don't use period information
685 //     if (fIsMC) period="[A-Z0-9]*";
686     //for MC use MC period information
687     if (fIsMC) period=fMCperiodTPC;
688 //pattern for the default entry (valid for all particles)
689     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
690     
691     //loop over entries and filter them
692     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
693       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
694       if (responseFunction==0x0) continue;
695       TString responseName=responseFunction->GetName();
696       
697       if (!reg.MatchB(responseName)) continue;
698       
699       TObjArray *arr=reg.MatchS(responseName);
700       TString particleName=arr->At(1)->GetName();
701       delete arr;
702       if (particleName.IsNull()) continue;
703       if (particleName=="ALL") grAll=responseFunction;
704       else {
705         //find particle id
706         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
707           TString particle=AliPID::ParticleName(ispec);
708           particle.ToUpper();
709           if ( particle == particleName ){
710             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
711             fTPCResponse.SetUseDatabase(kTRUE);
712             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
713             break;
714           }
715         }
716       }
717     }
718     
719     //set default response function to all particles which don't have a specific one
720     if (grAll){
721       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
722         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
723           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
724           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
725         }
726       }
727     }
728   }
729   
730   //
731   // Setup resolution parametrisation
732   //
733   
734   //default
735   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
736   
737   if (fRun>=122195){
738     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
739   }
740   if (fArrPidResponseMaster)
741   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
742   
743   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
744 }
745
746 //______________________________________________________________________________
747 void AliPIDResponse::SetTRDPidResponseMaster()
748 {
749   //
750   // Load the TRD pid params and references from the OADB
751   //
752   if(fTRDPIDParams) return;
753   AliOADBContainer contParams("contParams"); 
754
755   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
756   if(statusPars){
757     AliError("Failed initializing PID Params from OADB");
758   } else {
759     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
760     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
761     if(!fTRDPIDParams){
762       AliError(Form("TRD Params not found in run %d", fRun));
763     }
764   }
765
766   AliOADBContainer contRefs("contRefs");
767   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
768   if(statusRefs){
769     AliInfo("Failed Loading References for TRD");
770   } else {
771     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
772     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
773     if(!fTRDPIDReference){
774       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
775     }
776   }
777 }
778
779 //______________________________________________________________________________
780 void AliPIDResponse::InitializeTRDResponse(){
781   //
782   // Set PID Params and references to the TRD PID response
783   // 
784   fTRDResponse.SetPIDParams(fTRDPIDParams);
785   fTRDResponse.Load(fTRDPIDReference);
786   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
787     fTRDslicesForPID[0] = 0;
788     fTRDslicesForPID[1] = 7;
789   }
790 }
791
792 //_________________________________________________________________________
793 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
794   //
795   // Check whether track is identified as electron under a given electron efficiency hypothesis
796   //
797   Double_t probs[AliPID::kSPECIES];
798   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
799
800   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
801   // Take mean of the TRD momenta in the given tracklets
802   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
803   Int_t nmomenta = 0;
804   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
805     if(vtrack->GetTRDmomentum(iPl) > 0.){
806       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
807     }
808   }
809   p = TMath::Mean(nmomenta, trdmomenta);
810
811   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
812 }
813
814 //______________________________________________________________________________
815 void AliPIDResponse::SetEMCALPidResponseMaster()
816 {
817   //
818   // Load the EMCAL pid response functions from the OADB
819   //
820   TObjArray* fEMCALPIDParamsRun      = NULL;
821   TObjArray* fEMCALPIDParamsPass     = NULL;
822
823   if(fEMCALPIDParams) return;
824   AliOADBContainer contParams("contParams"); 
825
826   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
827   if(statusPars){
828     AliError("Failed initializing PID Params from OADB");
829   } 
830   else {
831     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
832
833     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
834     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
835     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
836
837     if(!fEMCALPIDParams){
838       AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
839       AliInfo("Will take the standard LHC11a pass2 instead ...");
840
841       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(144871));
842       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",2)));
843       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
844
845       if(!fEMCALPIDParams){
846         AliError(Form("DEFAULT EMCAL Params (LHC11a pass2) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));       
847       }
848     }
849   }
850 }
851
852 //______________________________________________________________________________
853 void AliPIDResponse::InitializeEMCALResponse(){
854   //
855   // Set PID Params to the EMCAL PID response
856   // 
857   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
858
859 }