bbb023706e4b50823c1ef998d85b8f37bfd2033c
[u/mrichter/AliRoot.git] / PWG4 / TwoPartCorr / AutoCorr.C
1 // $Id$
2
3 #include "AutoCorr.h"
4
5 ClassImp(AutoCorr)
6
7 Int_t AutoCorr::InitEventPools(Int_t depth, 
8                                Int_t nMultBins, Double_t multbin[], 
9                                Int_t nZvtxBins, Double_t zvtxbin[],
10                                Double_t ptMin, Double_t ptMax)
11 {
12   // First assign AutoCorr members
13   fNMultBins = nMultBins;
14   fNZvtxBins = nZvtxBins;
15
16   for (Int_t iM=0; iM<nMultBins; iM++) {
17     std::vector<EventPool*> evp;
18     for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
19       evp.push_back(new EventPool(depth, 
20                                   multbin[iM], multbin[iM+1], 
21                                   zvtxbin[iZ], zvtxbin[iZ+1],
22                                   ptMin, ptMax ));
23     }
24     fEvPool.push_back(evp);
25   }
26   
27   for (Int_t iM=0; iM<nMultBins; iM++) {
28     for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
29       fEvPool.at(iM).at(iZ)->SetMultBinIndex(iM);
30       fEvPool.at(iM).at(iZ)->SetZvtxBinIndex(iZ);
31     }
32   }
33   
34   bool print_this = false;
35   if (print_this) {
36     cout << "fEvPool outer size: " << fEvPool.size() << endl;
37     for (Int_t iM=0; iM<nMultBins; iM++) {
38       for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
39         if(fEvPool.at(iM).at(iZ)) {
40           cout << "multiplicity bin: " << iM;
41           cout << ", z-vertex bin: " << iZ;
42           fEvPool.at(iM).at(iZ)->PrintInfo();
43         }
44       }
45     }
46   }
47   
48   return 0;
49 }
50
51 EventPool* AutoCorr::GetEventPool(Int_t iMult, Int_t iZvtx) const
52 {
53   if (iMult < 0 || iMult >= fNMultBins) return 0x0;
54   if (iZvtx < 0 || iZvtx >= fNZvtxBins) return 0x0;
55   return fEvPool.at(iMult).at(iZvtx);
56 }
57
58 Int_t AutoCorr::UpdatePools(Int_t iEvent, const MyHeader* ev, TClonesArray* trk)
59 {
60   for (Int_t iM=0; iM<fNMultBins; iM++) {
61     for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
62       fEvPool.at(iM).at(iZ)->UpdatePool(iEvent, ev, trk);
63     }
64   }  
65   return 0;
66 }
67
68 Double_t AutoCorr::DeltaPhi(const MyPart &t1, const MyPart &t2,
69                             Double_t rangeMin, Double_t rangeMax) const
70 {
71   Double_t dphi = -999;
72   Double_t pi = TMath::Pi();
73   Double_t phia = t1.fPhi;  
74   Double_t phib = t2.fPhi;  
75   
76   if (phia < 0)         phia += 2*pi;
77   else if (phia > 2*pi) phia -= 2*pi;
78   if (phib < 0)         phib += 2*pi;
79   else if (phib > 2*pi) phib -= 2*pi;
80   dphi = phib - phia;
81   if (dphi < rangeMin)      dphi += 2*pi;
82   else if (dphi > rangeMax) dphi -= 2*pi;
83   
84   return dphi;
85 }
86
87 Double_t AutoCorr::DeltaEta(const MyPart &t1, const MyPart &t2) const
88 {
89   return t1.fEta - t2.fEta;
90 }
91
92 Bool_t AutoCorr::InBounds(Double_t val, Double_t min, Double_t max) const
93 {
94   if (val<min)
95     return 0;
96   if (val>max)
97     return 0;
98   return 1;
99 }
100
101 Bool_t AutoCorr::InBounds(Int_t val, Int_t min, Int_t max) const
102 {
103   if (val<min)
104     return 0;
105   if (val>max)
106     return 0;
107   return 1;
108 }
109
110 Bool_t AutoCorr::IsEventOk(const MyHeader &ev, Int_t minVc, 
111                            Int_t maxNTracklets, Double_t zMin, Double_t zMax,
112                            Int_t trbit) const
113 {
114   Bool_t VcOk = ev.fVc >= minVc;
115   Bool_t NTrackletsOK = ev.fNTracklets <= maxNTracklets;
116   Bool_t zOk = InBounds(ev.fVz, zMin, zMax);
117   Bool_t trOk = 1;
118   if (trbit) 
119     trOk = (ev.fTrClassMask>>trbit&1);
120   return (!ev.fIsPileupSPD && VcOk && NTrackletsOK && zOk && trOk);
121 }
122
123 Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax) const
124 {
125   return InBounds(t.fEta, etaMin, etaMax);    
126 }
127
128 Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax,
129                            Double_t ptMin, Double_t ptMax) const
130 {
131   Bool_t etaOk = InBounds(t.fEta, etaMin, etaMax);
132   Bool_t ptOk  = InBounds(t.fPt, ptMin, ptMax);    
133   return  etaOk && ptOk;
134 }
135
136 Bool_t AutoCorr::IsPairOk(const MyPart &t1, const MyPart &t2) const
137 {
138   Double_t deta = DeltaEta(t1, t2);
139   Double_t dphi = DeltaPhi(t1, t2);
140   Double_t dpmax = 0.03;
141   Double_t demax = 0.03;
142   Double_t dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
143   return (dr > 1);
144 }
145
146 Bool_t AutoCorr::IsMixedPairOk(const MyPart &t1, const MyPart &t2) const
147 {
148   Double_t deta = DeltaEta(t1, t2);
149   Double_t dphi = DeltaPhi(t1, t2);
150   Double_t dpmax = 0.03;
151   Double_t demax = 0.03;
152   Double_t dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
153   return (dr > 1);
154 }
155
156 Double_t AutoCorr::InvMass(const MyPart &p1, const MyPart &p2) const
157 {
158   Double_t px1 = p1.Px();
159   Double_t py1 = p1.Py();
160   Double_t pz1 = p1.Pz();
161   Double_t px2 = p2.Px();
162   Double_t py2 = p2.Py();
163   Double_t pz2 = p2.Pz();
164   Double_t pm1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
165   Double_t pm2 = TMath::Sqrt(px2*px2+py2*py2+pz1*pz2);
166   Double_t p12 = px1*px2+py1*py2+pz1*pz2;
167   Double_t m = TMath::Sqrt(pm1*pm2-p12);
168   return m;
169 }