Fixes for wrong use of const causing PW.CAST_TO_QUALIFIED_TYPE defect in Coverity
[u/mrichter/AliRoot.git] / PWGUD / DIFFRACTIVE / xsAndTwoProng / AliCDMesonUtils.cxx
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 //------------------------------------------------------------------------------
66 void AliCDMesonUtils::GetMassPtCtsOA(Int_t pid,
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 //------------------------------------------------------------------------------
105 void AliCDMesonUtils::GetPWAinfo(Int_t pid, const AliVTrack *trks[],
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 //------------------------------------------------------------------------------
173 Int_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 //------------------------------------------------------------------------------
207 void 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 //------------------------------------------------------------------------------
225 Int_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 //------------------------------------------------------------------------------
243 Int_t AliCDMesonUtils::GetCombPID(AliPIDResponse* pid, const AliVTrack* trks[],
244                                   Int_t mode, TH2* comb2trkPID /*= 0x0 */)
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 //------------------------------------------------------------------------------
270 Int_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 //------------------------------------------------------------------------------
293 void 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 //------------------------------------------------------------------------------
323 Bool_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         /*
335         //printf("AliCDMesonUtils active triggers: %s\n",
336         //       ESDEvent->GetHeader()->GetActiveTriggerInputs().Data());
337         //printf("AliCDMesonUtils fired triggers: %s\n",
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 //------------------------------------------------------------------------------
415 Bool_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 //------------------------------------------------------------------------------
452 void AliCDMesonUtils::DoVZEROStudy(const AliESDEvent *ESDEvent,
453                                    TObjArray* hists, Int_t run)
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");
491         if (!ent1) {
492                 printf("AliCDMesonUtils failed loading VZERO trigger entry from OCDB\n");
493                 return;
494         }
495         AliVZEROTriggerData *trigData = (AliVZEROTriggerData*)ent1->GetObject();
496         if (!trigData) {
497                 printf("AliCDMesonUtils failed loading VZERO trigger data from OCDB\n");
498                 return;
499         }
500
501         AliCDBEntry *ent2 = man->Get("VZERO/Calib/Data");
502         if (!ent2) {
503                 printf("AliCDMesonUtils failed loading VZERO calib entry from OCDB\n");
504                 return;
505         }
506         AliVZEROCalibData *calData = (AliVZEROCalibData*)ent2->GetObject();
507         if (!calData) {
508                 printf("AliCDMesonUtils failed loading VZERO calibration data from OCDB\n");
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 //------------------------------------------------------------------------------
536 Int_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 //------------------------------------------------------------------------------
566 void 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 //------------------------------------------------------------------------------
594 void 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 //------------------------------------------------------------------------------
640 void 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 //------------------------------------------------------------------------------
659 Int_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 //------------------------------------------------------------------------------
679 Int_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 //------------------------------------------------------------------------------
760 Int_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 //------------------------------------------------------------------------------
787 Int_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) {
809                         printf("AliCDMesonUtils - Estimated b resolution lower or equal zero!\n");
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 //------------------------------------------------------------------------------
838 Int_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
860 Int_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 //------------------------------------------------------------------------------
910 void AliCDMesonUtils::SPDLoadGeom(Int_t run)
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) {
937                 printf("AliCDMesonUtils failed loading geometry object\n");
938                 return;
939         }
940         AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
941         AliGeomManager::ApplyAlignObjsFromCDB("ITS");
942 }
943
944 //------------------------------------------------------------------------------
945 Bool_t AliCDMesonUtils::SPDLoc2Glo(Int_t id, const Double_t *loc,
946                                    Double_t *glo)
947 {
948         //
949         //SPDLoc2Glo, do not touch
950         //
951
952         static TGeoHMatrix mat;
953         Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
954         if (vid<0) {
955                 printf("AliCDMesonUtils Did not find module with such ID %d\n",id);
956                 return kFALSE;
957         }
958         AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
959         mat.LocalToMaster(loc,glo);
960         return kTRUE;
961 }
962
963
964 //------------------------------------------------------------------------------
965 Int_t AliCDMesonUtils::CheckChipEta(Int_t chipKey, TString scut,
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 //------------------------------------------------------------------------------
1057 void AliCDMesonUtils::GetNFO(const AliESDEvent *ESDEvent, TString etacut,
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 //--------------------------------------------------------------------------
1119 Int_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 //--------------------------------------------------------------------------
1150 void 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 //==========================================================================
1181 Int_t AliCDMesonUtils::GetPID(AliPIDResponse *pid, const AliVTrack *trk,
1182                               Int_t mode /* = 0 */)
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 //==========================================================================
1302 Int_t AliCDMesonUtils::GetPID(Int_t pdgCode)
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 //------------------------------------------------------------------------------
1317 Int_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 //------------------------------------------------------------------------------
1387 TLorentzVector AliCDMesonUtils::GetKinematics(const Double_t *pa,
1388                                               const Double_t *pb,
1389                                               Double_t ma,
1390                                               Double_t mb, Double_t& cts)
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 //------------------------------------------------------------------------------
1419 Double_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 }