Sqrt(2) added in the combined nsigma
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraAODPID.cxx
1 #include "AliSpectraAODPID.h"
2 #include "AliAODEvent.h"      
3 #include "TH1F.h"             
4 #include "TH2F.h"             
5 #include "TList.h"            
6 #include "AliAODTrack.h"      
7 #include "AliAODMCParticle.h" 
8 #include "AliPIDResponse.h"   
9 #include "AliAnalysisManager.h"
10 #include "AliInputEventHandler.h"
11 #include "AliSpectraAODTrackCuts.h"
12
13 ClassImp(AliSpectraAODPID)
14
15 AliSpectraAODPID::AliSpectraAODPID() : TNamed("PID", "PID object"), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0) {
16
17 }
18
19 AliSpectraAODPID::AliSpectraAODPID(AODPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0) {
20
21
22
23 }
24
25
26
27 void AliSpectraAODPID::FillQAHistos(AliSpectraAODHistoManager * hman, AliAODTrack * track, AliSpectraAODTrackCuts * trackCuts) {
28
29   // fill a bunch of QA histos
30
31
32   // Get PID response object, if needed
33   if(!fPIDResponse) {
34     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
35     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
36     fPIDResponse = inputHandler->GetPIDResponse();
37   }
38   
39   //Response
40   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);
41
42   hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
43
44
45   Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
46   Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); 
47   Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); 
48   Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;
49
50   if(track->Pt()>trackCuts->GetPtTOFMatching()){
51     hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
52     nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
53     nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); 
54     nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); 
55     
56     //TOF
57     hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
58     hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
59     hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
60     hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
61     hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
62     hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
63   }
64   
65   Double_t nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2);
66   Double_t nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2);
67   Double_t nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2);
68   
69   //TPC
70   hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
71   hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
72   hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
73   hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
74   hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));
75   hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));
76   //TPCTOF
77   hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
78   hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
79   hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
80   hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
81   hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
82   hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
83
84 }
85
86 Int_t AliSpectraAODPID::GetParticleSpecie(AliAODMCParticle * part) {
87   // return PID according to MC truth
88   switch(TMath::Abs(part->PdgCode())){
89   case 2212:
90     return kSpProton;
91     break;
92   case 321:
93     return kSpKaon;
94     break;
95   case 211:
96     return kSpPion;
97     break;
98   default:
99     return kSpUndefined;
100   } 
101 }
102
103
104 Int_t AliSpectraAODPID::GetParticleSpecie(AliAODTrack      * trk, AliSpectraAODTrackCuts * trackCuts) {
105   // return PID according to detectors
106   
107   // Get PID response object, if needed
108   if(!fPIDResponse) {
109     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
110     AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
111     fPIDResponse = inputHandler->GetPIDResponse();
112   }
113
114   if(!fPIDResponse) {
115     AliFatal("Cannot get pid response");
116     return 0;
117   }
118
119
120   // Compute nsigma for each hypthesis
121   AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
122
123   // --- TPC
124   Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));
125   Double_t nsigmaTPCkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); 
126   Double_t nsigmaTPCkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); 
127   // --- TOF
128   Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;
129   if(trk->Pt()>trackCuts->GetPtTOFMatching()){
130     nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
131     nsigmaTOFkKaon   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); 
132     nsigmaTOFkPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); 
133   }
134           
135   // --- combined
136   Double_t nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
137   Double_t nsigmaTPCTOFkKaon   = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
138   Double_t nsigmaTPCTOFkPion   = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
139
140
141   // select the nsigma to be used for the actual PID
142   Double_t nsigmaPion, nsigmaKaon, nsigmaProton;
143
144   switch (fPIDType) {
145   case kNSigmaTPC:
146     nsigmaProton  =  nsigmaTPCkProton;
147     nsigmaKaon    =  nsigmaTPCkKaon  ;
148     nsigmaPion    =  nsigmaTPCkPion  ;
149     break;
150   case kNSigmaTOF:
151     nsigmaProton  =  nsigmaTOFkProton;
152     nsigmaKaon    =  nsigmaTOFkKaon  ;
153     nsigmaPion    =  nsigmaTOFkPion  ;
154     break;
155   case kNSigmaTPCTOF:
156     nsigmaProton  =  nsigmaTPCTOFkProton;
157     nsigmaKaon    =  nsigmaTPCTOFkKaon  ;
158     nsigmaPion    =  nsigmaTPCTOFkPion  ;
159     break;
160   }
161
162   // guess the particle based on the smaller nsigma
163   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
164   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < fNSigmaPID)) return kSpKaon;
165   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < fNSigmaPID)) return kSpPion;
166   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) return kSpProton;
167
168   // else, return undefined
169   return kSpUndefined;
170
171 }
172
173 Long64_t AliSpectraAODPID::Merge(TCollection* list)
174 {
175   // Merging interface.
176   // Returns the number of merged objects (including this).
177
178   Printf("Merging");
179
180   if (!list)
181     return 0;
182
183   if (list->IsEmpty())
184     return 1;
185
186   TIterator* iter = list->MakeIterator();
187   TObject* obj;
188
189   // Actually, we don't do anything here...
190   // collections of all histograms
191   //  TList collections;
192
193   Int_t count = 0;
194
195   while ((obj = iter->Next())) {
196     AliSpectraAODPID* entry = dynamic_cast<AliSpectraAODPID*> (obj);
197     if (entry == 0) 
198       continue;
199
200     // TH1I * histo = entry->GetHistoCuts();      
201     // collections.Add(histo);
202     count++;
203   }
204   
205   //  fHistoCuts->Merge(&collections);
206   
207   delete iter;
208   Printf("OK");
209   return count+1;
210 }
211