]>
Commit | Line | Data |
---|---|---|
773d3f8c | 1 | #include "TPDGCode.h" |
2 | #include "TF1.h" | |
3 | #include "TH1F.h" | |
4 | #include "TH2F.h" | |
c7cf2032 | 5 | #include "TProfile.h" |
773d3f8c | 6 | #include "TGraph.h" |
7 | #include "TGraphErrors.h" | |
8 | ||
9 | #include <TClonesArray.h> | |
10 | #include <TObjArray.h> | |
11 | #include <TList.h> | |
12 | ||
c7cf2032 | 13 | // #include "AliPID.h" |
773d3f8c | 14 | #include "AliESDEvent.h" |
15 | #include "AliESDInputHandler.h" | |
16 | #include "AliTrackReference.h" | |
17 | ||
c7cf2032 | 18 | #include "AliAnalysisTask.h" |
19 | // #include "AliAnalysisManager.h" | |
773d3f8c | 20 | |
c7cf2032 | 21 | #include "AliTRDtrackerV1.h" |
773d3f8c | 22 | #include "AliTRDtrackV1.h" |
c7cf2032 | 23 | #include "AliTRDcluster.h" |
773d3f8c | 24 | #include "AliTRDReconstructor.h" |
25 | #include "AliCDBManager.h" | |
c7cf2032 | 26 | // #include "../Cal/AliTRDCalPID.h" |
773d3f8c | 27 | |
28 | ||
29 | #include "AliTRDpidChecker.h" | |
30 | #include "AliTRDtrackInfo/AliTRDtrackInfo.h" | |
31 | ||
32 | // calculate pion efficiency at 90% electron efficiency for 11 momentum bins | |
33 | // this task should be used with simulated data only | |
34 | ||
35 | ClassImp(AliTRDpidChecker) | |
36 | ||
37 | //________________________________________________________________________ | |
3d86166d | 38 | AliTRDpidChecker::AliTRDpidChecker() |
39 | :AliTRDrecoTask("PID", "PID Checker") | |
773d3f8c | 40 | ,fReconstructor(0x0) |
773d3f8c | 41 | { |
42 | // | |
43 | // Default constructor | |
44 | // | |
45 | ||
46 | fReconstructor = new AliTRDReconstructor(); | |
47 | fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam()); | |
773d3f8c | 48 | } |
49 | ||
50 | ||
51 | //________________________________________________________________________ | |
52 | AliTRDpidChecker::~AliTRDpidChecker() | |
53 | { | |
3d86166d | 54 | if(fReconstructor) delete fReconstructor; |
773d3f8c | 55 | } |
56 | ||
57 | ||
58 | //________________________________________________________________________ | |
59 | void AliTRDpidChecker::CreateOutputObjects() | |
60 | { | |
61 | // Create histograms | |
62 | // Called once | |
63 | ||
64 | OpenFile(0, "RECREATE"); | |
3d86166d | 65 | fContainer = new TObjArray(); |
c7cf2032 | 66 | fContainer -> Expand(kGraphNNerr + 1); |
773d3f8c | 67 | |
c7cf2032 | 68 | const Float_t epsilon = 1/(2*(kBins-1)); // get nice histos with bin center at 0 and 1 |
773d3f8c | 69 | |
70 | // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method | |
773d3f8c | 71 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ |
72 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 73 | fContainer->AddAt(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon), kLQlikelihood + iPart*AliTRDCalPID::kNMom + iMom); |
c4c5bbfb | 74 | if(fDebugLevel>=4) Printf("LQ Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kLQlikelihood + iPart*AliTRDCalPID::kNMom + iMom); |
773d3f8c | 75 | } |
76 | } | |
77 | ||
1a6d7c9a | 78 | // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method |
773d3f8c | 79 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ |
80 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 81 | fContainer->AddAt(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon), kNNlikelihood + iPart*AliTRDCalPID::kNMom + iMom); |
c4c5bbfb | 82 | if(fDebugLevel>=4) Printf("NN Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kNNlikelihood + iPart*AliTRDCalPID::kNMom + iMom); |
1a6d7c9a | 83 | } |
84 | } | |
85 | ||
86 | // histos of the dE/dx distribution for all 5 particle species and 11 momenta | |
87 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
88 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 89 | fContainer->AddAt(new TH1F(Form("dEdx%d_%d",iPart,iMom), Form("dEdx distribution for %d_%d", iPart, iMom), 200, 0, 10000), kdEdx + iPart*AliTRDCalPID::kNMom + iMom); |
773d3f8c | 90 | } |
91 | } | |
92 | ||
c7cf2032 | 93 | // histos of the pulse height distribution for all 5 particle species and 11 momenta |
94 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
95 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
96 | fContainer->AddAt(new TProfile(Form("PH%d_%d",iPart,iMom), Form("PH distribution for %d_%d", iPart, iMom), AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5), kPH + iPart*AliTRDCalPID::kNMom + iMom); | |
97 | } | |
98 | } | |
99 | ||
100 | // momentum distributions - absolute and in momentum bins | |
101 | fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum); | |
102 | fContainer->AddAt(new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5),kMomentumBin); | |
103 | ||
c4c5bbfb | 104 | |
105 | // TGraph of the pion efficiencies | |
106 | ||
107 | TGraph *gEffisLQ = new TGraph(AliTRDCalPID::kNMom); | |
108 | gEffisLQ -> SetLineColor(4); | |
109 | gEffisLQ -> SetMarkerColor(4); | |
110 | gEffisLQ -> SetMarkerStyle(29); | |
111 | ||
112 | TGraphErrors *gEffisErrLQ = new TGraphErrors(AliTRDCalPID::kNMom); | |
113 | gEffisErrLQ -> SetLineColor(4); | |
114 | gEffisErrLQ -> SetMarkerColor(4); | |
115 | gEffisErrLQ -> SetMarkerStyle(29); | |
116 | ||
117 | TGraph *gEffisNN = new TGraph(AliTRDCalPID::kNMom); | |
118 | gEffisNN -> SetLineColor(2); | |
119 | gEffisNN -> SetMarkerColor(2); | |
120 | gEffisNN -> SetMarkerStyle(29); | |
121 | ||
122 | TGraphErrors *gEffisErrNN = new TGraphErrors(AliTRDCalPID::kNMom); | |
123 | gEffisErrNN -> SetLineColor(2); | |
124 | gEffisErrNN -> SetMarkerColor(2); | |
125 | gEffisErrNN -> SetMarkerStyle(29); | |
126 | ||
127 | ||
128 | fContainer -> AddAt(gEffisLQ,kGraphLQ); | |
129 | fContainer -> AddAt(gEffisErrLQ,kGraphLQerr); | |
130 | fContainer -> AddAt(gEffisNN,kGraphNN); | |
131 | fContainer -> AddAt(gEffisErrNN,kGraphNNerr); | |
773d3f8c | 132 | } |
133 | ||
134 | //________________________________________________________________________ | |
135 | void AliTRDpidChecker::Exec(Option_t *) | |
136 | { | |
137 | // Main loop | |
138 | // Called for each event | |
139 | ||
140 | ||
141 | // if(!AliTracker::GetFieldMap()){ | |
142 | // AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG); | |
143 | // AliTracker::SetFieldMap(field, kTRUE); | |
144 | // } | |
145 | ||
773d3f8c | 146 | TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom]; |
147 | TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom]; | |
1a6d7c9a | 148 | TH1F *hdEdx[AliPID::kSPECIES][AliTRDCalPID::kNMom]; |
c7cf2032 | 149 | TProfile *hPH[AliPID::kSPECIES][AliTRDCalPID::kNMom]; |
773d3f8c | 150 | |
151 | for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){ | |
152 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 153 | hPIDLQ[iPart][iMom] = (TH1F*)fContainer->At(kLQlikelihood + iPart*AliTRDCalPID::kNMom+iMom); |
154 | hPIDNN[iPart][iMom] = (TH1F*)fContainer->At(kNNlikelihood + iPart*AliTRDCalPID::kNMom+iMom); | |
155 | hdEdx[iPart][iMom] = (TH1F*)fContainer->At(kdEdx + iPart*AliTRDCalPID::kNMom+iMom); | |
156 | hPH[iPart][iMom] = (TProfile*)fContainer->At(kPH + iPart*AliTRDCalPID::kNMom+iMom); | |
773d3f8c | 157 | } |
158 | } | |
c4c5bbfb | 159 | |
c7cf2032 | 160 | TH1F *hMom = (TH1F*)fContainer->At(kMomentum); |
161 | TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin); | |
1a6d7c9a | 162 | |
773d3f8c | 163 | Int_t labelsacc[10000]; |
164 | memset(labelsacc, 0, sizeof(Int_t) * 10000); | |
165 | ||
166 | Float_t mom; | |
167 | ULong_t status; | |
168 | Int_t nTRD = 0; | |
1a6d7c9a | 169 | Float_t *fdEdx; |
773d3f8c | 170 | |
171 | AliTRDtrackInfo *track = 0x0; | |
c7cf2032 | 172 | AliTRDtrackV1 *TRDtrack = 0x0; |
773d3f8c | 173 | AliTrackReference *ref = 0x0; |
174 | AliExternalTrackParam *esd = 0x0; | |
1a6d7c9a | 175 | |
5d6dc395 | 176 | AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer]; |
177 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) | |
1a6d7c9a | 178 | TRDtracklet[iChamb] = 0x0; |
179 | ||
c7cf2032 | 180 | AliTRDcluster *TRDcluster = 0x0; |
181 | ||
773d3f8c | 182 | for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){ |
183 | track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk); | |
184 | if(!track->HasESDtrack()) continue; | |
185 | status = track->GetStatus(); | |
186 | if(!(status&AliESDtrack::kTPCout)) continue; | |
187 | ||
188 | if(!(TRDtrack = track->GetTRDtrack())) continue; | |
189 | //&&(track->GetNumberOfClustersRefit() | |
1a6d7c9a | 190 | |
191 | // use only tracks that hit 6 chambers | |
5d6dc395 | 192 | if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue; |
c4c5bbfb | 193 | |
773d3f8c | 194 | ref = track->GetTrackRef(0); |
195 | esd = track->GetOuterParam(); | |
196 | mom = ref ? ref->P(): esd->P(); | |
197 | ||
198 | labelsacc[nTRD] = track->GetLabel(); | |
199 | nTRD++; | |
200 | ||
773d3f8c | 201 | |
202 | // set the 11 momentum bins | |
203 | Int_t iMomBin = -1; | |
204 | if(mom < .7) iMomBin = 0; | |
205 | else if(mom < .9) iMomBin = 1; | |
206 | else if(mom < 1.25) iMomBin = 2; | |
207 | else if(mom < 1.75) iMomBin = 3; | |
208 | else if(mom < 2.5) iMomBin = 4; | |
209 | else if(mom < 3.5) iMomBin = 5; | |
210 | else if(mom < 4.5) iMomBin = 6; | |
211 | else if(mom < 5.5) iMomBin = 7; | |
212 | else if(mom < 7.0) iMomBin = 8; | |
213 | else if(mom < 9.0) iMomBin = 9; | |
214 | else iMomBin = 10; | |
215 | ||
c7cf2032 | 216 | // fill momentum histo to have the momentum distribution |
217 | hMom -> Fill(mom); | |
218 | hMomBin -> Fill(iMomBin); | |
219 | ||
220 | ||
773d3f8c | 221 | // set the reconstructor |
222 | TRDtrack -> SetReconstructor(fReconstructor); | |
223 | ||
224 | ||
c7cf2032 | 225 | // if no monte carlo data available -> use TRDpid |
226 | if(!HasMCdata()){ | |
227 | fReconstructor -> SetOption("nn"); | |
228 | TRDtrack -> CookPID(); | |
229 | if(TRDtrack -> GetPID(0) > TRDtrack -> GetPID(1) + TRDtrack -> GetPID(2) + TRDtrack -> GetPID(3) + TRDtrack -> GetPID(4)){ | |
c4c5bbfb | 230 | track -> SetPDG(kElectron); |
231 | } | |
232 | else if(TRDtrack -> GetPID(4) > TRDtrack -> GetPID(2) && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(3) && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(1)){ | |
233 | track -> SetPDG(kProton); | |
234 | } | |
235 | else if(TRDtrack -> GetPID(3) > TRDtrack -> GetPID(1) && TRDtrack -> GetPID(3) > TRDtrack -> GetPID(2)){ | |
236 | track -> SetPDG(kKPlus); | |
237 | } | |
238 | else if(TRDtrack -> GetPID(1) > TRDtrack -> GetPID(2)){ | |
239 | track -> SetPDG(kMuonPlus); | |
240 | } | |
241 | else{ | |
242 | track -> SetPDG(kPiPlus); | |
c7cf2032 | 243 | } |
244 | } | |
245 | ||
246 | ||
247 | // calculate the probabilities for electron probability using 2-dim LQ, the deposited charge per chamber and the pulse height spectra and fill histograms | |
773d3f8c | 248 | fReconstructor -> SetOption("!nn"); |
249 | TRDtrack -> CookPID(); | |
c4c5bbfb | 250 | |
251 | if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] LQLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0)); | |
252 | ||
253 | ||
5d6dc395 | 254 | Float_t SumdEdx[AliTRDgeometry::kNlayer]; |
255 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ | |
1a6d7c9a | 256 | TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb); |
257 | SumdEdx[iChamb] = 0.; | |
258 | fdEdx = TRDtracklet[iChamb] -> GetdEdx(); | |
259 | SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2]; | |
260 | } | |
261 | ||
773d3f8c | 262 | switch(track->GetPDG()){ |
263 | case kElectron: | |
264 | case kPositron: | |
265 | hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
5d6dc395 | 266 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
1a6d7c9a | 267 | hdEdx[AliPID::kElectron][iMomBin] -> Fill(SumdEdx[iChamb]); |
5d6dc395 | 268 | for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){ |
269 | if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus))) | |
270 | continue; | |
271 | hPH[AliPID::kElectron][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus)); | |
272 | } | |
1a6d7c9a | 273 | } |
773d3f8c | 274 | break; |
275 | case kMuonPlus: | |
276 | case kMuonMinus: | |
277 | hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
5d6dc395 | 278 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
1a6d7c9a | 279 | hdEdx[AliPID::kMuon][iMomBin] -> Fill(SumdEdx[iChamb]); |
5d6dc395 | 280 | for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){ |
281 | if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus))) | |
282 | continue; | |
283 | hPH[AliPID::kMuon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus)); | |
284 | } | |
1a6d7c9a | 285 | } |
773d3f8c | 286 | break; |
287 | case kPiPlus: | |
288 | case kPiMinus: | |
289 | hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
5d6dc395 | 290 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
1a6d7c9a | 291 | hdEdx[AliPID::kPion][iMomBin] -> Fill(SumdEdx[iChamb]); |
5d6dc395 | 292 | for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){ |
293 | if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus))) | |
294 | continue; | |
295 | hPH[AliPID::kPion][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus)); | |
296 | } | |
1a6d7c9a | 297 | } |
773d3f8c | 298 | break; |
299 | case kKPlus: | |
300 | case kKMinus: | |
301 | hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
5d6dc395 | 302 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
1a6d7c9a | 303 | hdEdx[AliPID::kKaon][iMomBin] -> Fill(SumdEdx[iChamb]); |
5d6dc395 | 304 | for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){ |
305 | if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus))) | |
306 | continue; | |
307 | hPH[AliPID::kKaon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus)); | |
308 | } | |
1a6d7c9a | 309 | } |
773d3f8c | 310 | break; |
311 | case kProton: | |
312 | case kProtonBar: | |
1a6d7c9a | 313 | hPIDLQ[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0)); |
5d6dc395 | 314 | for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){ |
1a6d7c9a | 315 | hdEdx[AliPID::kProton][iMomBin] -> Fill(SumdEdx[iChamb]); |
5d6dc395 | 316 | for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){ |
317 | if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus))) | |
318 | continue; | |
319 | hPH[AliPID::kProton][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus)); | |
320 | } | |
1a6d7c9a | 321 | } |
773d3f8c | 322 | break; |
323 | } | |
324 | ||
325 | ||
326 | // calculate the probabilities and fill histograms for electrons using NN | |
327 | fReconstructor -> SetOption("nn"); | |
328 | TRDtrack->CookPID(); | |
c4c5bbfb | 329 | |
330 | ||
331 | if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] NNLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0)); | |
332 | ||
333 | ||
773d3f8c | 334 | switch(track->GetPDG()){ |
335 | case kElectron: | |
336 | case kPositron: | |
337 | hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
338 | break; | |
339 | case kMuonPlus: | |
340 | case kMuonMinus: | |
341 | hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
342 | break; | |
343 | case kPiPlus: | |
344 | case kPiMinus: | |
345 | hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
346 | break; | |
347 | case kKPlus: | |
348 | case kKMinus: | |
349 | hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0)); | |
350 | break; | |
351 | case kProton: | |
352 | case kProtonBar: | |
1a6d7c9a | 353 | hPIDNN[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0)); |
773d3f8c | 354 | break; |
355 | } | |
356 | } | |
357 | ||
3d86166d | 358 | PostData(0, fContainer); |
773d3f8c | 359 | } |
360 | ||
c4c5bbfb | 361 | //________________________________________________________ |
28efdace | 362 | void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t */*opt*/) |
c4c5bbfb | 363 | { |
364 | switch(ifig){ | |
365 | case 0: | |
366 | first = (Int_t)kGraphStart+1; last = first+1; | |
367 | break; | |
368 | case 1: | |
369 | first = (Int_t)kGraphStart+3; last = first+1; | |
370 | break; | |
371 | default: | |
372 | first = (Int_t)kGraphStart; last = first; | |
373 | break; | |
374 | } | |
375 | } | |
376 | ||
377 | ||
773d3f8c | 378 | //________________________________________________________________________ |
d85cd79c | 379 | Bool_t AliTRDpidChecker::PostProcess() |
773d3f8c | 380 | { |
c7cf2032 | 381 | // Draw result to the screen |
382 | // Called once at the end of the query | |
383 | ||
3d86166d | 384 | if (!fContainer) { |
773d3f8c | 385 | Printf("ERROR: list not available"); |
d85cd79c | 386 | return kFALSE; |
773d3f8c | 387 | } |
c7cf2032 | 388 | // return kTRUE; // testing protection |
773d3f8c | 389 | |
1a6d7c9a | 390 | // normalize the dE/dx histos |
391 | const Int_t kNSpectra = AliPID::kSPECIES*AliTRDCalPID::kNMom; | |
392 | TH1F *hdEdx[kNSpectra]; | |
393 | for(Int_t iHist = 0; iHist < kNSpectra; iHist++){ | |
c7cf2032 | 394 | hdEdx[iHist] = (TH1F*)fContainer->At(kdEdx + iHist); |
1a6d7c9a | 395 | if(hdEdx[iHist] -> GetEntries() > 0) |
396 | hdEdx[iHist] -> Scale(1./hdEdx[iHist] -> GetEntries()); | |
397 | else continue; | |
398 | } | |
399 | ||
400 | ||
773d3f8c | 401 | // container for the pion efficiencies and the errors |
402 | Double_t PionEffiLQ[AliTRDCalPID::kNMom], | |
403 | PionEffiNN[AliTRDCalPID::kNMom], | |
404 | PionEffiErrorLQ[AliTRDCalPID::kNMom], | |
405 | PionEffiErrorNN[AliTRDCalPID::kNMom]; | |
406 | ||
407 | ||
408 | // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ) | |
409 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 410 | PionEffiLQ[iMom] = GetPionEfficiency((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom)); |
411 | PionEffiErrorLQ[iMom] = GetError((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom)); | |
412 | if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]); | |
773d3f8c | 413 | } |
414 | ||
415 | // calculate the pion efficiencies and the errors for 90% electron efficiency (NN) | |
416 | for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){ | |
c7cf2032 | 417 | PionEffiNN[iMom] = GetPionEfficiency((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom)); |
418 | PionEffiErrorNN[iMom] = GetError((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom)); | |
419 | if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]); | |
773d3f8c | 420 | } |
421 | ||
422 | ||
423 | // create TGraph to plot the pion efficiencies | |
424 | TGraph *gEffisLQ=0x0, *gEffisNN=0x0; | |
425 | TGraphErrors *gEffisErrLQ=0x0, *gEffisErrNN=0x0; | |
426 | ||
c7cf2032 | 427 | gEffisLQ = (TGraph*)fContainer->At(kGraphLQ); |
428 | gEffisErrLQ = (TGraphErrors*)fContainer->At(kGraphLQerr); | |
429 | gEffisNN = (TGraph*)fContainer->At(kGraphNN); | |
430 | gEffisErrNN = (TGraphErrors*)fContainer->At(kGraphNNerr); | |
773d3f8c | 431 | |
432 | for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){ | |
433 | Float_t momentum = AliTRDCalPID::GetMomentum(iBin); | |
434 | gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]); | |
435 | gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]); | |
436 | gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]); | |
437 | ||
438 | gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]); | |
439 | gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]); | |
440 | gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]); | |
441 | } | |
442 | ||
443 | gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method"); | |
444 | gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method"); | |
445 | gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method"); | |
446 | gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method"); | |
c7cf2032 | 447 | |
c4c5bbfb | 448 | fNRefFigures = 2; |
c7cf2032 | 449 | return kTRUE; // testing protection |
d85cd79c | 450 | } |
451 | ||
c7cf2032 | 452 | |
d85cd79c | 453 | //________________________________________________________________________ |
454 | void AliTRDpidChecker::Terminate(Option_t *) | |
455 | { | |
456 | // Draw result to the screen | |
457 | // Called once at the end of the query | |
458 | ||
459 | fContainer = dynamic_cast<TObjArray*>(GetOutputData(0)); | |
460 | if (!fContainer) { | |
461 | Printf("ERROR: list not available"); | |
462 | return; | |
463 | } | |
773d3f8c | 464 | } |
465 | ||
466 | ||
467 | //________________________________________________________________________ | |
468 | Double_t AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){ | |
469 | ||
c7cf2032 | 470 | const Float_t Multi = 0.9; // electron efficiency |
773d3f8c | 471 | Int_t abin, bbin; |
472 | Double_t SumElecs, SumPions; // integrated sum of elecs and pions in the histos | |
473 | Double_t aBinSum, bBinSum; // content of a single bin in the histos | |
474 | ||
3d86166d | 475 | TH1F *Histo1 = (TH1F*)fContainer->At(Index1); // electron histo |
476 | TH1F *Histo2 = (TH1F*)fContainer->At(Index2); // pion histo | |
773d3f8c | 477 | |
478 | ||
c7cf2032 | 479 | if(Index1 >= kNNlikelihood) // print the correct momentum index for neural networks |
480 | Index1 = Index1 - kNNlikelihood; | |
773d3f8c | 481 | SumElecs = 0.; |
482 | if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){ | |
c7cf2032 | 483 | if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetPionEfficiency : Warning: Histo momentum intervall %d has no Entries!", Index1); |
773d3f8c | 484 | return -1.; |
485 | } | |
486 | Histo1 -> Scale(1./Histo1->GetEntries()); | |
487 | Histo2 -> Scale(1./Histo2->GetEntries()); | |
488 | ||
489 | ||
490 | // calculate threshold for pion efficiency | |
491 | for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){ | |
492 | aBinSum = 0; | |
493 | aBinSum = Histo1 -> GetBinContent(abin); | |
494 | if(!(aBinSum == 0)){ | |
495 | SumElecs = SumElecs + aBinSum; | |
496 | } | |
497 | if (SumElecs >= Multi){ | |
498 | break; | |
499 | } | |
500 | } | |
501 | ||
502 | ||
503 | // calculate pion efficiency | |
504 | SumPions = 0.; | |
505 | for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){ | |
506 | bBinSum = 0; | |
507 | bBinSum = Histo2 -> GetBinContent(bbin); | |
508 | if(!(bBinSum == 0)){ | |
509 | SumPions = SumPions + bBinSum; | |
510 | } | |
511 | } | |
512 | ||
513 | ||
514 | // print the electron efficiency and its cuts | |
c7cf2032 | 515 | if(fDebugLevel>=1) Printf("Cut for momentum intervall %d and electron efficiency of %f at: %5.3f", Index1, SumElecs, Histo1 -> GetBinCenter(abin)); |
516 | if(fDebugLevel>=1) Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries()); | |
773d3f8c | 517 | |
518 | ||
519 | // return the pion efficiency | |
520 | return SumPions; | |
521 | ||
522 | } | |
523 | ||
524 | ||
525 | //________________________________________________________________________ | |
526 | Double_t AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){ | |
527 | ||
528 | ||
c7cf2032 | 529 | // const Int_t kBins = 12000; // binning of the histos and eficiency calculation |
773d3f8c | 530 | const Float_t Multi = 0.9; // electron efficiency |
531 | Int_t abinE, bbinE, cbinE; | |
532 | Double_t SumElecsE[kBins], SumPionsE[kBins]; // array of the integrated sum in each bin | |
533 | Double_t aBinSumE, bBinSumE; // content of a single bin | |
534 | Double_t EleEffi, PioEffi; // electron and pion efficiency | |
535 | Bool_t bFirst = 1; // checks if threshold is crossed for the first time | |
536 | Double_t fError = 0.; // error of the pion efficiency | |
537 | ||
538 | ||
3d86166d | 539 | TH1F *Histo1 = (TH1F*)fContainer->At(Index1); // electron histo |
540 | TH1F *Histo2 = (TH1F*)fContainer->At(Index2); // pion histo | |
773d3f8c | 541 | |
c7cf2032 | 542 | if(Index1 > kNNlikelihood) // print the correct momentum index for neural networks |
543 | Index1 = Index1 - kNNlikelihood; | |
773d3f8c | 544 | if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){ |
c7cf2032 | 545 | if(fDebugLevel>=2) Printf("Warning: Histo momentum intervall %d has no Entries!", Index1); |
773d3f8c | 546 | return -1.; |
547 | } | |
548 | ||
549 | for(Int_t iBin = 0; iBin < kBins; iBin++){ | |
550 | SumElecsE[iBin] = 0.; | |
551 | SumPionsE[iBin] = 0.; | |
552 | } | |
553 | ||
554 | EleEffi = 0.; | |
555 | PioEffi = 0.; | |
556 | cbinE = -1; | |
557 | ||
558 | ||
559 | // calculate electron efficiency of each bin | |
560 | for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){ | |
561 | aBinSumE = 0; | |
562 | aBinSumE = Histo1 -> GetBinContent(abinE); | |
563 | ||
564 | SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE; | |
565 | if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){ | |
566 | bFirst = 0; | |
567 | cbinE = abinE; | |
568 | EleEffi = (SumElecsE[cbinE]); | |
569 | } | |
570 | } | |
571 | ||
572 | ||
573 | // calculate pion efficiency of each bin | |
574 | for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){ | |
575 | bBinSumE = 0; | |
576 | bBinSumE = Histo2 -> GetBinContent(bbinE); | |
577 | ||
578 | SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE; | |
579 | if(bbinE == cbinE){ | |
580 | PioEffi = (SumPionsE[cbinE]); | |
581 | } | |
582 | } | |
583 | ||
584 | ||
585 | // pion efficiency vs electron efficiency | |
586 | TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE); | |
587 | ||
588 | // use fit function to get derivate of the TGraph for error calculation | |
589 | TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", Multi-.05, Multi+.05); | |
590 | gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05); | |
c7cf2032 | 591 | if(fDebugLevel>=1) Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi)); |
773d3f8c | 592 | |
593 | // return the error of the pion efficiency | |
1a6d7c9a | 594 | if(((1.-PioEffi) < 0) || ((1.-EleEffi) < 0)){ |
c7cf2032 | 595 | if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetError : Warning: EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!"); |
1a6d7c9a | 596 | return -1; |
597 | } | |
773d3f8c | 598 | fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi))); |
599 | return fError; | |
600 | } | |
c7cf2032 | 601 |