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