]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/pid/AliAnalysisTaskPhiBayes.cxx
coverity
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskPhiBayes.cxx
CommitLineData
a8ad4709 1#include "AliAnalysisTaskPhiBayes.h"
2
3// ROOT includes
4#include <TMath.h>
5
6// AliRoot includes
7#include "AliInputEventHandler.h"
8#include "AliAODEvent.h"
9#include "AliAODVertex.h"
10#include "AliAODTrack.h"
11#include "AliESDtrack.h"
12#include "AliCentrality.h"
13#include "AliVHeader.h"
14#include "AliAODVZERO.h"
15#include "TFile.h"
16#include "AliOADBContainer.h"
17#include "TH2F.h"
18#include "TF1.h"
19#include "AliGenHijingEventHeader.h"
20#include "AliMCEvent.h"
21#include "AliAODMCHeader.h"
22#include "AliAODMCParticle.h"
23#include "TChain.h"
24#include "AliESDtrackCuts.h"
25#include "AliESDVertex.h"
26#include "AliEventplane.h"
27#include "AliAnalysisManager.h"
28#include "TRandom.h"
29
30
31Float_t AliAnalysisTaskPhiBayes::fPtPhiMin[AliAnalysisTaskPhiBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
32Float_t AliAnalysisTaskPhiBayes::fPtPhiMax[AliAnalysisTaskPhiBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
33
34// STL includes
35//#include <iostream>
36//using namespace std;
37
38ClassImp(AliAnalysisTaskPhiBayes)
39//_____________________________________________________________________________
40AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes():
41 AliAnalysisTaskSE(),
42 fVtxCut(10.0), // cut on |vertex| < fVtxCut
43 fEtaCut(0.8), // cut on |eta| < fEtaCut
44 fMinPt(0.15), // cut on pt > fMinPt
45 fIsMC(kFALSE),
46 fQAsw(kFALSE),
47 fNcluster(70),
48 fFilterBit(1),
49 fList(new TList()),
50 fList2(new TList()),
51 fList3(new TList()),
52 fCentrality(-1),
53 fPsi(0),
54 fPtPhi(0.),
55 fPhiPhi(0.),
56 fEtaPhi(0.),
57 fMassV0(0.),
58 fPtKp(0.),
59 fPhiKp(0.),
60 fEtaKp(0.),
61 fPtKn(0.),
62 fPhiKn(0.),
63 fEtaKn(0.),
64 fPidKp(0),
65 fPidKn(0),
66 fTOFTPCsignal(0),
67 fCombsignal(0),
68 fCutsDaughter(NULL),
69 fPIDCombined(NULL),
70 fContPid(NULL),
71 fContPid2(NULL),
72 fHmismTOF(0),
9c668341 73 fHchannelTOFdistr(0),
74 fTypeCol(2)
a8ad4709 75{
76 // Default constructor (should not be used)
77 fList->SetName("contPhiBayes1");
78 fList2->SetName("contPhiBayes2");
79 fList3->SetName("contPhiBayes3");
80
81 fList->SetOwner(kTRUE);
82 fList2->SetOwner(kTRUE);
83 fList3->SetOwner(kTRUE);
84
85 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
86 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
87
88 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
89 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
90}
91
92//______________________________________________________________________________
93AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes(const char *name):
94 AliAnalysisTaskSE(name),
95 fVtxCut(10.0), // cut on |vertex| < fVtxCut
96 fEtaCut(0.8), // cut on |eta| < fEtaCut
97 fMinPt(0.15), // cut on pt > fMinPt
98 fIsMC(kFALSE),
99 fQAsw(kFALSE),
100 fNcluster(70),
101 fFilterBit(1),
102 fList(new TList()),
103 fList2(new TList()),
104 fList3(new TList()),
105 fCentrality(-1),
106 fPsi(0),
107 fPtPhi(0.),
108 fPhiPhi(0.),
109 fEtaPhi(0.),
110 fMassV0(0.),
111 fPtKp(0.),
112 fPhiKp(0.),
113 fEtaKp(0.),
114 fPtKn(0.),
115 fPhiKn(0.),
116 fEtaKn(0.),
117 fPidKp(0),
118 fPidKn(0),
119 fTOFTPCsignal(0),
120 fCombsignal(0),
121 fCutsDaughter(NULL),
122 fPIDCombined(NULL),
123 fContPid(NULL),
124 fContPid2(NULL),
125 fHmismTOF(0),
9c668341 126 fHchannelTOFdistr(0),
127 fTypeCol(2)
a8ad4709 128{
129
130 DefineOutput(1, TList::Class());
131 DefineOutput(2, TList::Class());
132 DefineOutput(3, TList::Class());
133 DefineOutput(4, TList::Class());
134
135 // Output slot #1 writes into a TTree
136 fList->SetName("contPhiBayes1");
137 fList2->SetName("contPhiBayes2");
138 fList3->SetName("contPhiBayes3");
139
140 fList->SetOwner(kTRUE);
141 fList2->SetOwner(kTRUE);
142 fList3->SetOwner(kTRUE);
143
144 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
145 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
146
147 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
148 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
149}
150//_____________________________________________________________________________
151AliAnalysisTaskPhiBayes::~AliAnalysisTaskPhiBayes()
152{
153
154}
155
156//______________________________________________________________________________
157void AliAnalysisTaskPhiBayes::UserCreateOutputObjects()
158{
159
160 fPIDCombined=new AliPIDCombined;
161 fPIDCombined->SetDefaultTPCPriors();
162 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
163
164 Float_t invMmin = 0.985;
165 Float_t invMmax = 1.045;
166
167 const Int_t nBinPid = 14;
168
169 Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaKP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
170
171 Int_t binPid2[nBinPid] = {1/*ptPhi*/,8/*EtaKN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
172
173 fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
174 fContPid->SetTitleX("M_{#phi}");
175 fContPid->SetTitleY("centrality (%)");
176 fContPid->SetVarName(0,"p_{T}^#phi}");
177 fContPid->SetVarRange(0,0,10);
178 fContPid->SetVarName(1,"#eta^{#phi}");
179 fContPid->SetVarRange(1,-0.8,0.8);
180 fContPid->SetVarName(2,"p_{T}^{Kp}");
181 fContPid->SetVarRange(2,0.3,4.3);
182 fContPid->SetVarName(3,"p_{T}^{Kn}");
183 fContPid->SetVarRange(3,0.3,4.3);
184 fContPid->SetVarName(4,"BayesProb^{Kp}");
185 fContPid->SetVarRange(4,0,1.);
186 fContPid->SetVarName(5,"BayesProb^{Kn}");
187 fContPid->SetVarRange(5,0,1.);
188 fContPid->SetVarName(6,"isTOFmatch^{Kp}");
189 fContPid->SetVarRange(6,-0.5,1.5);
190 fContPid->SetVarName(7,"isTOFmatch^{Kn}");
191 fContPid->SetVarRange(7,-0.5,1.5);
192 fContPid->SetVarName(8,"isPhiTrue^{Kn}");
193 fContPid->SetVarRange(8,-0.5,1.5);
194 fContPid->SetVarName(9,"N#sigma^{Kp}");
195 fContPid->SetVarRange(9,1.25,6.25);
196 fContPid->SetVarName(10,"N#sigma^{Kn}");
197 fContPid->SetVarRange(10,1.25,6.25);
198 fContPid->SetVarName(11,"#Delta#phi^{Kp}");
199 fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
200 fContPid->SetVarName(12,"#Delta#phi^{Kn}");
201 fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
202 fContPid->SetVarName(13,"#Psi");
203 fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
204
205 fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
206 fContPid2->SetTitleX("M_{#phi}");
207 fContPid2->SetTitleY("centrality (%)");
208 fContPid2->SetVarName(0,"p_{T}^{#phi}");
209 fContPid2->SetVarRange(0,0,10);
210 fContPid2->SetVarName(1,"#eta^{#phi}");
211 fContPid2->SetVarRange(1,-0.8,0.8);
212 fContPid2->SetVarName(2,"p_{T}^{Kp}");
213 fContPid2->SetVarRange(2,0.3,4.3);
214 fContPid2->SetVarName(3,"p_{T}^{Kn}");
215 fContPid2->SetVarRange(3,0.3,4.3);
216 fContPid2->SetVarName(4,"BayesProb^{Kp}");
217 fContPid2->SetVarRange(4,0,1.);
218 fContPid2->SetVarName(5,"BayesProb^{Kn}");
219 fContPid2->SetVarRange(5,0,1.);
220 fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
221 fContPid2->SetVarRange(6,-0.5,1.5);
222 fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
223 fContPid2->SetVarRange(7,-0.5,1.5);
224 fContPid2->SetVarName(8,"isPhiTrue^{Kn}");
225 fContPid2->SetVarRange(8,-0.5,1.5);
226 fContPid2->SetVarName(9,"N#sigma^{Kp}");
227 fContPid2->SetVarRange(9,1.25,6.25);
228 fContPid2->SetVarName(10,"N#sigma^{Kn}");
229 fContPid2->SetVarRange(10,1.25,6.25);
230 fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
231 fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
232 fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
233 fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
234 fContPid2->SetVarName(13,"#Psi");
235 fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
236
237
238
239 const Int_t nDETsignal = 100; // mass
240 Double_t binDETsignal[nDETsignal+1];
241 for(Int_t i=0;i<nDETsignal+1;i++){
242 binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
243 }
244 const Int_t nDETsignal2 = 10; // centrality
245 Double_t binDETsignal2[nDETsignal2+1];
246 for(Int_t i=0;i<nDETsignal2+1;i++){
247 binDETsignal2[i] = i*100./ nDETsignal2;
248 }
249 fContPid->AddSpecies("phi",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
250 fContPid2->AddSpecies("phi2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
251
252 fList->Add(fContPid);
253 fList->Add(fContPid2);
254
255 hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
256 hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
257 hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
258 hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
259
260 hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
261 hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
262 hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
263 hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
264
265 fList2->Add(hMatching[0]);
266 fList2->Add(hMatching[1]);
267 fList2->Add(hMatching[2]);
268 fList2->Add(hMatching[3]);
269 fList2->Add(hTracking[0]);
270 fList2->Add(hTracking[1]);
271 fList2->Add(hTracking[2]);
272 fList2->Add(hTracking[3]);
273
274 fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
275 fList2->Add(fTOFTPCsignal);
276 fCombsignal = new TH1F("hCombsignal","Combined signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
277 fList2->Add(fCombsignal);
278
279 // QA plots
280 char name[100],title[400];
281 for(Int_t i=0;i < nCentrBin;i++){
282 snprintf(name,100,"hPiTPCQA_%i",i);
283 snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
284 fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
285 fList3->Add(fPiTPC[i]);
286
287 snprintf(name,100,"hKaTPCQA_%i",i);
288 snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
289 fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
290 fList3->Add(fKaTPC[i]);
291
292 snprintf(name,100,"hPrTPCQA_%i",i);
293 snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
294 fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
295 fList3->Add(fPrTPC[i]);
296
297 snprintf(name,100,"hElTPCQA_%i",i);
298 snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
299 fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
300 fList3->Add(fElTPC[i]);
301
302 snprintf(name,100,"hPiTOFQA_%i",i);
303 snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
304 fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
305 fList3->Add(fPiTOF[i]);
306
307 snprintf(name,100,"hKaTOFQA_%i",i);
308 snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
309 fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
310 fList3->Add(fKaTOF[i]);
311
312 snprintf(name,100,"hPrTOFQA_%i",i);
313 snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
314 fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
315 fList3->Add(fPrTOF[i]);
316
317 snprintf(name,100,"hElTOFQA_%i",i);
318 snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
319 fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
320 fList3->Add(fElTOF[i]);
321
322 }
323
324 // Post output data.
325 PostData(1, fList);
326 PostData(2, fList2);
327 PostData(3, fList3);
328}
329
330//______________________________________________________________________________
331void AliAnalysisTaskPhiBayes::UserExec(Option_t *)
332{
333 // Main loop
334 // Called for each event
335
336 fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
337 if(!fOutputAOD){
338 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
339 this->Dump();
340 return;
341 }
342
343 Float_t zvtx = GetVertex(fOutputAOD);
344
345
346
347 //Get the MC object
348 if(fIsMC){
349 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
350 if (!mcHeader) {
351 AliError("Could not find MC Header in AOD");
352 return;
353 }
354 }
355
356 if (TMath::Abs(zvtx) < fVtxCut) {
357
358 //Centrality
359 Float_t v0Centr = -10.;
360 Float_t trkCentr = -10.;
361 AliCentrality *centrality = fOutputAOD->GetCentrality();
362 if (centrality){
363 trkCentr = centrality->GetCentralityPercentile("V0M");
364 v0Centr = centrality->GetCentralityPercentile("TRK");
365 }
366
9c668341 367 if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
a8ad4709 368 fCentrality = v0Centr;
369 Analyze(fOutputAOD); // Do analysis!!!!
370
371 }
372 }
373
374}
375
376//________________________________________________________________________
377void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
378{
379 Int_t ntrack = aodEvent->GetNumberOfTracks();
380
381 fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
382 fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
383 fPidKp=0,fPidKn=0;
384 fMassV0=-1;
385
386 TLorentzVector phiV;
387
388 Double_t px,py,pz,E;
389
390 Float_t invMmin = 0.985;
391 Float_t invMmax = 1.045;
392
393 Int_t icentr = 8;
394 if(fCentrality < 0) icentr = 8;
395 else if(fCentrality < 10) icentr = 0;
396 else if(fCentrality < 20) icentr = 1;
397 else if(fCentrality < 30) icentr = 2;
398 else if(fCentrality < 40) icentr = 3;
399 else if(fCentrality < 50) icentr = 4;
400 else if(fCentrality < 60) icentr = 5;
401 else if(fCentrality < 70) icentr = 6;
402 else if(fCentrality < 80) icentr = 7;
403
404 Float_t addMismatchForMC = 0.005;
405 if(fCentrality < 50) addMismatchForMC += 0.01;
406 if(fCentrality < 20) addMismatchForMC += 0.02;
407
9c668341 408 if(fTypeCol == 0) addMismatchForMC = 0.005;
409 else if(fTypeCol == 1) addMismatchForMC *= 0.5;
410
a8ad4709 411 fPsi = 0;
412 /* Compute TPC EP */
413 Double_t Qx2 = 0, Qy2 = 0;
414 Double_t Qx3 = 0, Qy3 = 0;
415 for(Int_t iT = 0; iT < ntrack; iT++) {
416 AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
417
418 if (!aodTrack){
419 continue;
420 }
421
422 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
423
424 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
425 continue;
426
427 Double_t b[2] = {-99., -99.};
428 Double_t bCov[3] = {-99., -99., -99.};
429 if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
430 continue;
431
432 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
433 continue;
434
435 Qx2 += TMath::Cos(2*aodTrack->Phi());
436 Qy2 += TMath::Sin(2*aodTrack->Phi());
437 Qx3 += TMath::Cos(3*aodTrack->Phi());
438 Qy3 += TMath::Sin(3*aodTrack->Phi());
439
440 }
441
442 fPsi = TMath::ATan2(Qy2, Qx2)/2.;
443
444 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
445 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
446 AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
447
448 PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
449 PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
450 PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
451 PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
452
453 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
454
455 Double_t probP[10],probN[10];
456 Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
457
458
459 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
460 TClonesArray *mcArray = NULL;
461 if (mcHeader)
462 mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
463
464 Int_t nmc = 0;
465 if(mcArray)
466 nmc = mcArray->GetEntries();
467
468 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
469 AliAODTrack* track = aodEvent->GetTrack(i);
470
471 AliAODMCParticle *mcp = NULL;
472 Int_t pdg = 0;
473
474 if (!track){
475 continue;
476 }
477
478 Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
479
480 Int_t label = -1;
481 if(mcArray){
482 label = track->GetLabel();
483 if(label != -1 && label < nmc){
484 label = TMath::Abs(label);
485 mcp = (AliAODMCParticle*)mcArray->At(label);
486 pdg = TMath::Abs(mcp->GetPdgCode());
487 }
488 else
489 label = -1;
490 }
491 else{
492 /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
493 }
494
495 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
496 hTracking[0]->Fill(track->Pt(),fCentrality);
497 if(pdg == 211)
498 hTracking[1]->Fill(track->Pt(),fCentrality);
499 else if(pdg == 321)
500 hTracking[2]->Fill(track->Pt(),fCentrality);
501 else if(pdg == 2212)
502 hTracking[3]->Fill(track->Pt(),fCentrality);
503 else if(! mcp){ // Fill matching histo with the prob
504 hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
505 hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
506 hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
507 }
508 }
509
510 if(!tofMatch) continue;
511
512 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
513 hMatching[0]->Fill(track->Pt(),fCentrality);
514 if(pdg == 211)
515 hMatching[1]->Fill(track->Pt(),fCentrality);
516 else if(pdg == 321)
517 hMatching[2]->Fill(track->Pt(),fCentrality);
518 else if(pdg == 2212)
519 hMatching[3]->Fill(track->Pt(),fCentrality);
520 else if(! mcp){ // Fill matching histo with the prob
521 hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
522 hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
523 hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
524 }
525 }
526 }
527
007df888 528// Int_t pdg1 = -1;
529// Int_t pdg2 = -1;
a8ad4709 530
531
532 // start analysis phi
533 for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
534 AliAODTrack* KpTrack = aodEvent->GetTrack(i);
535
536 if (!KpTrack){
537 continue;
538 }
539
540 if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3 && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
541
e34b28fe 542 nSigmaTOF=5;
a8ad4709 543 nSigmaComb=5;
544 fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
545 fPidKp=0;
546
547 UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
548
549 Double_t oldpP[10];
550 fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
551
552 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
553 fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
554 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
555 fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
556 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
557 fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
558 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
559 fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
560
561 if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
562
563 Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
9c668341 564 /*
a8ad4709 565 if(mcArray){
566 Int_t labelK = TMath::Abs(KpTrack->GetLabel());
567 AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 568 pdg1 = TMath::Abs(mcp1->GetPdgCode());
a8ad4709 569 }
9c668341 570 */
a8ad4709 571
572 fPidKp = Int_t(probP[3]*100);
573
574 if(tofMatch1){
575 if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
576 // remove this tof hit
577 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
578 detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
579 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
580 fPidKp = Int_t(probP[4]*100);
581 tofMatch1=0;
582 }
583 else{
584 if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
585
586 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
587 fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
588 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
589 fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
590 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
591 fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
592 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
593 fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
594
595 nSigmaTOF = TMath::Abs(nSigmaTOF);
596
597 if(fIsMC){
598 Float_t mismAdd = addMismatchForMC;
599 if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
600
601 if(gRandom->Rndm() < mismAdd){
602 Float_t etaAbs = TMath::Abs(KpTrack->Eta());
603 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
604 channel = channel % 8736;
605 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
606
607 // generate random time
608 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
609 Double_t times[10];
610 KpTrack->GetIntegratedTimes(times);
611 nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
612 }
613 }
614
615 if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
616
617 if(nSigmaTOF < 2) fPidKp += 256;
618 else if(nSigmaTOF < 3) fPidKp += 512;
619 }
620 }
621
622 if(tofMatch1){
623 nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
624 if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
625 fCombsignal->Fill(nSigmaComb);
626 }
627 } else {
628 nSigmaComb = TMath::Abs(nSigmaTPC);
629 }
630
631 // use sigmaTOF instead of sigmaComb
e34b28fe 632 if(tofMatch1) nSigmaComb = nSigmaTOF;
a8ad4709 633
634 if(nSigmaComb < 2) nSigmaComb = 2;
635 else if(nSigmaComb < 3) nSigmaComb = 3;
636 else if(nSigmaComb < 5) nSigmaComb = 4.99;
637 else nSigmaComb = 6;
638
639 if(fPtKp > 4.299) fPtKp = 4.299;
640
641 for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
642 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
643
644 if (!KnTrack){
645 continue;
646 }
647
648 if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
649
650 px = KpTrack->Px() + KnTrack->Px();
651 py = KpTrack->Py() + KnTrack->Py();
652 pz = KpTrack->Pz() + KnTrack->Pz();
653 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
654 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
655
656 phiV.SetPxPyPzE(px,py,pz,E);
657 fMassV0 = phiV.M();
658
659 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
660
661 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
662 fPidKn=0;
663
664 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
665 Double_t oldpN[10];
666 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
667
668 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
669
670 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
671
e34b28fe 672 nSigmaTOF2=5;
a8ad4709 673 nSigmaComb2=5;
674
675 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 676 /*
a8ad4709 677 if(mcArray){
678 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
679 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 680 pdg2 = TMath::Abs(mcp2->GetPdgCode());
681 }
682 */
a8ad4709 683
684 fPidKn = Int_t(probN[3]*100);
685
686 if(tofMatch2){
687 if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
688 // remove this tof hit
689 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
690 detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
691 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
692 fPidKn = Int_t(probN[4]*100);
693 tofMatch2=0;
694 }
695 else{
696 if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
697
698 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
699
700 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
701
702 if(fIsMC){
703 Float_t mismAdd = addMismatchForMC;
704 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
705
706 if(gRandom->Rndm() < mismAdd){
707 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
708 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
709 channel = channel % 8736;
710 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
711
712 // generate random time
713 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
714 Double_t times[10];
715 KnTrack->GetIntegratedTimes(times);
e34b28fe 716 nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
a8ad4709 717 }
718 }
719
720 if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
721
722 if(nSigmaTOF2 < 2) fPidKn += 256;
723 else if(nSigmaTOF2 < 3) fPidKn += 512;
724 }
725 }
726
727 fPtPhi = phiV.Pt();
728 fEtaPhi = phiV.Eta();
729 fPhiPhi = phiV.Phi();
730
731 if(tofMatch2){
732 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
733 if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
734 fCombsignal->Fill(nSigmaComb2);
735 }
736 } else {
737 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
738 }
739
740 // use sigmaTOF instead of sigmaComb
e34b28fe 741 if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
a8ad4709 742
743 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
744 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
745 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
746 else nSigmaComb2 = 6;
747
748 Bool_t isTrue = kFALSE;
749
750 if(mcArray){
751 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
752 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
753
754 if(labelP > -1 && labelN > -1){
755 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
756 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
757
758 Int_t mP = partP->GetMother();
759 Int_t mN = partN->GetMother();
760 if(mP == mN && mP > -1){
761 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
762 Int_t pdgM = partM->GetPdgCode();
763 if(pdgM == 333) isTrue = kTRUE;
764 }
765 }
766
767 }
768
769 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
770 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
771
772 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
773 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
774 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
775 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
776
777 if(fPtKn > 4.299) fPtKn = 4.299;
778
779 Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
780
781
782 Int_t ipt = 0;
783 while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
784 ipt++;
785 }
786 ipt--;
787 if(ipt < 0) ipt = 0; // just to be sure
788
789 if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
790 if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
791 xTOfill[1] = KnTrack->Eta();
792 if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);// use tagging on positive track
793
794 }
795
796
797 }
798 } // end analysi phi
799
800}
801
802//_____________________________________________________________________________
803Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
804{
805
806 Float_t zvtx = -999;
807
808 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
809 if (!vtxAOD)
810 return zvtx;
811 if(vtxAOD->GetNContributors()>0)
812 zvtx = vtxAOD->GetZ();
813
814 return zvtx;
815}
816//_____________________________________________________________________________
817void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
818{
819 // Terminate loop
820 Printf("Terminate()");
821}
822//-------------------------------------------------------------------------
823Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
824 if(!fIsMC) return 1; // used only on MC
825
9c668341 826 if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 827 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
828
829 if(!(channel%20)) return 0; // 5% additional loss in MC
830 }
831
832 return 1;
833}