Changed DimuonCombinator into AliDimuCombinator to follow ALICE coding
[u/mrichter/AliRoot.git] / EVGEN / AliDimuCombinator.cxx
1 //
2 //
3 //
4 //
5 #include "AliDimuCombinator.h" 
6 #include "AliRun.h" 
7 #include "TRandom.h" 
8 //
9 ClassImp(AliDimuCombinator)
10 //
11 //                       Iterators
12 // 
13  GParticle* AliDimuCombinator::FirstMuon()
14      {
15          fimuon1=fimin1;
16          fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
17          while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
18              fimuon1++;
19              if (fimuon1 >= fimax1) {fmuon1=0; break;}
20              fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
21          }
22          return fmuon1;
23      }
24              
25  GParticle* AliDimuCombinator::FirstMuonSelected()
26      {
27          GParticle * muon=FirstMuon();
28          while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
29          return muon;
30      }
31              
32
33  GParticle* AliDimuCombinator::NextMuon()
34      {
35          fimuon1++;
36          if (fimuon1>=fNParticle) {fmuon1 = 0; return fmuon1;}
37          
38          fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
39          while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
40              fimuon1++;
41              if (fimuon1>=fimax1) {fmuon1 = 0; break;}
42              fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
43          }
44          return fmuon1;
45      }
46
47 GParticle* AliDimuCombinator::NextMuonSelected()
48      {
49          GParticle * muon=NextMuon();
50          while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
51          return muon;
52      }
53
54
55  void AliDimuCombinator::FirstPartner()
56      {
57          if (fimin1==fimin2) {
58              fimuon2=fimuon1+1;
59          } else {
60              fimuon2=fimin2;
61          }
62          if (fimuon2 >= fimax2) {fmuon2=0; return;}
63          fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
64          while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
65              fimuon2++;
66              if (fimuon2 >= fimax2) {fmuon2=0; break;}
67              fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
68          }
69      }
70 void AliDimuCombinator::FirstPartnerSelected()
71 {
72          FirstPartner();
73          while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
74 }
75
76
77  void AliDimuCombinator::NextPartner()
78      {
79          fimuon2++;
80          if (fimuon2>=fimax2) {fmuon2 = 0; return;}
81
82          
83          fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
84
85          while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
86              fimuon2++;
87              if (fimuon2>=fimax2) {fmuon2 = 0; break;}
88              fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
89          }
90
91      }
92
93 void AliDimuCombinator::NextPartnerSelected()
94 {
95          NextPartner();
96          while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
97 }
98
99
100  GParticle*  AliDimuCombinator::Partner()
101      {
102          return fmuon2;
103      }
104
105 void AliDimuCombinator::FirstMuonPair(GParticle* & muon1, GParticle* & muon2)
106      {
107          FirstMuon();
108          FirstPartner();
109          muon1=fmuon1;
110          muon2=fmuon2;   
111      }
112 void AliDimuCombinator::NextMuonPair(GParticle* & muon1, GParticle* & muon2)
113      {
114          NextPartner();
115          if (!Partner()) {
116              NextMuon();
117              FirstPartner();
118          }
119          muon1=fmuon1;
120          muon2=fmuon2;   
121      }
122 void AliDimuCombinator::FirstMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
123      {
124          FirstMuonSelected();
125          FirstPartnerSelected();
126          muon1=fmuon1;
127          muon2=fmuon2;   
128      }
129 void AliDimuCombinator::NextMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
130      {
131          NextPartnerSelected();
132          if (!Partner()) {
133              NextMuonSelected();
134              FirstPartnerSelected();
135          }
136          muon1=fmuon1;
137          muon2=fmuon2;   
138      }
139 void AliDimuCombinator::ResetRange()
140 {
141     fimin1=fimin2=0;
142     fimax1=fimax2=fNParticle;
143 }
144
145 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
146 {
147     fimin1=from;
148     fimax1=to;
149     if (fimax1 > fNParticle) fimax1=fNParticle;
150 }
151
152 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
153 {
154     fimin2=from;
155     fimax2=to;
156     if (fimax2 > fNParticle) fimax2=fNParticle;
157 }
158 //
159 //                       Selection
160 //
161
162 Bool_t AliDimuCombinator::Selected(GParticle* part)
163 {
164 // 
165 //
166     if (part==0) {return 0;}
167     
168     if (part->GetPT() > fPtMin && part->GetEta()>fEtaMin && part->GetEta()<fEtaMax) {
169         return 1;
170     } else {
171         return 0;
172     }
173     
174     
175 }
176
177 Bool_t AliDimuCombinator::Selected(GParticle* part1, GParticle* part2)
178 {
179      return Selected(part1)*Selected(part2);
180 }
181 //
182 //                       Kinematics
183 //
184 Float_t AliDimuCombinator::Mass(GParticle* part1, GParticle* part2)
185 {
186     Float_t px,py,pz,e;
187     px=part1->GetPx()+part2->GetPx();
188     py=part1->GetPy()+part2->GetPy();
189     pz=part1->GetPz()+part2->GetPz();    
190     e =part1->GetEnergy()+part2->GetEnergy();
191     Float_t p=px*px+py*py+pz*pz;
192     if (e*e < p) {
193         return -1; 
194     } else {
195         return TMath::Sqrt(e*e-p);
196     }
197 }
198
199 Float_t AliDimuCombinator::PT(GParticle* part1, GParticle* part2)
200 {
201     Float_t px,py;
202     px=part1->GetPx()+part2->GetPx();
203     py=part1->GetPy()+part2->GetPy();
204     return TMath::Sqrt(px*px+py*py);
205 }
206
207 Float_t AliDimuCombinator::Pz(GParticle* part1, GParticle* part2)
208 {
209     return part1->GetPz()+part2->GetPz();
210 }
211
212 Float_t AliDimuCombinator::Y(GParticle* part1, GParticle* part2)
213 {
214     Float_t pz,e;
215     pz=part1->GetPz()+part2->GetPz();
216     e =part1->GetEnergy()+part2->GetEnergy();
217     return 0.5*TMath::Log((e+pz)/(e-pz));
218 }
219 //                  Response
220 //
221 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
222 {
223     value+=gRandom->Gaus(0, width);
224 }
225 //              Weighting
226 // 
227
228 Float_t AliDimuCombinator::Decay_Prob(GParticle* part)
229 {
230     Float_t d, h, theta, CTau;
231     GParticle* parent = Parent(part);
232     Int_t ipar=Type(parent);
233     if (ipar==8 || ipar==9) {
234         CTau=780.4;
235     } else if (ipar==11 || ipar==12) {
236         CTau=370.9;
237     } else {
238         CTau=0;
239     }
240     
241     
242     Float_t GammaBeta=(parent->GetMomentum())/(parent->GetMass());
243 //
244 // this part is still very ALICE muon-arm specific
245 //
246     theta=parent->GetTheta();
247     h=90*TMath::Tan(theta);
248     
249     if (h<4) {
250         d=4/TMath::Sin(theta);
251     } else {
252         d=90/TMath::Cos(theta);
253     }
254     
255     if (CTau > 0) {
256         return 1-TMath::Exp(-d/CTau/GammaBeta);
257     } else {
258         return 1;
259     }
260 }
261
262 Float_t AliDimuCombinator::Weight(GParticle* part1, GParticle* part2)
263 {
264     Float_t wgt=(part1->GetWgt())*(part2->GetWgt());
265     
266     if (Correlated(part1, part2)) {
267         return wgt/(Parent(part1)->GetWgt())*fRate1;
268     } else {
269         return wgt*fRate1*fRate2;
270     }
271
272
273
274 Float_t AliDimuCombinator::Weight(GParticle* part)
275 {
276     return (part->GetWgt())*(Parent(part)->GetWgt())*fRate1;
277 }
278 Bool_t  AliDimuCombinator::Correlated(GParticle* part1, GParticle* part2)
279 {
280     if (Origin(part1) == Origin(part2)) {
281         return kTRUE;
282     } else {
283         return kFALSE;
284     }
285 }
286 GParticle* AliDimuCombinator::Parent(GParticle* part)
287 {
288     return (GParticle*) (fPartArray->UncheckedAt(part->GetParent()));
289 }
290
291 Int_t AliDimuCombinator::Origin(GParticle* part)
292 {
293     Int_t iparent= part->GetParent();
294     if (iparent < 0) return iparent;
295     Int_t ip;
296     while(1) {
297         ip=((GParticle*) fPartArray->UncheckedAt(iparent))->GetParent();
298         if (ip < 0) {
299             break;
300         } else {
301             iparent=ip;
302         }
303     }
304     return iparent;
305 }
306