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