EffC++ warnings corrected. (M. Lopez Noriega)
[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.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 // Jet kinematics reader 
18 // MC reader for jet analysis
19 // Author: Andreas Morsch (andreas.morsch@cern.ch)
20 //-------------------------------------------------------------------------
21
22 // From root ...
23 #include <TClonesArray.h>
24 #include <TPDGCode.h>
25 #include <TParticle.h>
26 #include <TParticlePDG.h>
27 //#include <TVector3.h>
28 #include <TLorentzVector.h>
29 #include <TSystem.h>
30 #include <TDatabasePDG.h>
31 #include <TRandom.h>
32 // From AliRoot ...
33 #include "AliJetKineReaderHeader.h"
34 #include "AliJetKineReader.h"
35 #include "AliRunLoader.h"
36 #include "AliStack.h"
37
38 ClassImp(AliJetKineReader)
39
40 AliJetKineReader::AliJetKineReader():
41   fRunLoader(0x0),
42   fMass(0),
43   fPdgC(0)
44 {
45   // Default constructor
46   fReaderHeader = 0x0;
47 }
48
49 //____________________________________________________________________________
50
51 AliJetKineReader::~AliJetKineReader()
52 {
53   // Destructor
54   delete fReaderHeader;
55 }
56
57 //____________________________________________________________________________
58
59 void AliJetKineReader::OpenInputFiles()
60 {
61   // Opens the input file using the run loader
62   const char* dirName = fReaderHeader->GetDirectory();
63   char path[256];
64   sprintf(path, "%s/galice.root",dirName);
65   fRunLoader = AliRunLoader::Open(path);
66   fRunLoader->LoadKinematics();
67   fRunLoader->LoadHeader(); 
68   
69   // get number of events
70   Int_t nMax = fRunLoader->GetNumberOfEvents();
71   printf("\nTotal number of events = %d", nMax);
72   
73   // set number of events in header
74   if (fReaderHeader->GetLastEvent() == -1)
75     fReaderHeader->SetLastEvent(nMax);
76   else {
77     Int_t nUsr = fReaderHeader->GetLastEvent();
78     fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
79   }
80 }
81
82 //____________________________________________________________________________
83
84 void AliJetKineReader::FillMomentumArray(Int_t event)
85 {
86   // Fill momentum array for event
87  
88   Int_t goodTrack = 0;
89   // Clear array
90   ClearArray();
91   // Get event from runloader
92   fRunLoader->GetEvent(event);
93   // Get the stack
94   AliStack* stack = fRunLoader->Stack();
95   // Number of primaries
96   Int_t nt = stack->GetNprimary();
97
98   // Get cuts set by user and header
99   Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
100   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
101   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
102   fAliHeader = fRunLoader->GetHeader();
103   
104   // temporary storage of signal and cut flags
105   Int_t* sflag  = new Int_t[nt];
106   Int_t* cflag  = new Int_t[nt];
107   
108   TLorentzVector p4;
109   // Loop over particles
110   for (Int_t it = 0; it < nt; it++) {
111     TParticle *part = stack->Particle(it);
112     Int_t   status  = part->GetStatusCode();
113     Int_t   pdg     = TMath::Abs(part->GetPdgCode());
114     Float_t pt      = part->Pt(); 
115     
116     // Skip non-final state particles, neutrinos 
117     // and particles with pt < pt_min 
118     if ((status != 1)            
119         || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
120     
121     Float_t p       = part->P();
122     Float_t p0      = p;
123     Float_t eta     = part->Eta();
124     Float_t phi     = part->Phi();
125     
126     // Fast simulation of TPC if requested
127     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
128       // Charged particles only
129       Float_t charge = 
130         TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
131       if (charge == 0)               continue;
132       // Simulate efficiency
133       if (!Efficiency(p0, eta, phi)) continue;
134       // Simulate resolution
135       p = SmearMomentum(4, p0);
136     } // End fast TPC
137     
138     // Fast simulation of EMCAL if requested
139     if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
140       Float_t charge = 
141         TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
142       // Charged particles only
143       if (charge != 0){
144         // Simulate efficiency
145         if (!Efficiency(p0, eta, phi)) continue;
146         // Simulate resolution
147         p = SmearMomentum(4, p0);
148       } // end "if" charged particles
149       // Neutral particles (exclude K0L, n, nbar)
150       if (pdg == kNeutron || pdg == kK0Long) continue;
151     } // End fast EMCAL
152     
153     fMass = part->GetCalcMass();
154     fPdgC = part->GetPdgCode();
155
156     // Fill momentum array
157     Float_t r  = p/p0;
158     Float_t px = r * part->Px(); 
159     Float_t py = r * part->Py(); 
160     Float_t pz = r * part->Pz();
161     Float_t m  = part->GetMass();
162     Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
163     p4 = TLorentzVector(px, py, pz, e);
164     if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
165     new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
166     
167     // all tracks are signal ... ?
168     sflag[goodTrack]=1;
169     cflag[goodTrack]=0;
170     if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
171     goodTrack++;
172   }
173
174   // set the signal flags
175   fSignalFlag.Set(goodTrack,sflag);
176   fCutFlag.Set(goodTrack,cflag);
177   printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
178
179 }
180
181
182 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
183 {
184   //  The relative momentum error, i.e. 
185   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
186   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
187   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
188   //
189   //  If we include information from TRD, b will be a factor 2/3 smaller.
190   //
191   //  ind = 1: high multiplicity
192   //  ind = 2: low  multiplicity
193   //  ind = 3: high multiplicity + TRD
194   //  ind = 4: low  multiplicity + TRD
195   
196   Float_t pSmeared;
197   Float_t a = 0.75;
198   Float_t b = 0.12;
199   
200   if (ind == 1) b = 0.12;
201   if (ind == 2) b = 0.08;
202   if (ind == 3) b = 0.12;    
203   if (ind == 4) b = 0.08;    
204   
205   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
206   pSmeared = p + gRandom->Gaus(0., sigma);
207   return pSmeared;
208 }
209
210
211 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
212 {
213   // Fast simulation of geometrical acceptance and tracking efficiency
214  
215   //  Tracking
216   Float_t eff = 0.99;
217   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
218   // Geometry
219   if (p > 0.5) {
220     phi *= 180. / TMath::Pi();
221     // Sector number 0 - 17
222     Int_t isec  = Int_t(phi / 20.);
223     // Sector centre
224     Float_t phi0 = isec * 20. + 10.;
225     Float_t phir = TMath::Abs(phi-phi0);
226     // 2 deg of dead space
227     if (phir > 9.) eff = 0.;
228   }
229
230   if (gRandom->Rndm() > eff) {
231     return kFALSE;
232   } else {
233     return kTRUE;
234   }
235
236 }
237