986305bd3ec1f8a2e1bc00e0dd9f092f82828771
[u/mrichter/AliRoot.git] / PWG1 / AliMCInfo.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
17 ///////////////////////////////////////////////////////////////////////////
18 /*
19
20 Origin: marian.ivanov@cern.ch
21 Container classes with MC infomation
22
23 The AliMCInfo contains the information about the particles properties 
24 during transportation throuch ALICE Detector
25
26 The base Information :
27 TParticle - fParticle             -  properties of the particle at creation point
28 AliTrackReference - fXXXRefernces -  TClonesArray of refernces in differnt detectors
29 fNXXXRef                          -  number of the track refernces in differnt detectors
30 AliTPCdigitRow    - fTPCRow       -  the map of the hitted rows - (will be repalced by TBits)
31 fRowsWith*                        - number of rows hitted by particle
32 fMCtracks         -               - number of turn over of the track inside of the TPC
33
34 ++++
35 some additional information usable for tree draw - TO SPEED UP tree queries
36 IMPORTANT FOR PROOF FAST PROTOTYPING ANALYSIS 
37                                       
38
39
40
41 */
42
43 #if !defined(__CINT__) || defined(__MAKECINT__)
44 #include <stdio.h>
45 //ROOT includes
46 #include "Rtypes.h"
47 #include "TClonesArray.h"
48 //ALIROOT includes
49 #include "AliTrackReference.h"
50 #include "AliMCInfo.h" 
51 #endif
52
53 //
54 // 
55
56 ClassImp(AliTPCdigitRow)
57 ClassImp(AliMCInfo)
58
59
60
61 ////////////////////////////////////////////////////////////////////////
62 AliMCInfo::AliMCInfo():
63   fTrackRef(),
64   fTrackRefOut(),
65   fTRdecay(),
66   fPrimPart(0),
67   fParticle(),
68   fMass(0),
69   fCharge(0),
70   fLabel(0),
71   fEventNr(),
72   fMCtracks(),
73   fPdg(0),
74   fTPCdecay(0),
75   fRowsWithDigitsInn(0),
76   fRowsWithDigits(0),
77   fRowsTrackLength(0),
78   fTPCtrackLength(-1),
79   fPrim(0),
80   fTPCRow(), 
81   fNTPCRef(0),                    // tpc references counter
82   fNITSRef(0),                    // ITS references counter
83   fNTRDRef(0),                    // TRD references counter
84   fNTOFRef(0),                    // TOF references counter
85   fTPCReferences(0),
86   fITSReferences(0),
87   fTRDReferences(0),
88   fTOFReferences(0)
89 {
90   //
91   // Default constructor
92   //
93   fTPCReferences  = new TClonesArray("AliTrackReference",10);
94   fITSReferences  = new TClonesArray("AliTrackReference",10);
95   fTRDReferences  = new TClonesArray("AliTrackReference",10);
96   fTOFReferences  = new TClonesArray("AliTrackReference",10);
97   fTRdecay.SetTrack(-1);
98   fCharge = 0;
99 }
100
101 AliMCInfo::AliMCInfo(const AliMCInfo& info):
102   TObject(),
103   fTrackRef(info.fTrackRef),
104   fTrackRefOut(info.fTrackRefOut),
105   fTRdecay(info.GetTRdecay()),
106   fPrimPart(info.fPrimPart),
107   fParticle(info.fParticle),
108   fMass(info.GetMass()),
109   fCharge(info.fCharge),
110   fLabel(info.fLabel),
111   fEventNr(info.fEventNr),
112   fMCtracks(info.fMCtracks),
113   fPdg(info.fPdg),
114   fTPCdecay(info.fTPCdecay),
115   fRowsWithDigitsInn(info.fRowsWithDigitsInn),
116   fRowsWithDigits(info.fRowsWithDigits),
117   fRowsTrackLength(info.fRowsTrackLength),
118   fTPCtrackLength(info.fTPCtrackLength),
119   fPrim(info.fPrim),
120   fTPCRow(info.fTPCRow), 
121   fNTPCRef(info.fNTPCRef),                    // tpc references counter
122   fNITSRef(info.fNITSRef),                    // ITS references counter
123   fNTRDRef(info.fNTRDRef),                    // TRD references counter
124   fNTOFRef(info.fNTOFRef),                    // TOF references counter
125   fTPCReferences(0),
126   fITSReferences(0),
127   fTRDReferences(0),
128   fTOFReferences(0)
129 {
130   //
131   // copy constructor
132   //
133   fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
134   fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
135   fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
136   fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
137 }
138
139
140 AliMCInfo& AliMCInfo::operator=(const AliMCInfo& info) { 
141   //
142   // Assignment operator
143   //
144   this->~AliMCInfo();
145   new (this) AliMCInfo(info);
146   return *this;
147 }
148
149
150 AliMCInfo::~AliMCInfo()
151 {
152   //
153   // Destructor of the class
154   //
155   if (fTPCReferences) {
156     delete fTPCReferences;
157   }
158   if (fITSReferences){
159     delete fITSReferences;
160   }
161   if (fTRDReferences){    
162     delete fTRDReferences;  
163   }
164   if (fTOFReferences){    
165     delete fTOFReferences;  
166   }
167
168 }
169
170
171
172 void AliMCInfo::Update()
173 {
174   //
175   // Update MC info
176   // Calculates some derived variables
177   //
178   //
179   fMCtracks =1;
180   //
181   // decay info
182   fTPCdecay=kFALSE;
183   if (fTRdecay.GetTrack()>0){
184     fDecayCoord[0] = fTRdecay.X();
185     fDecayCoord[1] = fTRdecay.Y();
186     fDecayCoord[2] = fTRdecay.Z();
187     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
188       fTPCdecay=kTRUE;     
189     }
190     else{
191       fDecayCoord[0] = 0;
192       fDecayCoord[1] = 0;
193       fDecayCoord[2] = 0;
194     }
195   }
196   //
197   //
198   //digits information update
199   fRowsWithDigits    = fTPCRow.RowsOn();    
200   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
201   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
202   //
203   //
204   // calculate primary ionization per cm  
205   if (fParticle.GetPDG()){
206     fMass = fParticle.GetMass();  
207     fCharge = fParticle.GetPDG()->Charge();
208     if (fTPCReferences->GetEntriesFast()>0){
209       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
210     }
211     if (fMass>0){
212       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
213                               fTrackRef.Py()*fTrackRef.Py()+
214                               fTrackRef.Pz()*fTrackRef.Pz());    
215       if (p>0.001){
216         Float_t betagama = p /fMass;
217         fPrim = TPCBetheBloch(betagama);
218       }else fPrim=0;
219     }
220   }else{
221     fMass =0;
222     fPrim =0;
223   }  
224 }
225
226
227
228
229 ////////////////////////////////////////////////////////////////////////
230 AliTPCdigitRow::AliTPCdigitRow()
231 {
232   Reset();
233 }
234 ////////////////////////////////////////////////////////////////////////
235 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
236 {
237   for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
238   return (*this);
239 }
240 ////////////////////////////////////////////////////////////////////////
241 void AliTPCdigitRow::SetRow(Int_t row) 
242 {
243   //
244   // set bit mask for given row
245   //
246   if (row >= 8*32) {
247     //    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
248     return;
249   }
250   Int_t iC = row/8;
251   Int_t iB = row%8;
252   SETBIT(fDig[iC],iB);
253 }
254
255 ////////////////////////////////////////////////////////////////////////
256 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
257 {
258 //
259 // return kTRUE if row is on
260 //
261   Int_t iC = row/8;
262   Int_t iB = row%8;
263   return TESTBIT(fDig[iC],iB);
264 }
265 ////////////////////////////////////////////////////////////////////////
266 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
267 {
268 //
269 // returns number of rows with a digit  
270 // count only rows less equal row number upto
271 //
272   Int_t total = 0;
273   for (Int_t i = 0; i<32; i++) {
274     for (Int_t j = 0; j < 8; j++) {
275       if (i*8+j > upto) return total;
276       if (TESTBIT(fDig[i],j))  total++;
277     }
278   }
279   return total;
280 }
281 ////////////////////////////////////////////////////////////////////////
282 void AliTPCdigitRow::Reset()
283 {
284 //
285 // resets all rows to zero
286 //
287   for (Int_t i = 0; i<32; i++) {
288     fDig[i] <<= 8;
289   }
290 }
291 ////////////////////////////////////////////////////////////////////////
292 Int_t AliTPCdigitRow::Last() const
293 {
294 //
295 // returns the last row number with a digit
296 // returns -1 if now digits 
297 //
298   for (Int_t i = 32-1; i>=0; i--) {
299     for (Int_t j = 7; j >= 0; j--) {
300       if TESTBIT(fDig[i],j) return i*8+j;
301     }
302   }
303   return -1;
304 }
305 ////////////////////////////////////////////////////////////////////////
306 Int_t AliTPCdigitRow::First() const
307 {
308 //
309 // returns the first row number with a digit
310 // returns -1 if now digits 
311 //
312   for (Int_t i = 0; i<32; i++) {
313     for (Int_t j = 0; j < 8; j++) {
314       if (TESTBIT(fDig[i],j)) return i*8+j;
315     }
316   }
317   return -1;
318 }
319
320
321 //_____________________________________________________________________________
322 Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
323 {
324   //
325   // Bethe-Bloch energy loss formula
326   //
327   const Double_t kp1=0.76176e-1;
328   const Double_t kp2=10.632;
329   const Double_t kp3=0.13279e-4;
330   const Double_t kp4=1.8631;
331   const Double_t kp5=1.9479;
332
333   Double_t dbg = (Double_t) bg;
334
335   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
336
337   Double_t aa = TMath::Power(beta,kp4);
338   Double_t bb = TMath::Power(1./dbg,kp5);
339
340   bb=TMath::Log(kp3+bb);
341   
342   return ((Float_t)((kp2-aa-bb)*kp1/aa));
343 }
344
345
346 void AliMCInfo::CalcTPCrows(TClonesArray * runArrayTR){
347   //
348   // Calculates the numebr of the track references for detectors
349   // In case of changing direction - curling tracks - the counter is not increasing
350   //
351   // Rough calculation 
352   // of the first and last point in the TPC  
353   //
354   fNTPCRef = 0;
355   fNITSRef = 0;
356   fNTRDRef = 0;
357   fNTOFRef = 0;  
358   Float_t tpcminRadius=250;
359   Float_t tpcmaxRadius=80;
360   Float_t dir=0;
361   Int_t nover=0;
362   
363   for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
364     //
365     AliTrackReference *ref = (AliTrackReference*)runArrayTR->At(iTrackRef);  
366     Float_t newdirection = (ref->X()*ref->Px()+ref->Y()*ref->Py()>0)? 1.:-1.; //inside or outside
367     if (dir*newdirection<0.5) {
368       nover++;
369       dir = newdirection;
370     }
371     //
372     if (ref->DetectorId()== AliTrackReference::kTRD){
373       tpcmaxRadius =250;
374       fNTRDRef++;
375     }
376     if (ref->DetectorId()== AliTrackReference::kITS){
377       fNITSRef++;
378       tpcminRadius =90;
379     }
380     if (ref->DetectorId()== AliTrackReference::kITS){
381       fNTOFRef++;
382       tpcmaxRadius =250;
383     }
384     //
385     if (ref->DetectorId()== AliTrackReference::kTPC){
386       fNTPCRef++;
387       if (ref->R()>tpcmaxRadius) tpcmaxRadius = ref->R();
388       if (ref->R()<tpcminRadius) tpcminRadius = ref->R();
389     }
390     if (ref->DetectorId()== AliTrackReference::kDisappeared){
391       if (TMath::Abs(ref->Z())<250 && TMath::Abs(ref->R()<250))
392         tpcmaxRadius = ref->R();
393     }
394   }
395   fTPCtrackLength = tpcmaxRadius-tpcminRadius;
396   fMCtracks=nover;
397 }