]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes for wrong use of const causing PW.CAST_TO_QUALIFIED_TYPE defect in Coverity
[u/mrichter/AliRoot.git] / PWGUD / DIFFRACTIVE / xsAndTwoProng / AliCDMesonUtils.cxx
CommitLineData
369562e9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//
16// AliCDMesonUtils
17// for
18// AliAnalysisTaskCDMeson
19//
20// Author:
21// Xianguo Lu <lu@physi.uni-heidelberg.de>
22// continued by
23// Felix Reidt <Felix.Reidt@cern.ch>
24
25#include <TH1.h>
26#include <TH2.h>
27#include <TGeoMatrix.h>
28#include <THnSparse.h>
29#include <TLorentzVector.h>
30#include <TRandom3.h>
31#include <TParticle.h>
32#include <TObjArray.h>
33
34#include "AliCDBManager.h"
35#include "AliCDBEntry.h"
36#include "AliCDBStorage.h"
37#include "AliESDEvent.h"
38#include "AliPIDResponse.h"
39#include "AliVTrack.h"
40#include "AliVParticle.h"
41#include "AliESDtrack.h"
42#include "AliESDtrackCuts.h"
43#include "AliESDFMD.h"
44#include "AliESDVZERO.h"
45#include "AliGeomManager.h"
46#include "AliITSAlignMille2Module.h"
47#include "AliITSsegmentationSPD.h"
48#include "AliMultiplicity.h"
49#include "AliPIDResponse.h"
50#include "AliSPDUtils.h"
51#include "AliTriggerAnalysis.h"
52#include "AliVZEROTriggerData.h"
53#include "AliVZEROCalibData.h"
54
55#include "AliAODTracklets.h"
56#include "AliAODEvent.h"
57
58#include "AliCDMesonBase.h"
59#include "AliCDMesonTracks.h"
60
61#include "AliCDMesonUtils.h"
62
63
64//==============================================================================
65//------------------------------------------------------------------------------
34e8fe1d 66void AliCDMesonUtils::GetMassPtCtsOA(Int_t pid,
369562e9 67 const TVector3* momenta[], Double_t & mass,
68 Double_t &pt, Double_t &cts, Double_t &oa)
69{
70 //
71 // Get Mass Pt Cts OA
72 //
73
74 Double_t tmpp1[3], tmpp2[3];
75 momenta[0]->GetXYZ(tmpp1);
76 momenta[1]->GetXYZ(tmpp2);
77
78 Double_t masstrk = -999;
79 if(pid==AliCDMesonBase::kBinPionE || pid==AliCDMesonBase::kBinPion ||
80 pid==AliCDMesonBase::kBinSinglePion)
81 masstrk = AliPID::ParticleMass(AliPID::kPion);
82 else if(pid==AliCDMesonBase::kBinKaonE || pid==AliCDMesonBase::kBinKaon ||
83 pid==AliCDMesonBase::kBinSingleKaon)
84 masstrk = AliPID::ParticleMass(AliPID::kKaon);
85 else if(pid==AliCDMesonBase::kBinProtonE || pid==AliCDMesonBase::kBinProton ||
86 pid==AliCDMesonBase::kBinSingleProton)
87 masstrk = AliPID::ParticleMass(AliPID::kProton);
88 else if(pid==AliCDMesonBase::kBinElectronE ||
89 pid==AliCDMesonBase::kBinElectron ||
90 pid==AliCDMesonBase::kBinSingleElectron)
91 masstrk = AliPID::ParticleMass(AliPID::kElectron);
92 else
93 masstrk = AliPID::ParticleMass(AliPID::kPion);
94
95 const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masstrk, masstrk,
96 cts);
97 mass = sumv.M();
98 pt = sumv.Pt();
99
100 oa = GetOA(tmpp1, tmpp2);
101}
102
103
104//------------------------------------------------------------------------------
34e8fe1d 105void AliCDMesonUtils::GetPWAinfo(Int_t pid, const AliVTrack *trks[],
369562e9 106 Float_t& theta, Float_t& phi, Float_t& mass,
107 Float_t momentum[])
108{
109 // transforms the coordiante system to the helicity frame and determines
110 // invariant mass and 3-vector of the two track system
111 //
112
113 Double_t tmpp1[3], tmpp2[3]; // 3-vectors of the daughter tracks
114 trks[0]->GetPxPyPz(tmpp1);
115 trks[1]->GetPxPyPz(tmpp2);
116
117 // determine mass of the daugther tracks
118 Double_t masstrk = -999;
119 if(pid==AliCDMesonBase::kBinPion || pid==AliCDMesonBase::kBinSinglePion)
120 masstrk = AliPID::ParticleMass(AliPID::kPion);
121 else if(pid==AliCDMesonBase::kBinKaon || pid==AliCDMesonBase::kBinSingleKaon)
122 masstrk = AliPID::ParticleMass(AliPID::kKaon);
123 else if(pid==AliCDMesonBase::kBinProton ||
124 pid==AliCDMesonBase::kBinSingleProton)
125 masstrk = AliPID::ParticleMass(AliPID::kProton);
126 else if(pid==AliCDMesonBase::kBinElectron ||
127 pid==AliCDMesonBase::kBinSingleElectron)
128 masstrk = AliPID::ParticleMass(AliPID::kElectron);
129 else
130 masstrk = AliPID::ParticleMass(AliPID::kPion);
131
132 TLorentzVector va, vb; // 4-vectors of the two tracks
133 va.SetXYZM(tmpp1[0], tmpp1[1], tmpp1[2], masstrk);
134 vb.SetXYZM(tmpp2[0], tmpp2[1], tmpp2[2], masstrk);
135
136 TLorentzVector sumv = va+vb; // 4-vector of the pair
137 TVector3 sumXZ(sumv.X(), 0, sumv.Z()); // its projection to the xz-plane
138
139 // mass and momentum
140 mass = sumv.M();
141 momentum[0] = sumv.X();
142 momentum[1] = sumv.Y();
143 momentum[2] = sumv.Z();
144
145 // coordinate axis vectors
146 const TVector3 x(1.,0,0);
147 const TVector3 y(0,1.,0);
148 const TVector3 z(0,0,1.);
149
150 const Double_t alpha = -z.Angle(sumXZ);
151
152 sumXZ.Rotate(alpha, y);
153 sumv.Rotate(alpha, y);
154 va.Rotate(alpha, y);
155 vb.Rotate(alpha, y);
156
157 const Double_t beta = z.Angle(sumv.Vect());
158 sumv.Rotate(beta, x);
159 va.Rotate(beta, x);
160 vb.Rotate(beta, x);
161
162
163 const TVector3 bv = -sumv.BoostVector();
164 va.Boost(bv);
165 vb.Boost(bv);
166
167 // angles
168 theta = va.Theta();
169 phi = va.Phi();
170}
171
172//------------------------------------------------------------------------------
173Int_t AliCDMesonUtils::GetEventType(const AliESDEvent *ESDEvent)
174{
175 // checks of which type a event is:
176 // beam-beam interaction (I), beam-gas (A/C), empty (E)
177 //
178
179 TString firedTriggerClasses = ESDEvent->GetFiredTriggerClasses();
180
181 if (firedTriggerClasses.Contains("CINT1A-ABCE-NOPF-ALL")) { // A
182 return AliCDMesonBase::kBinEventA;
183 }
184 if (firedTriggerClasses.Contains("CINT1C-ABCE-NOPF-ALL")) { // C
185 return AliCDMesonBase::kBinEventC;
186 }
187 if (firedTriggerClasses.Contains("CINT1B-ABCE-NOPF-ALL")) { // I
188 return AliCDMesonBase::kBinEventI;
189 }
190 if (firedTriggerClasses.Contains("CINT1-E-NOPF-ALL")) { // E
191 return AliCDMesonBase::kBinEventE;
192 }
193 if (firedTriggerClasses.Contains("CDG5-E")) { // E
194 return AliCDMesonBase::kBinEventE;
195 }
196 if (firedTriggerClasses.Contains("CDG5-I")) { // I
197 return AliCDMesonBase::kBinEventI;
198 }
199 if (firedTriggerClasses.Contains("CDG5-AC")) { // AC
200 return AliCDMesonBase::kBinEventAC;
201 }
202 return AliCDMesonBase::kBinEventUnknown;
203}
204
205
206//------------------------------------------------------------------------------
207void AliCDMesonUtils::SwapTrack(const AliVTrack* trks[])
208{
209 //
210 // swap two esdtracks, needed since only the properties of one a stored and
211 // AliRoot sorts them
212 //
213
214 TRandom3 tmprd(0);
215 if(tmprd.Rndm()>0.5)
216 return;
217
218 const AliVTrack* tmpt = trks[0];
219 trks[0]=trks[1];
220 trks[1]=tmpt;
221}
222
223
224//------------------------------------------------------------------------------
225Int_t AliCDMesonUtils::GetCombCh(const AliVTrack * trks[])
226{
227 //
228 // get combination of charges
229 //
230
231 const Double_t ch1 = trks[0]->Charge();
232 const Double_t ch2 = trks[1]->Charge();
233
234 if(ch1*ch2<0){
235 return AliCDMesonBase::kBinPM;
236 }
237 else
238 return AliCDMesonBase::kBinPPMM;
239}
240
241
242//------------------------------------------------------------------------------
243Int_t AliCDMesonUtils::GetCombPID(AliPIDResponse* pid, const AliVTrack* trks[],
34e8fe1d 244 Int_t mode, TH2* comb2trkPID /*= 0x0 */)
369562e9 245{
246 //
247 // Get PID for a two track system (data)
248 //
249
250 if (!pid){
251 return AliCDMesonBase::kBinPIDUnknown;
252 }
253
254 // get PID for the single tracks
255 Int_t kpid[2];
256 for(Int_t ii=0; ii<2; ii++){
257 kpid[ii] = GetPID(pid, trks[ii], mode);
258 }
259
260 // PID QA histogram
261 if (comb2trkPID) {
262 comb2trkPID->Fill(kpid[0], kpid[1]);
263 }
264
265 return CombinePID(kpid);
266}
267
268
269//------------------------------------------------------------------------------
270Int_t AliCDMesonUtils::GetCombPID(const TParticle* particles[],
271 TH2 *comb2trkPID /* = 0x0 */)
272{
273 //
274 // Get PID for a two track system (MC)
275 //
276
277 // get PID for the single tracks
278 Int_t kpid[2];
279 for(Int_t ii=0; ii<2; ii++){
280 kpid[ii] = GetPID(particles[ii]->GetPdgCode());
281 }
282
283 // PID QA histogram
284 if (comb2trkPID) {
285 comb2trkPID->Fill(kpid[0], kpid[1]);
286 }
287
288 return CombinePID(kpid);
289}
290
291
292//------------------------------------------------------------------------------
293void AliCDMesonUtils::GetSPDTrackletMult(const AliESDEvent *ESDEvent,
294 Int_t& sum, Int_t& forwardA,
295 Int_t& forwardC, Int_t& central)
296{
297 // obtain the multiplicity for eta < -.9 (forwardC), eta > 0.9 (fowardA), in
298 // the central barrel (central) and its sum from SPD tracklets with the
299 // AliMultiplicity class
300 // code simular to PWGPP/ITS/AliAnalysisTaskSPD.cxx
301
302 sum = forwardA = forwardC = central = 0; // initialize values
303
304 const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
305 for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
306 iTracklet++) {
307 float_t eta = mult->GetEta(iTracklet);
308 if (eta < -0.9) {
309 forwardC++;
310 }
311 else if (eta < 0.9) {
312 central++;
313 }
314 else {
315 forwardA++;
316 }
317 sum++;
318 }
319}
320
321
322//------------------------------------------------------------------------------
323Bool_t AliCDMesonUtils::CutEvent(const AliESDEvent *ESDEvent, TH1 *hspd,
324 TH1 *hpriv, TH1 *hpriVtxPos, TH1 *hpriVtxDist,
325 TH2 *hfo, TH1* hfochans, Int_t &kfo,
326 Int_t &nip, Int_t &nop, TH1 *hpriVtxX,
327 TH1 *hpriVtxY, TH1 *hpriVtxZ)
328{
329 //
330 // CutEvent
331 //
332
333 AliTriggerAnalysis triggerAnalysis;
334 /*
649e06c4 335 //printf("AliCDMesonUtils active triggers: %s\n",
369562e9 336 // ESDEvent->GetHeader()->GetActiveTriggerInputs().Data());
649e06c4 337 //printf("AliCDMesonUtils fired triggers: %s\n",
369562e9 338 // ESDEvent->GetHeader()->GetFiredTriggerInputs().Data());
339 //*/
340 //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
341 //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
342 // returns the number of fired chips in the SPD
343 // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
344 // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
345
346
347 // SPD FastOR cut
348 const Int_t fastORhw = triggerAnalysis.SPDFiredChips(ESDEvent, 1);
349 if (hspd) hspd->Fill(fastORhw);
350 /*
351 if(fastORhw<1)
352 return kFALSE;
353 */
354
355 if (hfochans) {
356 const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
357
358 for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
359 if (mult->TestFastOrFiredChips(iChipKey)) {
360 hfochans->Fill((Double_t)iChipKey);
361 }
362 }
363 }
364
365 if(kfo){ // spd trigger study
366 Int_t nfoctr[10];
367 GetNFO(ESDEvent, "[1.0]", nfoctr, 0x0, 0x0);
368 nip = nfoctr[kInnerPixel];
369 nop = nfoctr[kOuterPixel];
370
371 if(AliCDMesonBase::GetGapBin("V0", GetV0(ESDEvent)) ==
372 AliCDMesonBase::kBinDG) {
373 hfo->Fill(nip, nop);
374 }
375
376 if(nop<2)
377 kfo = 1;
378 else
379 kfo = nop;
380
381 if(kfo>=10)
382 kfo=9;
383 }
384
385 // collision vertex cut
386 // A cut in XY is implicitly done during the reconstruction by constraining
387 // the vertex to the beam diamond.
388
389 // Primary vertex
390 Bool_t kpr0 = kTRUE;
391 const AliESDVertex *vertex = ESDEvent->GetPrimaryVertexTracks();
392 if(vertex->GetNContributors()<1) {
393 // SPD vertex
394 vertex = ESDEvent->GetPrimaryVertexSPD();
395 if(vertex->GetNContributors()<1) {
396 // NO VERTEX, SKIP EVENT
397 kpr0 = kFALSE;
398 }
399 }
400 const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ()) < 10.);
401 // 10 is the common value, unit: cm
402 if (hpriv) hpriv->Fill(kpriv);
403 if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
404 if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
405 if(!kpriv)
406 return kFALSE;
407
408 if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
409 if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
410 if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
411 return kTRUE;
412}
413
414//------------------------------------------------------------------------------
415Bool_t AliCDMesonUtils::CutEvent(const AliAODEvent *AODEvent, TH1 *hpriv,
416 TH1 *hpriVtxX, TH1 *hpriVtxY, TH1 *hpriVtxZ,
417 TH1 *hpriVtxPos, TH1 *hpriVtxDist)
418{
419 //
420 // Cut Event for AOD Events, to be combined with the ESD Track Cut
421 //
422
423 // TODO: no idea about fast or yet, to be thought of
424
425 // Primary vertex
426 Bool_t kpr0 = kTRUE;
427 const AliAODVertex *vertex = AODEvent->GetPrimaryVertex();
428 if(vertex->GetNContributors()<1) {
429 // SPD vertex
430 vertex = AODEvent->GetPrimaryVertexSPD();
431 if(vertex->GetNContributors()<1) {
432 // NO GOOD VERTEX, SKIP EVENT
433 kpr0 = kFALSE;
434 }
435 }
436 const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ())<10.);
437 // 10 is the common value, unit: cm
438 hpriv->Fill(kpriv);
439 if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
440 if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
441 if(!kpriv)
442 return kFALSE;
443
444 if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
445 if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
446 if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
447 return kTRUE;
448}
449
450
451//------------------------------------------------------------------------------
452void AliCDMesonUtils::DoVZEROStudy(const AliESDEvent *ESDEvent,
34e8fe1d 453 TObjArray* hists, Int_t run)
369562e9 454{
455 //
456 //
457 // IMPORTANT: the order of the histograms here and in
458 // AliCDMesonBase::GetHistVZEROStudies(..) has to match
459
460 const AliESDVZERO* esdV0 = ESDEvent->GetVZEROData();
461 if (!esdV0) {
462 Printf("ERROR: esd V0 not available");
463 return;
464 }
465
466 // determine trigger decision
467 AliTriggerAnalysis triggerAnalysis;
468 //const Bool_t khw = kFALSE;
469
470 //Float_t v0a = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
471 // AliTriggerAnalysis::kV0BB) ? 1. : 0.;
472 //Float_t v0c = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
473 // AliTriggerAnalysis::kV0BB) ? 1. : 0.;
474
475 // obtain OCDB objects
476 AliCDBManager *man = AliCDBManager::Instance();
477 TString cdbpath;
478 if (man->IsDefaultStorageSet()) {
479 const AliCDBStorage *dsto = man->GetDefaultStorage();
480 cdbpath = TString(dsto->GetBaseFolder());
481 }
482 else { //should not be used!
483 man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
484 cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
485 }
486 man->SetSpecificStorage("VZERO/Trigger/Data",cdbpath);
487 man->SetSpecificStorage("VZERO/Calib/Data",cdbpath);
488 man->SetRun(run);
489
490 AliCDBEntry *ent1 = man->Get("VZERO/Trigger/Data");
649e06c4 491 if (!ent1) {
65711644 492 printf("AliCDMesonUtils failed loading VZERO trigger entry from OCDB\n");
493 return;
649e06c4 494 }
369562e9 495 AliVZEROTriggerData *trigData = (AliVZEROTriggerData*)ent1->GetObject();
496 if (!trigData) {
649e06c4 497 printf("AliCDMesonUtils failed loading VZERO trigger data from OCDB\n");
369562e9 498 return;
499 }
500
501 AliCDBEntry *ent2 = man->Get("VZERO/Calib/Data");
65711644 502 if (!ent2) {
503 printf("AliCDMesonUtils failed loading VZERO calib entry from OCDB\n");
504 return;
505 }
369562e9 506 AliVZEROCalibData *calData = (AliVZEROCalibData*)ent2->GetObject();
507 if (!calData) {
649e06c4 508 printf("AliCDMesonUtils failed loading VZERO calibration data from OCDB\n");
369562e9 509 return;
510 }
511
512 // fill histograms
513 Int_t pmtHist = 0;
514 Int_t iPMT = 0;
515 for (iPMT = 0; iPMT < 64; ++iPMT) {
516 Int_t board = AliVZEROCalibData::GetBoardNumber(iPMT);
517 Int_t channel = AliVZEROCalibData::GetFEEChannelNumber(iPMT);
518 ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
519 trigData->GetPedestalCut(0, board, channel));
520 // calData->GetDiscriThr(iPMT));
521 }
522 pmtHist = iPMT;
523 for (iPMT = 0; iPMT < 64; ++iPMT) {
524 ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
525 esdV0->GetMultiplicity(iPMT));
526 }
527 /*
528 // not used in 2010 for pp
529 ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeA(), v0a);
530 ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeC(), v0c);
531 */
532}
533
534
535//------------------------------------------------------------------------------
536Int_t AliCDMesonUtils::GetGapConfig(const AliESDEvent *ESDEvent,
537 TH2 *hitMapSPDinner, TH2 *hitMapSPDouter,
538 TH2 *hitMapSPDtrklt, TH2 *hitMapFMDa,
539 TH2 *hitMapFMDc, TH1 **fmdSums,
540 TH2 *TPCGapDCAaSide, TH2 *TPCGapDCAcSide)
541{
542 //
543 // GetGapConfigAndTracks
544 //
545 // retrieves the gap configuration of a track and returns it as
546 // an bit vector
547 // kBaseLine ensures, that this event is valid
548 // + is equivalent to | in this case
549 Int_t gapConfig = AliCDMesonBase::kBitBaseLine + GetV0(ESDEvent)
550 + GetFMD(ESDEvent, hitMapFMDa, hitMapFMDc, fmdSums)
551 + GetSPD(ESDEvent, hitMapSPDinner, hitMapSPDouter, hitMapSPDtrklt);
552 if (gapConfig == AliCDMesonBase::kBitBaseLine) {
553 gapConfig += GetTPC(ESDEvent, TPCGapDCAaSide, TPCGapDCAcSide);
554 }
555 else {
556 gapConfig += GetTPC(ESDEvent, 0x0, 0x0);
557 }
558 if (GetFastORmultiplicity(ESDEvent) > 0) {
559 gapConfig += AliCDMesonBase::kBitCentAct;
560 }
561 return gapConfig; // + GetZDC(ESDEvent);
562}
563
564
565//------------------------------------------------------------------------------
566void AliCDMesonUtils::FillEtaPhiMap(const AliVEvent *event,
567 const AliCDMesonTracks* tracks, TH2 *map,
568 TH2 *map_c)
569{
570 //
571 // Fills the eta phi information about all tracks and about the tracks surving
572 // the track cuts (CutTrack) into two separate histograms
573 //
574
575 // all tracks
576 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); ++itrk) {
577 if (AliVParticle *trk = event->GetTrack(itrk)) {
578 if (map) {
579 map->Fill(trk->Eta(), trk->Phi());
580 }
581 }
582 }
583
584 // tracks that survived the cuts
585 for (Int_t itrk = 0; itrk < tracks->GetCombinedTracks(); ++itrk) {
586 if (map_c) {
587 map_c->Fill(tracks->GetTrack(itrk)->Eta(), tracks->GetTrack(itrk)->Phi());
588 }
589 }
590}
591
592
593//------------------------------------------------------------------------------
594void AliCDMesonUtils::GetMultFMD(const AliESDEvent *ESDEvent, Int_t& fmdA,
595 Int_t& fmdC, Float_t *fmdSums)
596{
597 //
598 // Multiplicity seen by FMD
599 //
600 // WARNING: this function is only working with a modified AliRoot so far
601
602#ifdef STD_ALIROOT
603 fmdA = FMDHitCombinations(ESDEvent, 0);
604 fmdC = FMDHitCombinations(ESDEvent, 1);
605#else
606 AliTriggerAnalysis triggerAnalysis;
607 triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
608 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide, &fmdA);
609 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide, &fmdC);
610#endif
611
612 if (fmdSums) {
613 const AliESDFMD* fmdData =
614 (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
615 for (UShort_t det=1; det<=3; det++) {
616 Int_t nRings = (det == 1 ? 1 : 2);
617 for (UShort_t ir = 0; ir < nRings; ir++) {
618 Char_t ring = (ir == 0 ? 'I' : 'O');
619 UShort_t nsec = (ir == 0 ? 20 : 40);
620 UShort_t nstr = (ir == 0 ? 512 : 256);
621 for (UShort_t sec =0; sec < nsec; sec++) {
622 for (UShort_t strip = 0; strip < nstr; strip++) {
623 Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
624 if (mult == AliESDFMD::kInvalidMult) continue;
625
626 if (det == 1 && ring == 'I') fmdSums[0] += mult;
627 else if (det == 2 && ring == 'I') fmdSums[1] += mult;
628 else if (det == 2 && ring == 'O') fmdSums[2] += mult;
629 else if (det == 3 && ring == 'I') fmdSums[3] += mult;
630 else if (det == 3 && ring == 'O') fmdSums[4] += mult;
631 }
632 }
633 }
634 }
635 }
636}
637
638
639//------------------------------------------------------------------------------
640void AliCDMesonUtils::GetMultSPD(const AliESDEvent * ESDEvent, Int_t& spdIA,
641 Int_t& spdIC, Int_t& spdOA, Int_t& spdOC)
642{
643 //
644 // Retrieves the multiplicity seen by the SPD FastOR
645 //
646
647 Int_t nfoctr[10];
648 GetNFO(ESDEvent, "]0.9[", nfoctr, 0x0, 0x0);
649
650 spdIA = nfoctr[kIPA];
651 spdIC = nfoctr[kIPC];
652 spdOA = nfoctr[kOPA];
653 spdOC = nfoctr[kOPC];
654}
655
656
657//==============================================================================
658//------------------------------------------------------------------------------
659Int_t AliCDMesonUtils::GetV0(const AliESDEvent * ESDEvent)
660{
661 //
662 //GetV0
663 //
664
665 AliTriggerAnalysis triggerAnalysis;
666 const Bool_t khw = kFALSE;
667 const Bool_t v0A =
668 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
669 AliTriggerAnalysis::kV0BB);
670 const Bool_t v0C =
671 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
672 AliTriggerAnalysis::kV0BB);
673
674 return v0A * AliCDMesonBase::kBitV0A + v0C * AliCDMesonBase::kBitV0C;
675}
676
677
678//------------------------------------------------------------------------------
679Int_t AliCDMesonUtils::GetFMD(const AliESDEvent *ESDEvent, TH2 *hitMapFMDa,
680 TH2 *hitMapFMDc, TH1 **fmdSums)
681{
682 //
683 // GetFMD
684 //
685
686 AliTriggerAnalysis triggerAnalysis;
687 triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
688 const Bool_t fmdA =
689 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide);
690 const Bool_t fmdC =
691 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide);
692
693 //printf("FR - GetFMD\n");
694
695 // prepartions for a charge summation algorithm
696 Bool_t hitMaps = (Bool_t)(hitMapFMDa && hitMapFMDc);
697 Bool_t calcSum = fmdSums ?
698 (fmdSums[0] && fmdSums[1] && fmdSums[2] && fmdSums[3] && fmdSums[4]) :
699 kFALSE; // are the histograms defined?
700 Float_t sum[] = { 0., 0., 0., 0., 0. };
701 // summed multiplicity in the FMD detectors
702
703 // hit map generation
704 if (hitMaps || calcSum) {
705 const AliESDFMD* fmdData =
706 (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
707
708 for (UShort_t det=1; det<=3;det++) {
709 Int_t nRings = (det == 1 ? 1 : 2);
710 for (UShort_t ir = 0; ir < nRings; ir++) {
711 Char_t ring = (ir == 0 ? 'I' : 'O');
712 UShort_t nsec = (ir == 0 ? 20 : 40);
713 UShort_t nstr = (ir == 0 ? 512 : 256);
714 for (UShort_t sec =0; sec < nsec; sec++) {
715 for (UShort_t strip = 0; strip < nstr; strip++) {
716 Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
717 if (mult == AliESDFMD::kInvalidMult) continue;
718
719 if (calcSum) {
720 if (det == 1 && ring == 'I') sum[0] += mult;
721 else if (det == 2 && ring == 'I') sum[1] += mult;
722 else if (det == 2 && ring == 'O') sum[2] += mult;
723 else if (det == 3 && ring == 'I') sum[3] += mult;
724 else if (det == 3 && ring == 'O') sum[4] += mult;
725 }
726
727 if (hitMaps) { // care about hit map specific information
728 const Float_t eta = fmdData->Eta(det,ring,sec,strip);
729 const Float_t phi =
730 fmdData->Phi(det,ring,sec,strip) / 180. * TMath::Pi();
731 //printf("FR - GetFMD: %f %f %f\n", eta, phi, mult);
732 if (eta != AliESDFMD::kInvalidEta) {
733 if ((-3.5 < eta) && (eta < -1.5)) {
734 hitMapFMDc->Fill(eta, phi, mult); // use mult as weight
735 }
736 else if ((1.5 < eta) && (eta < 5.5)) {
737 hitMapFMDa->Fill(eta, phi, mult); // use mult as weight
738 }
739 }
740 }
741 }
742 }
743 }
744 }
745 }
746
747 if (calcSum) {
748 //printf("DEBUG -- SUM(%f,%f,%f,%f,%f)\n", sum[0], sum[1], sum[2], sum[3],
749 // sum[4]);
750 for (UInt_t i = 0; i < 5; i++) { //
751 fmdSums[i]->Fill(sum[i]);
752 }
753 }
754
755 return fmdA * AliCDMesonBase::kBitFMDA + fmdC * AliCDMesonBase::kBitFMDC;
756}
757
758
759//------------------------------------------------------------------------------
760Int_t AliCDMesonUtils::GetSPD(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
761 TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt)
762{
763 //
764 // GetSPD
765 //
766
767 Int_t nfoctr[10];
768 GetNFO(ESDEvent, "]0.9[", nfoctr, hitMapSPDinner, hitMapSPDouter);
769 // get multiplicity from fastOR and fill corresponding hit maps
770
771 if (hitMapSPDtrklt) FillSPDtrkltMap(ESDEvent, hitMapSPDtrklt);
772 // fill tracklet hit map
773
774 const Int_t ipA = nfoctr[kIPA]; // inner layer A side
775 const Int_t ipC = nfoctr[kIPC]; // inner layer C side
776 const Int_t opA = nfoctr[kOPA]; // outer layer A side
777 const Int_t opC = nfoctr[kOPC]; // outer layer C side
778
779 const Bool_t spdA = ipA + opA; // A side hit?
780 const Bool_t spdC = ipC + opC; // C side hit?
781
782 return spdA * AliCDMesonBase::kBitSPDA + spdC * AliCDMesonBase::kBitSPDC;
783}
784
785
786//------------------------------------------------------------------------------
787Int_t AliCDMesonUtils::GetTPC(const AliESDEvent * ESDEvent, TH2 *TPCGapDCAaSide,
788 TH2 *TPCGapDCAcSide)
789{
790 //
791 //GetTPC
792 //
793
794 const Double_t etacut = 0.9;
795 Int_t nA = 0;
796 Int_t nC = 0;
797
798 AliESDtrackCuts cuts;
799 cuts.SetMaxDCAToVertexXY(0.1);
800 cuts.SetMaxDCAToVertexZ(2.);
801
802 for(Int_t itrack = 0; itrack < ESDEvent->GetNumberOfTracks(); itrack++){
803 const AliESDtrack* esdtrack = ESDEvent->GetTrack(itrack);
804 Float_t b[2];
805 Float_t bCov[3];
806 esdtrack->GetImpactParameters(b,bCov);
807
808 if (bCov[0]<=0 || bCov[2]<=0) {
649e06c4 809 printf("AliCDMesonUtils - Estimated b resolution lower or equal zero!\n");
369562e9 810 bCov[0]=0;
811 bCov[2]=0;
812 }
813
814 Float_t dcaToVertexXY = b[0];
815 Float_t dcaToVertexZ = b[1];
816 if (esdtrack->Eta() > etacut) {
817 if (cuts.AcceptTrack(esdtrack)) nA++;
818 if (TPCGapDCAaSide) {
819 TPCGapDCAaSide->Fill(dcaToVertexXY, dcaToVertexZ);
820 }
821 }
822 else if (esdtrack->Eta() < -etacut) {
823 if (cuts.AcceptTrack(esdtrack)) nC++;
824 if (TPCGapDCAcSide) {
825 TPCGapDCAcSide->Fill(dcaToVertexXY, dcaToVertexZ);
826 }
827 }
828 }
829
830 const Bool_t tpcA = nA;
831 const Bool_t tpcC = nC;
832
833 return tpcA * AliCDMesonBase::kBitTPCA + tpcC * AliCDMesonBase::kBitTPCC;
834}
835
836
837//------------------------------------------------------------------------------
838Int_t AliCDMesonUtils::GetZDC(const AliESDEvent * ESDEvent)
839{
840 //
841 //GetZDC
842 //
843
844 const Int_t qa = ESDEvent->GetESDZDC()->GetESDQuality();
845 Bool_t zdcA = kFALSE, zdcC = kFALSE;
846 for(Int_t ii=0; ii<6; ii++){
847 if(qa & (1<<ii)){
848 if(ii<4) zdcA = kTRUE;
849 else zdcC = kTRUE;
850 }
851 }
852
853 return zdcA * AliCDMesonBase::kBitZDCA + zdcC * AliCDMesonBase::kBitZDCC;
854}
855
856
857//==============================================================================
858//------------------------------------------------------------------------------
859#ifdef STD_ALIROOT
860Int_t AliCDMesonUtils::FMDHitCombinations(const AliESDEvent* ESDEvent,
861 Int_t side)
862{
863 //
864 // copy of the FMDHitCombinations function originating from AliTriggerAnalysis
865 //
866 // side == 0 -> A side, side > 0 -> C side
867
868 // workaround for AliESDEvent::GetFMDData is not const!
869 const AliESDFMD* fmdData =
870 (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
871 if (!fmdData)
872 {
873 puts("AliESDFMD not available");
874 return -1;
875 }
876
877 Int_t detFrom = (side == 0) ? 1 : 3;
878 Int_t detTo = (side == 0) ? 2 : 3;
879
880 Int_t triggers = 0;
881 Float_t totalMult = 0;
882 for (UShort_t det=detFrom;det<=detTo;det++) {
883 Int_t nRings = (det == 1 ? 1 : 2);
884 for (UShort_t ir = 0; ir < nRings; ir++) {
885 Char_t ring = (ir == 0 ? 'I' : 'O');
886 UShort_t nsec = (ir == 0 ? 20 : 40);
887 UShort_t nstr = (ir == 0 ? 512 : 256);
888 for (UShort_t sec =0; sec < nsec; sec++) {
889 for (UShort_t strip = 0; strip < nstr; strip++) {
890 Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
891 if (mult == AliESDFMD::kInvalidMult) continue;
892
893 if (mult > 0.3) // fFMDLowCut
894 totalMult = totalMult + mult;
895 else {
896 if (totalMult > 0.5) // fFMDHitCut
897 triggers++;
898 }
899 }
900 }
901 }
902 }
903 return triggers;
904}
905#endif
906
907
908//==============================================================================
909//------------------------------------------------------------------------------
34e8fe1d 910void AliCDMesonUtils::SPDLoadGeom(Int_t run)
369562e9 911{
912 // method to get the gGeomanager
913 // it is called at the CreatedOutputObject stage
914 // to comply with the CAF environment
915
916 AliCDBManager *man = AliCDBManager::Instance();
917
918 TString cdbpath;
919 if (man->IsDefaultStorageSet()) {
920 const AliCDBStorage *dsto = man->GetDefaultStorage();
921 cdbpath = TString(dsto->GetBaseFolder());
922 //printf("default was set to: %s\n", cdbpath.Data());
923 }
924 else { //should not be used!
925 // man->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
926 // would be needed on grid
927 man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
928 cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
929 }
930
931 man->SetSpecificStorage("ITS/Align/Data",cdbpath);
932 man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
933 man->SetRun(run);
934
935 AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
936 if (!obj) {
649e06c4 937 printf("AliCDMesonUtils failed loading geometry object\n");
369562e9 938 return;
939 }
940 AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
649e06c4 941 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
369562e9 942}
943
944//------------------------------------------------------------------------------
34e8fe1d 945Bool_t AliCDMesonUtils::SPDLoc2Glo(Int_t id, const Double_t *loc,
649e06c4 946 Double_t *glo)
369562e9 947{
948 //
949 //SPDLoc2Glo, do not touch
950 //
951
952 static TGeoHMatrix mat;
953 Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
954 if (vid<0) {
649e06c4 955 printf("AliCDMesonUtils Did not find module with such ID %d\n",id);
369562e9 956 return kFALSE;
957 }
958 AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
959 mat.LocalToMaster(loc,glo);
960 return kTRUE;
961}
962
963
964//------------------------------------------------------------------------------
34e8fe1d 965Int_t AliCDMesonUtils::CheckChipEta(Int_t chipKey, TString scut,
369562e9 966 const Double_t vtxPos[],
967 TH2 *hitMapSPDinner, TH2 *hitMapSPDouter)
968{
969 //
970 //CheckChipEta
971 //
972
973 // retrieves the position in eta for a given chip and applies the cut
974 // results:
975 // 0 <= out of range
976 // -1 <= negative pseudo-rapidity position, in range (C-Side)
977 // 1 <= positive pseudo-rapidity position, in range (A-Side)
978 //
979 // scut: "[0.9" or "]0.9", only 3 digits for the value!!
980
981
982 const Bool_t kincl = (scut[0] == '[');
983 const TString cutval = scut(1,3);
984 const Double_t etacut = fabs(cutval.Atof());
985
986 //no eta cut, save time
987 if(kincl && etacut>=2)
988 return kTRUE;
989
990 Int_t etaside = 1;
991 //------------------------------- NOT TO TOUCH ------------------------>>
992 UInt_t module=999, offchip=999;
993 AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
994 UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
995 if(hs<2) offchip = 4 - offchip; // inversion in the inner layer...
996
997 const Int_t col[]={
998 hs<2? 0 : 31,
999 hs<2? 31 : 0,
1000 hs<2? 31 : 0,
1001 hs<2? 0 : 31};
1002 const Int_t aa[]={0, 0, 255, 255};
1003 const AliITSsegmentationSPD seg;
1004
1005 for(Int_t ic=0; ic<4; ic++){
1006 Float_t localchip[3]={0.,0.,0.};
1007 seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]);
1008 // local coordinate of the chip center
1009 //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
1010 const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
1011 Double_t glochip[3]={0.,0.,0.};
1012 if(!SPDLoc2Glo(module,local,glochip)){
1013 return kFALSE;
1014 }
1015
1016 //-------------------------------------------------------------------<<
1017
1018 const TVector3 pos(glochip[0]-vtxPos[0], glochip[1]-vtxPos[1],
1019 glochip[2]-vtxPos[2]);
1020 //pos.Print();
1021
1022 if (chipKey < 400) { // inner SPD layer
1023 if (hitMapSPDinner) {
1024 Double_t phi = pos.Phi(); // output in the range -Pi +Pi
1025 if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
1026 const Double_t eta = pos.Eta();
1027 hitMapSPDinner->Fill(eta, phi);
1028 }
1029 }
1030 else {
1031 if (hitMapSPDouter) { // outer SPD layer
1032 Double_t phi = pos.Phi(); // output in the range -Pi +Pi
1033 if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
1034 const Double_t eta = pos.Eta();
1035 hitMapSPDouter->Fill(eta, phi);
1036 }
1037 }
1038
1039 if( kincl && fabs(pos.Eta()) > etacut)
1040 return kFALSE;
1041
1042 if(!kincl){
1043 if(fabs(pos.Eta()) < etacut)
1044 return kFALSE;
1045 else if(pos.Eta()<0)
1046 etaside = -1;
1047 else
1048 etaside = 1;
1049 }
1050 }
1051
1052 return etaside;
1053}
1054
1055
1056//------------------------------------------------------------------------------
34e8fe1d 1057void AliCDMesonUtils::GetNFO(const AliESDEvent *ESDEvent, TString etacut,
369562e9 1058 Int_t ctr[], TH2 *hitMapSPDinner,
1059 TH2 *hitMapSPDouter)
1060{
1061 //
1062 // GetNFO
1063 //
1064 // analyzes the SPD fastOR for a given eta range and returns
1065 // an array with the number of hits in:
1066
1067 Int_t ninner=0; // inner layer
1068 Int_t nouter=0; // outer layer
1069 Int_t ipA = 0; // inner layer A side
1070 Int_t ipC = 0; // inner layer C side
1071 Int_t opA = 0; // outer layer A side
1072 Int_t opC = 0; // outer layer C side
1073
1074 const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
1075
1076 // position of the primary vertex
1077 Double_t tmp[3] = { 0., 0., 0. };
1078 ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
1079 Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
1080
1081
1082 for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
1083 if(mult->TestFastOrFiredChips(iChipKey)){
1084 // here you check if the FastOr bit is 1 or 0
1085 const Int_t iseta = CheckChipEta(iChipKey, etacut, vtxPos, hitMapSPDinner,
1086 hitMapSPDouter);
1087 if(iseta==0)
1088 continue;
1089
1090 if(iChipKey<400) {
1091 ninner++; // here you count the FastOr bits in the inner layer
1092 if(iseta>0)
1093 ipA ++;
1094 else
1095 ipC ++;
1096 }
1097 else {
1098 nouter++; // here you count the FastOr bits in the outer layer
1099 if(iseta>0)
1100 opA ++;
1101 else
1102 opC ++;
1103 }
1104 }
1105 }
1106
1107 ctr[kInnerPixel]= ninner;
1108 ctr[kOuterPixel]= nouter;
1109 ctr[kIPA]= ipA;
1110 ctr[kIPC]= ipC;
1111 ctr[kOPA]= opA;
1112 ctr[kOPC]= opC;
1113
1114 return;
1115}
1116
1117
1118//--------------------------------------------------------------------------
1119Int_t AliCDMesonUtils::GetFastORmultiplicity(const AliESDEvent* ESDEvent)
1120{
1121 // determine the number of fired fastOR chips in both layers within
1122 // -0.9 < eta < 0.9
1123 //
1124
1125 const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
1126
1127 // position of the primary vertex
1128 Double_t tmp[3] = { 0., 0., 0. };
1129 ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
1130 Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
1131
1132 Int_t multiplicity = 0;
1133
1134 for (Int_t iChipKey=0; iChipKey < 1200; iChipKey++) {
1135 if(mult->TestFastOrFiredChips(iChipKey)){
1136 // here you check if the FastOr bit is 1 or 0
1137 const Int_t iseta = CheckChipEta(iChipKey, "[0.9]", vtxPos, 0x0, 0x0);
1138 if(iseta==0)
1139 continue;
1140 else
1141 ++multiplicity;
1142 }
1143 }
1144
1145 return multiplicity;
1146}
1147
1148
1149//--------------------------------------------------------------------------
1150void AliCDMesonUtils::FillSPDtrkltMap(const AliVEvent* event,
1151 TH2 *hitMapSPDtrklt)
1152{
1153 //
1154 // fill eta phi map of SPD tracklets
1155 //
1156
1157 if (hitMapSPDtrklt) {
1158 if (TString(event->ClassName()) == "AliESDEvent") {
1159 const AliMultiplicity *mult = ((AliESDEvent*)event)->GetMultiplicity();
1160 for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
1161 iTracklet++) {
1162 Double_t eta = mult->GetEta(iTracklet);
1163 Double_t phi = mult->GetPhi(iTracklet);
1164 hitMapSPDtrklt->Fill(eta, phi);
1165 }
1166 }
1167 else if (TString(event->ClassName()) == "AliAODEvent") {
1168 const AliAODTracklets* mult = ((AliAODEvent*)event)->GetTracklets();
1169 for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
1170 iTracklet++) {
1171 Double_t eta = -TMath::Log(TMath::Tan(mult->GetTheta(iTracklet)/2.));
1172 Double_t phi = mult->GetPhi(iTracklet);
1173 hitMapSPDtrklt->Fill(eta, phi);
1174 }
1175 }
1176 }
1177}
1178
1179
1180//==========================================================================
1181Int_t AliCDMesonUtils::GetPID(AliPIDResponse *pid, const AliVTrack *trk,
34e8fe1d 1182 Int_t mode /* = 0 */)
369562e9 1183{
1184 // determines PID for ESDs and AODs
1185 //
1186 //
1187
1188 if (!pid) return AliCDMesonBase::kBinPIDUnknown; // no pid available
1189
1190 Double_t tpcpion = -999., tpckaon = -999., tpcproton = -999.,
1191 tpcelectron = -999.;
1192 Double_t tofpion = -999., tofkaon = -999., tofproton = -999.,
1193 tofelectron = -999.;
1194 Double_t itspion = -999., itskaon = -999., itsproton = -999.,
1195 itselectron = -999.;
1196
1197
1198 // check whether the track was pure ITS standalone track (has not left TPC)
1199 const Bool_t kits = !(trk->GetStatus() & AliESDtrack::kTPCout);
1200
1201 if (kits) { // do ITS pid
1202 const Double_t nin = 3.;
1203
1204 itspion = pid->NumberOfSigmasITS( trk, AliPID::kPion );
1205 itskaon = pid->NumberOfSigmasITS( trk, AliPID::kKaon );
1206 itsproton = pid->NumberOfSigmasITS( trk, AliPID::kProton );
1207 itselectron = pid->NumberOfSigmasITS( trk, AliPID::kElectron );
1208
1209 if(fabs(itspion) < nin) // inclusion only
1210 return AliCDMesonBase::kBinPion;
1211 else if(fabs(itskaon) < nin)
1212 return AliCDMesonBase::kBinKaon;
1213 else if(fabs(itsproton) < nin)
1214 return AliCDMesonBase::kBinProton;
1215 else if(fabs(itselectron) < nin)
1216 return AliCDMesonBase::kBinElectron;
1217 else{
1218 return AliCDMesonBase::kBinPIDUnknown;
1219 }
1220 }
1221
1222 tpcpion = pid->NumberOfSigmasTPC( trk, AliPID::kPion );
1223 tpckaon = pid->NumberOfSigmasTPC( trk, AliPID::kKaon );
1224 tpcproton = pid->NumberOfSigmasTPC( trk, AliPID::kProton );
1225 tpcelectron = pid->NumberOfSigmasTPC( trk, AliPID::kElectron );
1226
1227 // derive information, whether tof pid is available
1228 const Bool_t ka = !(trk->GetStatus() & AliESDtrack::kTOFmismatch);
1229 const Bool_t kb = (trk->GetStatus() & AliESDtrack::kTOFpid);
1230 const Bool_t ktof = ka && kb;
1231
1232 // retrieve TOF information if available
1233 if(ktof){
1234 tofpion = pid->NumberOfSigmasTOF(trk, AliPID::kPion);
1235 tofkaon = pid->NumberOfSigmasTOF(trk, AliPID::kKaon);
1236 tofproton = pid->NumberOfSigmasTOF(trk, AliPID::kProton);
1237 tofelectron = pid->NumberOfSigmasTOF(trk, AliPID::kElectron);
1238 }
1239
1240 if (mode == 0) {
1241 // 2010 cuts for PID
1242 const Double_t nin = 3.;
1243
1244 // inclusion + exclusion (either TPC or TOF)
1245 if((fabs(tpcpion) < nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
1246 && fabs(tpcelectron) > nin) ||
1247 (fabs(tofpion) < nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
1248 && fabs(tofelectron) > nin))
1249 return AliCDMesonBase::kBinPionE;
1250 else if((fabs(tpcpion) > nin && fabs(tpckaon) < nin && fabs(tpcproton) > nin
1251 && fabs(tpcelectron) > nin) ||
1252 (fabs(tofpion) > nin && fabs(tofkaon) < nin && fabs(tofproton) > nin
1253 && fabs(tofelectron) > nin))
1254 return AliCDMesonBase::kBinKaonE;
1255 else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) < nin
1256 && fabs(tpcelectron) > nin) ||
1257 (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) < nin
1258 && fabs(tofelectron) > nin))
1259 return AliCDMesonBase::kBinProtonE;
1260 else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
1261 && fabs(tpcelectron) < nin) ||
1262 (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
1263 && fabs(tofelectron) < nin))
1264 return AliCDMesonBase::kBinElectronE;
1265 else if(fabs(tpcpion) < nin && fabs(tofpion) < nin) // inclusion (TPC + TOF)
1266 return AliCDMesonBase::kBinPion;
1267 else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
1268 return AliCDMesonBase::kBinKaon;
1269 else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
1270 return AliCDMesonBase::kBinProton;
1271 else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
1272 return AliCDMesonBase::kBinElectron;
1273 else{
1274 return AliCDMesonBase::kBinPIDUnknown;
1275 }
1276 }
1277 else if (mode == 1) {
1278 // 2011 cuts for PID in LHC11f
1279 // TPC: [-3,5] sigma (pion)
1280 // TOF: 3 sigma for all,
1281 // only Pion is tuned!
1282 // ONLY INCLUSION CUTS NO EXCLUSION
1283 const Double_t nin = 3.;
1284
1285 if(tpcpion < 4. && tpcpion > -2. && fabs(tofpion) < -3. )
1286 return AliCDMesonBase::kBinPion;
1287 else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
1288 return AliCDMesonBase::kBinKaon;
1289 else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
1290 return AliCDMesonBase::kBinProton;
1291 else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
1292 return AliCDMesonBase::kBinElectron;
1293 else{
1294 return AliCDMesonBase::kBinPIDUnknown;
1295 }
1296 }
1297 return AliCDMesonBase::kBinPIDUnknown;
1298}
1299
1300
1301//==========================================================================
34e8fe1d 1302Int_t AliCDMesonUtils::GetPID(Int_t pdgCode)
369562e9 1303{
1304 //
1305 // determine particle type based on PDG code
1306 //
1307
1308 if (TMath::Abs(pdgCode) == 211) return AliCDMesonBase::kBinPionE;
1309 else if (TMath::Abs(pdgCode) == 321) return AliCDMesonBase::kBinKaonE;
1310 else if (TMath::Abs(pdgCode) == 2212) return AliCDMesonBase::kBinProtonE;
1311 else if (TMath::Abs(pdgCode) == 11) return AliCDMesonBase::kBinElectronE;
1312 else return AliCDMesonBase::kBinPIDUnknown;
1313}
1314
1315
1316//------------------------------------------------------------------------------
1317Int_t AliCDMesonUtils::CombinePID(const Int_t pid[])
1318{
1319 //
1320 // combine the PID result
1321 //
1322
1323 // determine return value
1324 if (pid[0] == pid[1]) { // same result for both tracks
1325 return pid[0];
1326 }
1327 // one track identified with exclusion the other only without
1328 else if ((pid[0] == AliCDMesonBase::kBinPionE &&
1329 pid[1] == AliCDMesonBase::kBinPion) ||
1330 (pid[1] == AliCDMesonBase::kBinPionE &&
1331 pid[0] == AliCDMesonBase::kBinPion)) {
1332 return AliCDMesonBase::kBinPion;
1333 }
1334 else if ((pid[0] == AliCDMesonBase::kBinKaonE &&
1335 pid[1] == AliCDMesonBase::kBinKaon) ||
1336 (pid[1] == AliCDMesonBase::kBinKaonE &&
1337 pid[0] == AliCDMesonBase::kBinKaon)) {
1338 return AliCDMesonBase::kBinKaon;
1339 }
1340 else if ((pid[0] == AliCDMesonBase::kBinProtonE &&
1341 pid[1] == AliCDMesonBase::kBinProton) ||
1342 (pid[1] == AliCDMesonBase::kBinProtonE &&
1343 pid[0] == AliCDMesonBase::kBinProton)) {
1344 return AliCDMesonBase::kBinProton;
1345 }
1346 else if ((pid[0] == AliCDMesonBase::kBinElectronE &&
1347 pid[1] == AliCDMesonBase::kBinElectron) ||
1348 (pid[1] == AliCDMesonBase::kBinElectronE &&
1349 pid[0] == AliCDMesonBase::kBinElectron)) {
1350 return AliCDMesonBase::kBinElectron;
1351 }
1352 // one track identified and one not
1353 else if (((pid[0] == AliCDMesonBase::kBinPionE ||
1354 pid[0] == AliCDMesonBase::kBinPion) &&
1355 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1356 (pid[1] == AliCDMesonBase::kBinPion &&
1357 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1358 return AliCDMesonBase::kBinSinglePion;
1359 }
1360 else if (((pid[0] == AliCDMesonBase::kBinKaonE ||
1361 pid[0] == AliCDMesonBase::kBinKaon)&&
1362 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1363 (pid[1] == AliCDMesonBase::kBinKaon &&
1364 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1365 return AliCDMesonBase::kBinSingleKaon;
1366 }
1367 else if (((pid[0] == AliCDMesonBase::kBinProtonE ||
1368 pid[0] == AliCDMesonBase::kBinProton) &&
1369 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1370 (pid[1] == AliCDMesonBase::kBinProton &&
1371 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1372 return AliCDMesonBase::kBinSingleProton;
1373 }
1374 else if (((pid[0] == AliCDMesonBase::kBinElectronE ||
1375 pid[0] == AliCDMesonBase::kBinElectron) &&
1376 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1377 (pid[1] == AliCDMesonBase::kBinElectron &&
1378 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1379 return AliCDMesonBase::kBinSingleElectron;
1380 }
1381 else
1382 return AliCDMesonBase::kBinPIDUnknown;
1383}
1384
1385
1386//------------------------------------------------------------------------------
1387TLorentzVector AliCDMesonUtils::GetKinematics(const Double_t *pa,
1388 const Double_t *pb,
34e8fe1d 1389 Double_t ma,
1390 Double_t mb, Double_t& cts)
369562e9 1391{
1392 //
1393 //get kinematics, cts = cos(theta_{#star})
1394 //
1395
1396 TLorentzVector va, vb;
1397 va.SetXYZM(pa[0], pa[1], pa[2], ma);
1398 vb.SetXYZM(pb[0], pb[1], pb[2], mb);
1399 const TLorentzVector sumv = va+vb;
1400
1401 const TVector3 bv = -sumv.BoostVector();
1402
1403 va.Boost(bv);
1404 vb.Boost(bv);
1405
1406 // 3-vectors in the restframe of the mother particle
1407 const TVector3 pra = va.Vect();
1408 const TVector3 prb = vb.Vect();
1409 const TVector3 diff = pra - prb; // their difference
1410
1411 cts = (diff.Mag2()-pra.Mag2()-prb.Mag2()) / (-2.*pra.Mag()*prb.Mag());
1412 // cosine theta star, calculated according to the law-of-cosines
1413
1414 return sumv;
1415}
1416
1417
1418//------------------------------------------------------------------------------
1419Double_t AliCDMesonUtils::GetOA(const Double_t *pa, const Double_t *pb)
1420{
1421 //
1422 //cosOpeningAngle
1423 //
1424
1425 TVector3 va, vb;
1426 va.SetXYZ(pa[0], pa[1], pa[2]);
1427 vb.SetXYZ(pb[0], pb[1], pb[2]);
1428
1429 const TVector3 ua = va.Unit();
1430 const TVector3 ub = vb.Unit();
1431
1432 return ua.Dot(ub);
1433}