Corrections in the destructors (J. Thaeder)
[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 <TLorentzVector.h>
28 #include <TSystem.h>
29 #include <TDatabasePDG.h>
30 #include <TRandom.h>
31 // From AliRoot ...
32 #include "AliAODJet.h"
33 #include "AliPDG.h"
34 #include "AliJetKineReaderHeader.h"
35 #include "AliJetKineReader.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliStack.h"
39 #include "AliHeader.h"
40 #include "AliGenPythiaEventHeader.h"
41 #include "AliLog.h"
42
43 ClassImp(AliJetKineReader)
44
45 AliJetKineReader::AliJetKineReader():
46     AliJetReader(),  
47     fAliHeader(0),
48     fMCEvent(0),
49     fGenJets(0)
50 {
51   // Default constructor
52 }
53
54 //____________________________________________________________________________
55
56 AliJetKineReader::~AliJetKineReader()
57 {
58   // Destructor
59     if (fAliHeader) {
60         delete fAliHeader;
61         fAliHeader = 0;
62     }
63 }
64
65 //____________________________________________________________________________
66
67 void AliJetKineReader::OpenInputFiles()
68 {
69   // Opens the input file using the run loader
70   // Not used anymore   
71 }
72
73 //____________________________________________________________________________
74
75 Bool_t AliJetKineReader::FillMomentumArray()
76 {
77   // Fill momentum array for event
78     Int_t goodTrack = 0;
79     // Clear array
80     // temporary storage of signal and cut flags
81     Int_t* sflag  = new Int_t[50000];
82     Int_t* cflag  = new Int_t[50000];
83   
84     ClearArray();
85     // Get the stack
86     AliStack* stack = fMCEvent->Stack();
87     // Number of primaries
88     Int_t nt = stack->GetNprimary();
89       
90     // Get cuts set by user and header
91     Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
92     Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
93     Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
94     fAliHeader = fMCEvent->Header();
95       
96     
97     TLorentzVector p4;
98     // Loop over particles
99     for (Int_t it = 0; it < nt; it++) {
100         TParticle *part = stack->Particle(it);
101         Int_t   status  = part->GetStatusCode();
102         Int_t   pdg     = TMath::Abs(part->GetPdgCode());
103         Float_t pt      = part->Pt(); 
104           
105         // Skip non-final state particles, neutrinos 
106         // and particles with pt < pt_min 
107         if ((status != 1)            
108             || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
109         
110         Float_t p       = part->P();
111         Float_t p0      = p;
112         Float_t eta     = part->Eta();
113         Float_t phi     = part->Phi();
114         
115         // Fast simulation of TPC if requested
116         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
117             // Charged particles only
118             Float_t charge = 
119                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
120             if (charge == 0)               continue;
121             // Simulate efficiency
122             if (!Efficiency(p0, eta, phi)) continue;
123             // Simulate resolution
124             p = SmearMomentum(4, p0);
125         } // End fast TPC
126           
127           // Fast simulation of EMCAL if requested
128         if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
129             Float_t charge = 
130                 TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
131             // Charged particles only
132             if (charge != 0){
133                 // Simulate efficiency
134                 if (!Efficiency(p0, eta, phi)) continue;
135                 // Simulate resolution
136                 p = SmearMomentum(4, p0);
137             } // end "if" charged particles
138               // Neutral particles (exclude K0L, n, nbar)
139             if (pdg == kNeutron || pdg == kK0Long) continue;
140         } // End fast EMCAL
141         
142           // Fill momentum array
143         Float_t r  = p/p0;
144         Float_t px = r * part->Px(); 
145         Float_t py = r * part->Py(); 
146         Float_t pz = r * part->Pz();
147         TParticlePDG *pPDG = part->GetPDG();
148         Float_t m = 0;
149         if(!pPDG){
150           // this is very rare...
151           // Is avoided by AliPDG::AddParticlesToPdgDataBase();
152           // but this should be called only once... (how in proof?)
153           // Calucluate mass with unsmeared momentum values
154           m  = part->Energy()*part->Energy() - 
155             (px * px + py * py + pz * pz)/r/r; 
156           if(m>0)m = TMath::Sqrt(m);
157           else m = 0;
158           AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));
159
160         }
161         else{
162           m  = pPDG->Mass();
163         }
164         Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
165         p4 = TLorentzVector(px, py, pz, e);
166         if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
167         new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
168         
169         // all tracks are signal ... ?
170         sflag[goodTrack]=1;
171         cflag[goodTrack]=0;
172         if (pt > ptMin) cflag[goodTrack] = 1; // track surviving pt cut
173         goodTrack++;
174     } // track loop
175
176 // set the signal flags
177     fSignalFlag.Set(goodTrack, sflag);
178     fCutFlag.Set(goodTrack, cflag);
179     AliInfo(Form(" Number of good tracks %d \n", goodTrack));
180     
181     delete[] sflag;
182     delete[] cflag;
183     
184     return kTRUE;
185 }
186
187
188 Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
189 {
190   //  The relative momentum error, i.e. 
191   //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
192   //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
193   //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
194   //
195   //  If we include information from TRD, b will be a factor 2/3 smaller.
196   //
197   //  ind = 1: high multiplicity
198   //  ind = 2: low  multiplicity
199   //  ind = 3: high multiplicity + TRD
200   //  ind = 4: low  multiplicity + TRD
201   
202   Float_t pSmeared;
203   Float_t a = 0.75;
204   Float_t b = 0.12;
205   
206   if (ind == 1) b = 0.12;
207   if (ind == 2) b = 0.08;
208   if (ind == 3) b = 0.12;    
209   if (ind == 4) b = 0.08;    
210   
211   Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
212   pSmeared = p + gRandom->Gaus(0., sigma);
213   return pSmeared;
214 }
215
216
217 Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
218 {
219   // Fast simulation of geometrical acceptance and tracking efficiency
220  
221   //  Tracking
222   Float_t eff = 0.99;
223   if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
224   // Geometry
225   if (p > 0.5) {
226     phi *= 180. / TMath::Pi();
227     // Sector number 0 - 17
228     Int_t isec  = Int_t(phi / 20.);
229     // Sector centre
230     Float_t phi0 = isec * 20. + 10.;
231     Float_t phir = TMath::Abs(phi-phi0);
232     // 2 deg of dead space
233     if (phir > 9.) eff = 0.;
234   }
235
236   if (gRandom->Rndm() > eff) {
237     return kFALSE;
238   } else {
239     return kTRUE;
240   }
241
242 }
243
244 TClonesArray*  AliJetKineReader::GetGeneratedJets()
245 {
246     //
247     // Get generated jets from mc header
248     //
249     if (!fGenJets) {
250         fGenJets = new TClonesArray("AliAODJets", 100);
251     } else {
252         fGenJets->Clear();
253     }
254     
255     
256     AliHeader* alih = GetAliHeader(); 
257     if (alih == 0) return 0;
258     AliGenEventHeader * genh = alih->GenEventHeader();
259     if (genh == 0) return 0;
260
261     Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
262     for (Int_t i = 0; i < nj; i++) {
263         Float_t p[4];
264         ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
265         new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
266     }
267
268     return fGenJets;
269 }
270
271 void AliJetKineReader::SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* mc)
272 {
273     // Connect the input event
274     fMCEvent = (AliMCEvent*) mc;
275 }