First attempt to use AliAnalisysTask (Marian)
[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()
141 {
142   //
143   // Destructor of the class
144   //
145   if (fTPCReferences) {
146     delete fTPCReferences;
147   }
148   if (fITSReferences){
149     delete fITSReferences;
150   }
151   if (fTRDReferences){    
152     delete fTRDReferences;  
153   }
154   if (fTOFReferences){    
155     delete fTOFReferences;  
156   }
157
158 }
159
160
161
162 void AliMCInfo::Update()
163 {
164   //
165   // Update MC info
166   // Calculates some derived variables
167   //
168   //
169   fMCtracks =1;
170   //
171   // decay info
172   fTPCdecay=kFALSE;
173   if (fTRdecay.GetTrack()>0){
174     fDecayCoord[0] = fTRdecay.X();
175     fDecayCoord[1] = fTRdecay.Y();
176     fDecayCoord[2] = fTRdecay.Z();
177     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
178       fTPCdecay=kTRUE;     
179     }
180     else{
181       fDecayCoord[0] = 0;
182       fDecayCoord[1] = 0;
183       fDecayCoord[2] = 0;
184     }
185   }
186   //
187   //
188   //digits information update
189   fRowsWithDigits    = fTPCRow.RowsOn();    
190   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
191   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
192   //
193   //
194   // calculate primary ionization per cm  
195   if (fParticle.GetPDG()){
196     fMass = fParticle.GetMass();  
197     fCharge = fParticle.GetPDG()->Charge();
198     if (fTPCReferences->GetEntriesFast()>0){
199       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
200     }
201     if (fMass>0){
202       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
203                               fTrackRef.Py()*fTrackRef.Py()+
204                               fTrackRef.Pz()*fTrackRef.Pz());    
205       if (p>0.001){
206         Float_t betagama = p /fMass;
207         fPrim = TPCBetheBloch(betagama);
208       }else fPrim=0;
209     }
210   }else{
211     fMass =0;
212     fPrim =0;
213   }  
214 }
215
216
217
218
219 ////////////////////////////////////////////////////////////////////////
220 AliTPCdigitRow::AliTPCdigitRow()
221 {
222   Reset();
223 }
224 ////////////////////////////////////////////////////////////////////////
225 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
226 {
227   for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
228   return (*this);
229 }
230 ////////////////////////////////////////////////////////////////////////
231 void AliTPCdigitRow::SetRow(Int_t row) 
232 {
233   //
234   // set bit mask for given row
235   //
236   if (row >= 8*32) {
237     //    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
238     return;
239   }
240   Int_t iC = row/8;
241   Int_t iB = row%8;
242   SETBIT(fDig[iC],iB);
243 }
244
245 ////////////////////////////////////////////////////////////////////////
246 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
247 {
248 //
249 // return kTRUE if row is on
250 //
251   Int_t iC = row/8;
252   Int_t iB = row%8;
253   return TESTBIT(fDig[iC],iB);
254 }
255 ////////////////////////////////////////////////////////////////////////
256 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
257 {
258 //
259 // returns number of rows with a digit  
260 // count only rows less equal row number upto
261 //
262   Int_t total = 0;
263   for (Int_t i = 0; i<32; i++) {
264     for (Int_t j = 0; j < 8; j++) {
265       if (i*8+j > upto) return total;
266       if (TESTBIT(fDig[i],j))  total++;
267     }
268   }
269   return total;
270 }
271 ////////////////////////////////////////////////////////////////////////
272 void AliTPCdigitRow::Reset()
273 {
274 //
275 // resets all rows to zero
276 //
277   for (Int_t i = 0; i<32; i++) {
278     fDig[i] <<= 8;
279   }
280 }
281 ////////////////////////////////////////////////////////////////////////
282 Int_t AliTPCdigitRow::Last() const
283 {
284 //
285 // returns the last row number with a digit
286 // returns -1 if now digits 
287 //
288   for (Int_t i = 32-1; i>=0; i--) {
289     for (Int_t j = 7; j >= 0; j--) {
290       if TESTBIT(fDig[i],j) return i*8+j;
291     }
292   }
293   return -1;
294 }
295 ////////////////////////////////////////////////////////////////////////
296 Int_t AliTPCdigitRow::First() const
297 {
298 //
299 // returns the first row number with a digit
300 // returns -1 if now digits 
301 //
302   for (Int_t i = 0; i<32; i++) {
303     for (Int_t j = 0; j < 8; j++) {
304       if (TESTBIT(fDig[i],j)) return i*8+j;
305     }
306   }
307   return -1;
308 }
309
310
311 //_____________________________________________________________________________
312 Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
313 {
314   //
315   // Bethe-Bloch energy loss formula
316   //
317   const Double_t kp1=0.76176e-1;
318   const Double_t kp2=10.632;
319   const Double_t kp3=0.13279e-4;
320   const Double_t kp4=1.8631;
321   const Double_t kp5=1.9479;
322
323   Double_t dbg = (Double_t) bg;
324
325   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
326
327   Double_t aa = TMath::Power(beta,kp4);
328   Double_t bb = TMath::Power(1./dbg,kp5);
329
330   bb=TMath::Log(kp3+bb);
331   
332   return ((Float_t)((kp2-aa-bb)*kp1/aa));
333 }
334
335
336 void AliMCInfo::CalcTPCrows(TClonesArray * runArrayTR){
337   //
338   // Calculates the numebr of the track references for detectors
339   // In case of changing direction - curling tracks - the counter is not increasing
340   //
341   // Rough calculation 
342   // of the first and last point in the TPC  
343   //
344   fNTPCRef = 0;
345   fNITSRef = 0;
346   fNTRDRef = 0;
347   fNTOFRef = 0;  
348   Float_t tpcminRadius=250;
349   Float_t tpcmaxRadius=80;
350   Float_t dir=0;
351   Int_t nover=0;
352   
353   for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
354     //
355     AliTrackReference *ref = (AliTrackReference*)runArrayTR->At(iTrackRef);  
356     Float_t newdirection = (ref->X()*ref->Px()+ref->Y()*ref->Py()>0)? 1.:-1.; //inside or outside
357     if (dir*newdirection<0.5) {
358       nover++;
359       dir = newdirection;
360     }
361     //
362     if (ref->DetectorId()== AliTrackReference::kTRD){
363       tpcmaxRadius =250;
364       fNTRDRef++;
365     }
366     if (ref->DetectorId()== AliTrackReference::kITS){
367       fNITSRef++;
368       tpcminRadius =90;
369     }
370     if (ref->DetectorId()== AliTrackReference::kITS){
371       fNTOFRef++;
372       tpcmaxRadius =250;
373     }
374     //
375     if (ref->DetectorId()== AliTrackReference::kTPC){
376       fNTPCRef++;
377       if (ref->R()>tpcmaxRadius) tpcmaxRadius = ref->R();
378       if (ref->R()<tpcminRadius) tpcminRadius = ref->R();
379     }
380     if (ref->DetectorId()== AliTrackReference::kDisappeared){
381       if (TMath::Abs(ref->Z())<250 && TMath::Abs(ref->R()<250))
382         tpcmaxRadius = ref->R();
383     }
384   }
385   fTPCtrackLength = tpcmaxRadius-tpcminRadius;
386   fMCtracks=nover;
387 }