]>
Commit | Line | Data |
---|---|---|
6279d59e | 1 | // |
2 | // Implementation file for implementation of data analysis aft 900 GeV | |
3 | // | |
4 | // Author: A. Pulvirenti | |
5 | // | |
6 | ||
7 | #include "Riostream.h" | |
4c06ce33 | 8 | #include <iomanip> |
6279d59e | 9 | |
10 | #include "TH1.h" | |
5faf5a07 | 11 | #include "TH3.h" |
12 | #include "TFile.h" | |
6279d59e | 13 | #include "TTree.h" |
14 | #include "TParticle.h" | |
15 | #include "TRandom.h" | |
16 | #include "TLorentzVector.h" | |
7356f978 | 17 | #include "TDatabasePDG.h" |
6279d59e | 18 | |
19 | #include "AliLog.h" | |
20 | #include "AliESDpid.h" | |
21 | #include "AliESDEvent.h" | |
22 | #include "AliESDVertex.h" | |
23 | #include "AliESDtrack.h" | |
5faf5a07 | 24 | #include "AliESDtrackCuts.h" |
6279d59e | 25 | #include "AliStack.h" |
26 | #include "AliMCEvent.h" | |
27 | #include "AliTOFT0maker.h" | |
28 | #include "AliTOFcalib.h" | |
29 | #include "AliCDBManager.h" | |
4c06ce33 | 30 | #include "AliITSPIDResponse.h" |
7356f978 | 31 | #include "AliCFContainer.h" |
32 | #include "AliVEventHandler.h" | |
33 | #include "AliESDInputHandler.h" | |
34 | #include "AliMultiInputEventHandler.h" | |
35 | #include "AliAnalysisManager.h" | |
6279d59e | 36 | |
5faf5a07 | 37 | #include "AliRsnEvent.h" |
38 | #include "AliRsnTarget.h" | |
39 | #include "AliRsnDaughter.h" | |
5faf5a07 | 40 | |
6279d59e | 41 | #include "AliRsnAnalysisPhi7TeV.h" |
42 | ||
43 | //__________________________________________________________________________________________________ | |
5faf5a07 | 44 | AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name, Bool_t isMC) : |
7356f978 | 45 | AliAnalysisTaskSE(name), |
46 | fIsMC(isMC), | |
47 | fAddTPC(kTRUE), | |
48 | fAddITS(kFALSE), | |
49 | fCheckITS(kTRUE), | |
50 | fCheckTPC(kTRUE), | |
51 | fCheckTOF(kTRUE), | |
52 | fStep(1000), | |
53 | fPhiMass(1.019455), | |
54 | fKaonMass(0.49677), | |
55 | fMaxVz(10.0), | |
56 | fITSNSigma(3.0), | |
57 | fTPCMomentumThreshold(0.350), | |
58 | fTOFNSigma(3.0), | |
59 | fESDtrackCutsTPC(0), | |
60 | fESDtrackCutsITS(0), | |
61 | fESDpid(0), | |
62 | fOutList(0x0), | |
63 | fCFunlike(0), | |
64 | fCFlikePP(0), | |
65 | fCFlikeMM(0), | |
66 | fCFtrues(0), | |
67 | fCFkaons(0), | |
68 | fHEvents(0x0) | |
6279d59e | 69 | { |
70 | // | |
71 | // Constructor | |
72 | // | |
73 | ||
7356f978 | 74 | fTPCNSigma[0] = 5.0; |
75 | fTPCNSigma[1] = 3.0; | |
76 | ||
77 | fVertexX[0] = fVertexX[1] = 0; | |
78 | fVertexY[0] = fVertexZ[1] = 0; | |
79 | fVertexZ[0] = fVertexZ[1] = 0; | |
80 | ||
81 | DefineOutput(1, TList::Class()); | |
6279d59e | 82 | } |
83 | ||
84 | //__________________________________________________________________________________________________ | |
85 | AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) : | |
7356f978 | 86 | AliAnalysisTaskSE(copy), |
87 | fIsMC(copy.fIsMC), | |
88 | fAddTPC(copy.fAddTPC), | |
89 | fAddITS(copy.fAddITS), | |
90 | fCheckITS(copy.fCheckITS), | |
91 | fCheckTPC(copy.fCheckTPC), | |
92 | fCheckTOF(copy.fCheckTOF), | |
93 | fStep(copy.fStep), | |
94 | fPhiMass(copy.fPhiMass), | |
95 | fKaonMass(copy.fKaonMass), | |
96 | fMaxVz(copy.fMaxVz), | |
97 | fITSNSigma(copy.fITSNSigma), | |
98 | fTPCMomentumThreshold(copy.fTPCMomentumThreshold), | |
99 | fTOFNSigma(copy.fTOFNSigma), | |
100 | fESDtrackCutsTPC(copy.fESDtrackCutsTPC), | |
101 | fESDtrackCutsITS(copy.fESDtrackCutsITS), | |
102 | fESDpid(copy.fESDpid), | |
103 | fOutList(0x0), | |
104 | fCFunlike(0), | |
105 | fCFlikePP(0), | |
106 | fCFlikeMM(0), | |
107 | fCFtrues(0), | |
108 | fCFkaons(0), | |
109 | fHEvents(0x0) | |
6279d59e | 110 | { |
111 | // | |
112 | // Copy constructor | |
113 | // | |
5faf5a07 | 114 | |
7356f978 | 115 | fTPCNSigma[0] = copy.fTPCNSigma[0]; |
116 | fTPCNSigma[1] = copy.fTPCNSigma[1]; | |
117 | ||
118 | fVertexX[0] = fVertexX[1] = 0; | |
119 | fVertexY[0] = fVertexZ[1] = 0; | |
120 | fVertexZ[0] = fVertexZ[1] = 0; | |
6279d59e | 121 | } |
122 | ||
123 | //__________________________________________________________________________________________________ | |
124 | AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7TeV& copy) | |
125 | { | |
126 | // | |
127 | // Assignment operator | |
128 | // | |
aaf7af8d | 129 | if (this == ©) |
130 | return *this; | |
7356f978 | 131 | fIsMC = copy.fIsMC; |
132 | fAddTPC = copy.fAddTPC; | |
133 | fAddITS = copy.fAddITS; | |
134 | fCheckITS = copy.fCheckITS; | |
135 | fCheckTPC = copy.fCheckTPC; | |
136 | fCheckTOF = copy.fCheckTOF; | |
137 | fStep = copy.fStep; | |
138 | fPhiMass = copy.fPhiMass; | |
139 | fKaonMass = copy.fKaonMass; | |
140 | fMaxVz = copy.fMaxVz; | |
141 | fITSNSigma = copy.fITSNSigma; | |
142 | fTPCNSigma[0] = copy.fTPCNSigma[0]; | |
143 | fTPCNSigma[1] = copy.fTPCNSigma[1]; | |
144 | fTPCMomentumThreshold = copy.fTPCMomentumThreshold; | |
145 | fTOFNSigma = copy.fTOFNSigma; | |
146 | fESDtrackCutsTPC = copy.fESDtrackCutsTPC; | |
147 | fESDtrackCutsITS = copy.fESDtrackCutsITS; | |
148 | fESDpid = copy.fESDpid; | |
149 | ||
150 | return (*this); | |
6279d59e | 151 | } |
152 | ||
153 | //__________________________________________________________________________________________________ | |
154 | AliRsnAnalysisPhi7TeV::~AliRsnAnalysisPhi7TeV() | |
155 | { | |
156 | // | |
157 | // Destructor | |
5faf5a07 | 158 | // |
159 | } | |
160 | ||
6279d59e | 161 | //__________________________________________________________________________________________________ |
162 | void AliRsnAnalysisPhi7TeV::UserCreateOutputObjects() | |
163 | { | |
164 | // | |
165 | // Create the output data container | |
166 | // | |
a3dd4b4d | 167 | |
7356f978 | 168 | // create output list |
169 | OpenFile(1); | |
170 | fOutList = new TList; | |
171 | ||
172 | // event related histograms | |
173 | fHEvents = new TH1I("hEvents", "Event details", kVertexTypes, 0, kVertexTypes); | |
174 | fVertexX[0] = new TH1F("hVertexTRKX", "X position of primary vertex (tracks)", 200, -2, 2); | |
175 | fVertexY[0] = new TH1F("hVertexTRKY", "Y position of primary vertex (tracks)", 200, -2, 2); | |
176 | fVertexZ[0] = new TH1F("hVertexTRKZ", "Z position of primary vertex (tracks)", 400, -20, 20); | |
177 | fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)" , 200, -2, 2); | |
178 | fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)" , 200, -2, 2); | |
179 | fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)" , 400, -20, 20); | |
180 | ||
181 | fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks"); | |
182 | fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD"); | |
183 | fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks"); | |
184 | fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD"); | |
185 | fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex"); | |
186 | fHEvents->GetXaxis()->SetBinLabel(6, "Empty event"); | |
187 | ||
188 | // define axes: | |
189 | // pairs [0] = inv. mass | |
190 | // [1] = inv. mass resolution | |
191 | // [2] = pt | |
192 | // [3] = rec. multiplicity | |
193 | // [4] = MC multiplicity (if available) | |
194 | // | |
195 | // kaons [0] = pt | |
196 | // [1] = pTPC | |
197 | // [2] = DCAr | |
198 | // [3] = DCAz | |
199 | // [4] = signalITS | |
200 | // [5] = signalTPC | |
201 | // [6] = betaTOF | |
202 | // [7] = rec. multiplicity | |
203 | // [8] = MC multiplicity (if available) | |
204 | Double_t minIM = 0.9, maxIM = 1.4; | |
205 | Double_t minIMRes = -5.0, maxIMRes = 5.0; | |
206 | Double_t minMom = 0.0, maxMom = 5.0; | |
207 | Double_t minDCAr = 0.0, maxDCAr = 1.0; | |
208 | Double_t minDCAz = 0.0, maxDCAz = 5.0; | |
209 | Double_t minITS = 0.0, maxITS = 800.0; // raw signal | |
210 | Double_t minTPC = 0.0, maxTPC = 800.0; // raw signal | |
211 | Double_t minTOF = 0.0, maxTOF = 1.0; // beta = length / TOF time / c | |
212 | Double_t mult[] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 500.}; | |
213 | Int_t nIM = 500; | |
214 | Int_t nIMRes = 100; | |
215 | Int_t nMom = 50; | |
216 | Int_t nDCAr = 100; | |
217 | Int_t nDCAz = 100; | |
218 | Int_t nITS = 800; | |
219 | Int_t nTPC = 800; | |
220 | Int_t nTOF = 100; | |
221 | Int_t nMult = sizeof(mult) / sizeof(mult[0]) - 1; | |
222 | Int_t nBins1[9] = {nMom, nMom, nDCAr, nDCAz, nITS, nTPC, nTOF, nMult, nMult}; | |
223 | Int_t nBins2[5] = {nIM , nIMRes, nMom, nMult, nMult}; | |
224 | ||
225 | // containers for analysis pairs have 2 steps (quality, quality+PID) | |
226 | // and those for efficiency have 3 (add full MC) | |
227 | // REMINDER: to use variable bins, call container->SetBinLimits(iaxis, Double_t *array); | |
228 | fCFunlike = new AliCFContainer("CFUnlike", "", 2, 5, nBins2); | |
229 | fCFlikePP = new AliCFContainer("CFLikePP", "", 2, 5, nBins2); | |
230 | fCFlikeMM = new AliCFContainer("CFLikeMM", "", 2, 5, nBins2); | |
231 | fCFtrues = new AliCFContainer("CFtrues" , "", 3, 5, nBins2); | |
232 | fCFkaons = new AliCFContainer("CFkaons" , "", 3, 9, nBins1); | |
233 | ||
234 | // create a unique temp array to all pair containers | |
235 | AliCFContainer *CFpair[4]; | |
236 | CFpair[0] = fCFunlike; | |
237 | CFpair[1] = fCFlikePP; | |
238 | CFpair[2] = fCFlikeMM; | |
239 | CFpair[3] = fCFtrues; | |
240 | ||
241 | // define axes: | |
242 | // pairs [0] = inv. mass | |
243 | // [1] = inv. mass resolution | |
244 | // [2] = pt | |
245 | // [3] = rec. multiplicity | |
246 | // [4] = MC multiplicity (if available) | |
247 | // | |
248 | // kaons [0] = pt | |
249 | // [1] = pTPC | |
250 | // [2] = DCAr | |
251 | // [3] = DCAz | |
252 | // [4] = signalITS | |
253 | // [5] = signalTPC | |
254 | // [6] = betaTOF | |
255 | // [7] = rec. multiplicity | |
256 | // [8] = MC multiplicity (if available) | |
257 | Int_t i; | |
258 | for (i = 0; i < 4; i++) { | |
259 | CFpair[i]->SetBinLimits(0, minIM , maxIM ); | |
260 | CFpair[i]->SetBinLimits(1, minIMRes, maxIMRes); | |
261 | CFpair[i]->SetBinLimits(2, minMom , maxMom ); | |
262 | CFpair[i]->SetBinLimits(3, mult); | |
263 | CFpair[i]->SetBinLimits(4, mult); | |
264 | } | |
265 | fCFkaons->SetBinLimits(0, minMom , maxMom ); | |
266 | fCFkaons->SetBinLimits(1, minMom , maxMom ); | |
267 | fCFkaons->SetBinLimits(2, minDCAr, maxDCAr); | |
268 | fCFkaons->SetBinLimits(3, minDCAz, maxDCAz); | |
269 | fCFkaons->SetBinLimits(4, minITS , maxITS ); | |
270 | fCFkaons->SetBinLimits(5, minTPC , maxTPC ); | |
271 | fCFkaons->SetBinLimits(6, minTOF , maxTOF ); | |
272 | fCFkaons->SetBinLimits(7, mult); | |
273 | fCFkaons->SetBinLimits(8, mult); | |
274 | ||
275 | // add all outputs to the list | |
276 | fOutList->Add(fHEvents); | |
277 | fOutList->Add(fVertexX[0]); | |
278 | fOutList->Add(fVertexY[0]); | |
279 | fOutList->Add(fVertexZ[0]); | |
280 | fOutList->Add(fVertexX[1]); | |
281 | fOutList->Add(fVertexY[1]); | |
282 | fOutList->Add(fVertexZ[1]); | |
283 | fOutList->Add(fCFunlike); | |
284 | fOutList->Add(fCFlikePP); | |
285 | fOutList->Add(fCFlikeMM); | |
286 | fOutList->Add(fCFtrues); | |
287 | fOutList->Add(fCFkaons); | |
288 | ||
289 | // update histogram container | |
290 | PostData(1, fOutList); | |
6279d59e | 291 | } |
292 | ||
293 | //__________________________________________________________________________________________________ | |
294 | void AliRsnAnalysisPhi7TeV::UserExec(Option_t *) | |
295 | { | |
296 | // | |
297 | // Main execution function. | |
298 | // Fills the fHEvents data member with the following legenda: | |
299 | // 0 -- event OK, prim vertex with tracks | |
300 | // 1 -- event OK, prim vertex with SPD | |
301 | // 2 -- event OK but vz large | |
302 | // 3 -- event bad | |
303 | // | |
304 | ||
7356f978 | 305 | static Int_t evNum = 0; |
306 | evNum++; | |
307 | if (evNum % fStep == 0) AliInfo(Form("Processing event #%7d", evNum)); | |
308 | ||
309 | // retrieve ESD event and related stack (if available) | |
310 | AliESDEvent *esd = static_cast<AliESDEvent*>(fInputEvent); | |
311 | AliMCEvent *mc = static_cast<AliMCEvent*>(fMCEvent); | |
312 | ||
313 | // if the AliESDpid object is not assigned, do this now | |
314 | // find an AliESDpid from input handler | |
315 | if (!fESDpid) { | |
316 | AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); | |
317 | AliVEventHandler *evh = mgr->GetInputEventHandler(); | |
318 | AliESDInputHandler *esdin = 0x0; | |
319 | if (evh->IsA() == AliMultiInputEventHandler::Class()) { | |
320 | AliMultiInputEventHandler *multh = static_cast<AliMultiInputEventHandler*>(evh); | |
321 | AliInputEventHandler *input = multh->GetFirstInputEventHandler(); | |
322 | if (input->IsA() == AliESDInputHandler::Class()) { | |
323 | esdin = static_cast<AliESDInputHandler*>(input); | |
324 | } | |
325 | } | |
326 | else if (evh->IsA() == AliESDInputHandler::Class()) { | |
327 | esdin = static_cast<AliESDInputHandler*>(evh); | |
328 | } | |
329 | if (!esdin) { | |
330 | AliFatal("Required ESD input handler"); | |
331 | return; | |
332 | } | |
333 | fESDpid = esdin->GetESDpid(); | |
334 | } | |
a3dd4b4d | 335 | |
7356f978 | 336 | // check the event |
337 | EVertexType eval = EventEval(esd); | |
338 | fHEvents->Fill((Int_t)eval); | |
a3dd4b4d | 339 | |
7356f978 | 340 | // if the event is good for analysis, process it |
341 | if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex) { | |
342 | ProcessESD(esd, mc); | |
343 | if (mc) ProcessMC (esd, mc); | |
344 | } | |
a3dd4b4d | 345 | |
7356f978 | 346 | // update histogram container |
347 | PostData(1, fOutList); | |
6279d59e | 348 | } |
349 | ||
350 | //__________________________________________________________________________________________________ | |
351 | void AliRsnAnalysisPhi7TeV::Terminate(Option_t *) | |
352 | { | |
353 | // | |
354 | // Terminate | |
355 | // | |
356 | } | |
357 | ||
4c06ce33 | 358 | //__________________________________________________________________________________________________ |
7356f978 | 359 | AliRsnAnalysisPhi7TeV::EVertexType AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd) |
4c06ce33 | 360 | { |
361 | // | |
362 | // Checks if the event is good for analysis. | |
363 | // Returns one of the flag values defined in the header | |
364 | // | |
365 | ||
7356f978 | 366 | // reject empty events |
367 | Int_t ntracks = esd->GetNumberOfTracks(); | |
368 | if (!ntracks) return kEmptyEvent; | |
369 | ||
370 | // get the best primary vertex: | |
371 | // first try the one with tracks | |
372 | const AliESDVertex *vTrk = esd->GetPrimaryVertexTracks(); | |
373 | const AliESDVertex *vSPD = esd->GetPrimaryVertexSPD(); | |
374 | Double_t vzTrk = 1000000.0; | |
375 | Double_t vzSPD = 1000000.0; | |
376 | Int_t ncTrk = -1; | |
377 | Int_t ncSPD = -1; | |
378 | if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors(); | |
379 | if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors(); | |
380 | if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv()); | |
381 | if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv()); | |
382 | if (vTrk && ncTrk > 0) { | |
383 | // fill the histograms | |
384 | fVertexX[0]->Fill(vTrk->GetXv()); | |
385 | fVertexY[0]->Fill(vTrk->GetYv()); | |
386 | fVertexZ[0]->Fill(vTrk->GetZv()); | |
387 | ||
388 | // check VZ position | |
389 | if (vzTrk <= fMaxVz) | |
390 | return kGoodTracksPrimaryVertex; | |
391 | else | |
392 | return kFarTracksPrimaryVertex; | |
393 | } else if (vSPD && ncSPD > 0) { | |
394 | // fill the histograms | |
395 | fVertexX[1]->Fill(vSPD->GetXv()); | |
396 | fVertexY[1]->Fill(vSPD->GetYv()); | |
397 | fVertexZ[1]->Fill(vSPD->GetZv()); | |
398 | ||
399 | // check VZ position | |
400 | if (vzSPD <= fMaxVz) | |
401 | return kGoodSPDPrimaryVertex; | |
402 | else | |
403 | return kFarSPDPrimaryVertex; | |
404 | } else | |
405 | return kNoGoodPrimaryVertex; | |
4c06ce33 | 406 | } |
407 | ||
6279d59e | 408 | //__________________________________________________________________________________________________ |
409 | void AliRsnAnalysisPhi7TeV::ProcessESD | |
7356f978 | 410 | (AliESDEvent *esd, AliMCEvent *mc) |
6279d59e | 411 | { |
412 | // | |
413 | // This function works with the ESD object | |
7356f978 | 414 | // and fills all containers |
415 | // | |
416 | ||
417 | // some utility arrays, declared static | |
418 | // to avoid creting/destroying at each event | |
419 | static TArrayC quality; // will be 1 or 0 depending if cut is passed or not | |
420 | static TArrayC pid; // will be 1 or 0 depending if cut is passed or not | |
421 | ||
422 | // get stack, if possible | |
423 | AliStack *stack = 0x0; | |
424 | if (mc) stack = mc->Stack(); | |
425 | ||
426 | // get number of tracks, which is used for setting arrays | |
427 | Int_t itrack, ntracks = esd->GetNumberOfTracks(); | |
428 | quality.Set(ntracks); | |
429 | pid .Set(ntracks); | |
430 | ||
431 | // loop on tracks and sets the flags | |
432 | for (itrack = 0; itrack < ntracks; itrack++) { | |
433 | ||
434 | // get track | |
435 | AliESDtrack *track = esd->GetTrack(itrack); | |
436 | ||
437 | // check cuts | |
438 | quality[itrack] = (Char_t)OkQuality(track); | |
439 | pid [itrack] = (Char_t)OkPID(track); | |
440 | } | |
441 | ||
442 | // now, we can fill the histograms | |
443 | Int_t i1, i2, label; | |
444 | Double_t var[5]; | |
445 | TLorentzVector p1[2], p2[2], sum[2]; | |
446 | AliESDtrack *tr1, *tr2; | |
447 | TParticle *pr1, *pr2; | |
448 | ETrackType type1, type2; | |
449 | ||
450 | // multiplicity is defined for whole event | |
451 | var[3] = AliESDtrackCuts::GetReferenceMultiplicity(esd, kTRUE); | |
452 | var[4] = (mc ? mc->GetNumberOfTracks() : 0.0); | |
453 | ||
454 | for (i1 = 0; i1 < ntracks; i1++) { | |
455 | ||
456 | // check only quality | |
457 | if (!quality[i1]) continue; | |
458 | ||
459 | // skip tracks which are not of a good type | |
460 | tr1 = esd->GetTrack(i1); | |
461 | type1 = TrackType(tr1); | |
462 | if (type1 >= kTrackTypes) continue; | |
463 | ||
464 | // skip neutral | |
465 | tr1 = esd->GetTrack(i1); | |
466 | if (tr1->Charge() == 0) continue; | |
467 | ||
468 | // try to retrieve the MC | |
469 | pr1 = 0x0; | |
470 | label = TMath::Abs(tr1->GetLabel()); | |
471 | if (stack) { | |
472 | if (label >= 0 && label < stack->GetNtrack()) { | |
473 | pr1 = stack->Particle(label); | |
474 | } | |
0d73200d | 475 | } |
7356f978 | 476 | |
477 | // initialize 4-momenta | |
478 | p1[0].SetXYZM(tr1->Px(), tr1->Py(), tr1->Pz(), fKaonMass); | |
479 | if (pr1) p1[1].SetXYZM(pr1->Px(), pr1->Py(), pr1->Pz(), fKaonMass); | |
480 | ||
481 | for (i2 = i1+1; i2 < ntracks; i2++) { | |
482 | ||
483 | // check only quality | |
484 | if (!quality[i2]) continue; | |
485 | ||
486 | // skip tracks which are not of a good type | |
487 | tr2 = esd->GetTrack(i2); | |
488 | type2 = TrackType(tr2); | |
489 | if (type2 >= kTrackTypes) continue; | |
490 | ||
491 | // skip neutrals | |
492 | tr2 = esd->GetTrack(i2); | |
493 | if (tr2->Charge() == 0) continue; | |
494 | ||
495 | // try to retrieve MC | |
496 | pr2 = 0x0; | |
497 | label = TMath::Abs(tr2->GetLabel()); | |
498 | if (stack) { | |
499 | if (label >= 0 && label < stack->GetNtrack()) { | |
500 | pr2 = stack->Particle(label); | |
501 | } | |
502 | } | |
503 | ||
504 | // initialize 4-momenta | |
505 | p2[0].SetXYZM(tr2->Px(), tr2->Py(), tr2->Pz(), fKaonMass); | |
506 | if (pr2) p1[1].SetXYZM(pr2->Px(), pr2->Py(), pr2->Pz(), fKaonMass); | |
507 | ||
508 | // sum 4-momenta | |
509 | for (Int_t i = 0; i < 2; i++) sum[i] = p1[i] + p2[i]; | |
510 | ||
511 | // compute values | |
512 | // [0] = inv. mass | |
513 | // [1] = inv. mass resolution | |
514 | // [2] = pt | |
515 | // [3] = rec. multiplicity | |
516 | // [4] = MC multiplicity (if available) | |
517 | var[0] = sum[0].M(); | |
518 | if (mc) var[1] = (sum[1].M() - sum[0].M()) / sum[1].M(); | |
519 | var[2] = sum[0].Perp(); | |
520 | ||
521 | // fill containers | |
522 | AliCFContainer *cf = 0x0; | |
523 | if (tr1->Charge() != tr2->Charge()) | |
524 | cf = fCFunlike; | |
525 | else if (tr1->Charge() > 0) | |
526 | cf = fCFlikePP; | |
527 | else | |
528 | cf = fCFlikeMM; | |
529 | ||
530 | // fill steps: 0 = quality onl1, 1 = quality + PID | |
531 | cf->Fill(var, 0); | |
532 | if (pid[i1] && pid[i2]) cf->Fill(var, 1); | |
0d73200d | 533 | } |
7356f978 | 534 | } |
5faf5a07 | 535 | } |
6279d59e | 536 | |
7356f978 | 537 | //__________________________________________________________________________________________________ |
538 | void AliRsnAnalysisPhi7TeV::ProcessMC(AliESDEvent *esd, AliMCEvent *mc) | |
5faf5a07 | 539 | { |
540 | // | |
7356f978 | 541 | // Searches true pairs and fill efficiency for kaons |
542 | // | |
543 | ||
544 | // get stack, if possible | |
545 | AliStack *stack = 0x0; | |
546 | if (mc) stack = mc->Stack(); else return; | |
547 | ||
548 | // count tracks | |
549 | Int_t ipart, itrack, npart = stack->GetNprimary(); | |
550 | AliESDtrack *matched[2] = {0x0, 0x0}; | |
551 | Double_t var1[10], var2[5]; | |
552 | Float_t b[2], bcov[3]; | |
553 | ||
554 | var1[7] = var2[3] = AliESDtrackCuts::GetReferenceMultiplicity(esd, kTRUE); | |
555 | var1[8] = var2[4] = npart; | |
556 | ||
557 | // loop on kaons | |
558 | for (ipart = 0; ipart < npart; ipart++) { | |
559 | TParticle *part = stack->Particle(ipart); | |
560 | if (!stack->IsPhysicalPrimary(ipart)) continue; | |
561 | if (TMath::Abs(part->GetPdgCode()) != 321) continue; | |
562 | ||
563 | // search best reconstructed track | |
564 | // which matches that particle ID | |
565 | // and checks it against cuts | |
566 | matched[0] = 0x0; | |
567 | itrack = BestMatchedTrack(esd, ipart); | |
568 | if (itrack >= 0) matched[0] = esd->GetTrack(itrack); | |
569 | ||
570 | // compute all variables fo single-tracks | |
571 | // [0] = pt | |
572 | // [1] = pTPC | |
573 | // [2] = DCAr | |
574 | // [3] = DCAz | |
575 | // [4] = signalITS | |
576 | // [5] = signalTPC | |
577 | // [6] = betaTOF | |
578 | // [7] = rec. multiplicity | |
579 | // [8] = MC multiplicity (if available) | |
580 | var1[0] = part->Pt(); | |
581 | var1[1] = -1.0; | |
582 | for (Int_t i = 0; i < 6; i++) var1[i] = -1.0; | |
583 | if (matched[0]) { | |
584 | matched[0]->GetImpactParameters(b, bcov); | |
585 | if (IsTPC(matched[0])) var1[1] = matched[0]->GetInnerParam()->P(); | |
586 | var1[2] = TMath::Abs((Double_t)b[0]); | |
587 | var1[3] = TMath::Abs((Double_t)b[1]); | |
588 | var1[4] = matched[0]->GetITSsignal(); | |
589 | var1[5] = matched[0]->GetTPCsignal(); | |
590 | if (MatchTOF(matched[0])) { | |
591 | var1[6] = matched[0]->GetIntegratedLength(); | |
592 | var1[6] /= matched[0]->GetTOFsignal(); | |
593 | var1[6] /= 0.299792458E8; | |
594 | } | |
595 | } | |
596 | ||
597 | // fill container | |
598 | fCFkaons->Fill(var1, 0); | |
599 | ||
600 | // replace MC momentum with rec | |
601 | if (matched[0]) { | |
602 | var1[0] = matched[0]->Pt(); | |
603 | if (OkQuality(matched[0])) { | |
604 | fCFkaons->Fill(var1, 1); | |
605 | if (OkPID(matched[0])) fCFkaons->Fill(var1, 2); | |
606 | } | |
607 | } | |
608 | } | |
609 | ||
610 | // loop on phi's | |
611 | Int_t id1, id2; | |
612 | TParticle *d1, *d2; | |
613 | TLorentzVector p1[2], p2[2], sum[2]; | |
614 | for (ipart = 0; ipart < npart; ipart++) { | |
615 | TParticle *part = stack->Particle(ipart); | |
616 | ||
617 | // check that particle is a phi | |
618 | // and that has decayed into kaons | |
619 | if (TMath::Abs(part->GetPdgCode()) != 333) continue; | |
620 | id1 = part->GetFirstDaughter(); | |
621 | id2 = part->GetLastDaughter(); | |
622 | d1 = d2 = 0x0; | |
623 | if (id1 >= 0 && id1 < npart) d1 = stack->Particle(id1); | |
624 | if (id2 >= 0 && id2 < npart) d2 = stack->Particle(id2); | |
625 | if (!d1 || !d2) continue; | |
626 | if (TMath::Abs(d1->GetPdgCode()) != 321) continue; | |
627 | if (TMath::Abs(d2->GetPdgCode()) != 321) continue; | |
628 | ||
629 | // get matched tracks | |
630 | matched[0] = matched[1] = 0x0; | |
631 | itrack = BestMatchedTrack(esd, id1); | |
632 | if (itrack >= 0) matched[0] = esd->GetTrack(itrack); | |
633 | itrack = BestMatchedTrack(esd, id2); | |
634 | if (itrack >= 0) matched[1] = esd->GetTrack(itrack); | |
635 | ||
636 | // compute values | |
637 | // [0] = inv. mass | |
638 | // [1] = inv. mass resolution | |
639 | // [2] = pt | |
640 | // [3] = rec. multiplicity | |
641 | // [4] = MC multiplicity (if available) | |
642 | p1[0].SetXYZM(d1->Px(), d1->Py(), d1->Pz(), fKaonMass); | |
643 | p2[0].SetXYZM(d2->Px(), d2->Py(), d2->Pz(), fKaonMass); | |
644 | p1[1].SetXYZM(0.0, 0.0, 0.0, 0.0); | |
645 | p2[1].SetXYZM(0.0, 0.0, 0.0, 0.0); | |
646 | if (matched[0]) p1[1].SetXYZM(matched[0]->Px(), matched[0]->Py(), matched[0]->Pz(), fKaonMass); | |
647 | if (matched[1]) p2[1].SetXYZM(matched[1]->Px(), matched[1]->Py(), matched[1]->Pz(), fKaonMass); | |
648 | for (Int_t i = 0; i < 2; i++) sum[i] = p1[i] + p2[i]; | |
649 | var2[0] = sum[0].M(); | |
650 | var2[1] = (sum[0].M() - sum[1].M()) / sum[0].M(); | |
651 | var2[2] = sum[0].Perp(); | |
652 | ||
653 | // fill container | |
654 | fCFtrues->Fill(var2, 0); | |
655 | ||
656 | // replace MC momenta with rec | |
657 | if (matched[0] && matched[1]) { | |
658 | var2[0] = sum[1].M(); | |
659 | var2[2] = sum[1].Perp(); | |
660 | if (OkQuality(matched[0]) && OkQuality(matched[1])) { | |
661 | fCFtrues->Fill(var2, 1); | |
662 | if (OkPID(matched[0]) && OkPID(matched[1])) fCFtrues->Fill(var2, 2); | |
663 | } | |
664 | } | |
665 | } | |
5faf5a07 | 666 | } |
6279d59e | 667 | |
5faf5a07 | 668 | //______________________________________________________________________________ |
7356f978 | 669 | Bool_t AliRsnAnalysisPhi7TeV::OkPIDITS(AliESDtrack *track, AliPID::EParticleType pid) |
5faf5a07 | 670 | { |
671 | // | |
672 | // Check ITS particle identification with 3sigma cut | |
673 | // | |
4c06ce33 | 674 | |
7356f978 | 675 | // reject not ITS standalone tracks |
676 | if (!IsITS(track)) return kFALSE; | |
677 | ||
678 | // count PID layers and reject if they are too few | |
679 | Int_t k, nITSpidLayers = 0; | |
680 | UChar_t itsCluMap = track->GetITSClusterMap(); | |
681 | for (k = 2; k < 6; k++) if (itsCluMap & (1 << k)) ++nITSpidLayers; | |
682 | if (nITSpidLayers < 3) { | |
683 | AliDebug(AliLog::kDebug + 2, "Rejecting track with too few ITS pid layers"); | |
684 | return kFALSE; | |
685 | } | |
686 | ||
687 | // check the track type (ITS+TPC or ITS standalone) | |
688 | // and reject it if it is of none of the allowed types | |
689 | Bool_t isSA = kFALSE; | |
690 | if (IsTPC(track)) isSA = kFALSE; | |
691 | else if (IsITS(track)) isSA = kTRUE; | |
692 | else { | |
693 | AliWarning("Track is neither ITS+TPC nor ITS standalone"); | |
694 | return kFALSE; | |
695 | } | |
696 | ||
697 | // create the PID response object and compute nsigma | |
698 | AliITSPIDResponse &itsrsp = fESDpid->GetITSResponse(); | |
699 | Double_t mom = track->P(); | |
700 | Double_t nSigma = itsrsp.GetNumberOfSigmas(mom, track->GetITSsignal(), pid, nITSpidLayers, isSA); | |
701 | ||
702 | // evaluate the cut | |
703 | Bool_t ok = (TMath::Abs(nSigma) <= fITSNSigma); | |
704 | ||
705 | // debug message | |
706 | AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", nSigma, fITSNSigma, (ok ? "passed" : "failed"))); | |
707 | ||
708 | // outcome | |
709 | return ok; | |
5faf5a07 | 710 | } |
4c06ce33 | 711 | |
5faf5a07 | 712 | //______________________________________________________________________________ |
7356f978 | 713 | Bool_t AliRsnAnalysisPhi7TeV::OkPIDTPC(AliESDtrack *track, AliPID::EParticleType pid) |
5faf5a07 | 714 | { |
715 | // | |
716 | // Check TPC particle identification with {3|5}sigmacut, | |
717 | // depending on the track total momentum. | |
718 | // | |
4c06ce33 | 719 | |
7356f978 | 720 | // reject not TPC tracks |
721 | if (!IsTPC(track)) return kFALSE; | |
722 | if (!track->GetInnerParam()) return kFALSE; | |
6279d59e | 723 | |
7356f978 | 724 | // setup TPC PID response |
725 | AliTPCPIDResponse &tpcrsp = fESDpid->GetTPCResponse(); | |
a3dd4b4d | 726 | |
7356f978 | 727 | // get momentum and number of sigmas and choose the reference band |
728 | Double_t mom = track->GetInnerParam()->P(); | |
729 | Double_t nSigma = tpcrsp.GetNumberOfSigmas(mom, track->GetTPCsignal(), track->GetTPCsignalN(), pid); | |
730 | Double_t maxNSigma = fTPCNSigma[0]; | |
731 | if (mom > fTPCMomentumThreshold) maxNSigma = fTPCNSigma[1]; | |
a3dd4b4d | 732 | |
7356f978 | 733 | // evaluate the cut |
734 | Bool_t ok = (TMath::Abs(nSigma) <= maxNSigma); | |
a3dd4b4d | 735 | |
7356f978 | 736 | // debug message |
737 | AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", nSigma, maxNSigma, (ok ? "passed" : "failed"))); | |
a3dd4b4d | 738 | |
7356f978 | 739 | // outcome |
740 | return ok; | |
6279d59e | 741 | } |
742 | ||
5faf5a07 | 743 | //______________________________________________________________________________ |
7356f978 | 744 | Bool_t AliRsnAnalysisPhi7TeV::OkPIDTOF(AliESDtrack *track, AliPID::EParticleType pid) |
0d73200d | 745 | { |
746 | // | |
5faf5a07 | 747 | // Check TOF particle identification if matched there. |
0d73200d | 748 | // |
749 | ||
7356f978 | 750 | // reject not TOF-matched tracks |
751 | if (!MatchTOF(track)) return kFALSE; | |
a3dd4b4d | 752 | |
7356f978 | 753 | // setup TOF PID response |
754 | AliTOFPIDResponse &tofrsp = fESDpid->GetTOFResponse(); | |
a3dd4b4d | 755 | |
7356f978 | 756 | // get info for computation |
757 | Double_t momentum = track->P(); | |
758 | Double_t time = track->GetTOFsignal(); | |
759 | Double_t timeint[AliPID::kSPECIES]; | |
760 | tofrsp.GetStartTime(momentum); | |
761 | track->GetIntegratedTimes(timeint); | |
5faf5a07 | 762 | |
7356f978 | 763 | // check the cut |
764 | Double_t timeDiff = time - timeint[(Int_t)pid]; | |
765 | Double_t sigmaRef = tofrsp.GetExpectedSigma(momentum, timeint[(Int_t)pid], AliPID::ParticleMass(pid)); | |
766 | Double_t nSigma = timeDiff / sigmaRef; | |
a3dd4b4d | 767 | |
7356f978 | 768 | // evaluate the cut |
769 | Bool_t ok = (TMath::Abs(nSigma) <= fTOFNSigma); | |
a3dd4b4d | 770 | |
7356f978 | 771 | // debug message |
772 | AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- nsigma = %f -- cut %s", nSigma, fTOFNSigma, (ok ? "passed" : "failed"))); | |
a3dd4b4d | 773 | |
7356f978 | 774 | // outcome |
775 | return ok; | |
5faf5a07 | 776 | } |
777 | ||
778 | //______________________________________________________________________________ | |
7356f978 | 779 | Int_t AliRsnAnalysisPhi7TeV::MatchedTrack(AliESDEvent *esd, Int_t label, Int_t &npassed, Int_t start) |
6279d59e | 780 | { |
781 | // | |
7356f978 | 782 | // Starting from 'start' searches a track in ESD which has label 'label', |
783 | // and if it is found, records how many cuts it passes: | |
784 | // 0 = none | |
785 | // 1 = only quality | |
786 | // 2 = quality and PID | |
787 | // Returns the ID of that track (-1 if not found) | |
788 | // | |
789 | ||
790 | Int_t it, ntracks = esd->GetNumberOfTracks(); | |
791 | AliESDtrack *track = 0x0; | |
792 | ||
793 | for (it = start; it < ntracks; it++) { | |
794 | track = esd->GetTrack(it); | |
795 | if (TMath::Abs(track->GetLabel() != label)) continue; | |
796 | ||
797 | npassed = 0; | |
798 | if (OkQuality(track)) npassed++; | |
799 | if (OkPID (track)) npassed++; | |
800 | return it; | |
801 | } | |
802 | ||
803 | return -1; | |
804 | } | |
a3dd4b4d | 805 | |
7356f978 | 806 | //______________________________________________________________________________ |
807 | Int_t AliRsnAnalysisPhi7TeV::BestMatchedTrack(AliESDEvent *esd, Int_t label) | |
808 | { | |
809 | // | |
810 | // Searches the best track which matches the specified label | |
811 | // where 'best' means the one which passes most cuts | |
812 | // | |
813 | ||
814 | // first attempt: | |
815 | // if it fails, there are no matched tracks | |
816 | Int_t ncuts = 0; | |
fbcdae2f | 817 | Int_t itrack = MatchedTrack(esd, label, ncuts); |
7356f978 | 818 | /*if (itrack < 0) return -1; |
819 | ||
820 | // if it succeeds, use it as a start for a loop | |
821 | Int_t bestCuts = ncuts; | |
822 | Int_t start = itrack + 1; | |
823 | for (;;) { | |
824 | it = MatchedTrack(esd, label, ncuts, start); | |
825 | if (it < 0) break; | |
826 | ||
827 | if (ncuts > bestCuts) { | |
828 | bestCuts = ncuts; | |
829 | itrack = it; | |
830 | start = it + 1; | |
831 | } | |
832 | }*/ | |
833 | ||
834 | return itrack; | |
6279d59e | 835 | } |