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