Parameterization at 8 TeV done by Ara and Vardanush
[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEcalJets.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 /* $Id$ */
17 //
18 // Realisation of AliGenReader to be used with AliGenExtFile
19 // It reads Pythia Jet events from a ntuple like event structure.
20 // The event format is defined in Init()
21 // NextEvent() is used to loop over events and NextParticle() to loop over particles.  
22 // Author: andreas.morsch@cern.ch
23 //
24 #include <TFile.h>
25 #include <TParticle.h>
26 #include <TDatabasePDG.h>
27 #include <TTree.h>
28
29 #include "AliGenReaderEcalJets.h"
30
31 ClassImp(AliGenReaderEcalJets)
32
33
34 AliGenReaderEcalJets::AliGenReaderEcalJets():
35     fNcurrent(0),
36     fNparticle(0),
37     fNev(0),
38     fNpart(0),
39     fNjet(0),     
40     fNsjet(0),
41     fNpjet(0),
42     fTreeNtuple(0) 
43 {
44 // Default constructor
45     for (Int_t i = 0; i < 200; i++) {
46         if (i < 2) {
47             fX[i]     = 0.;    
48             fXtyp[i]  = 0; 
49         } else if (i < 10) {
50             fJet[i]   = 0.;  
51             fJeta[i]  = 0.; 
52             fJphi[i]  = 0.; 
53             fJset[i]  = 0.; 
54             fJseta[i] = 0.;
55             fJsphi[i] = 0.;
56             fJpet[i]  = 0.; 
57             fJpeta[i] = 0.;
58             fJpphi[i] = 0.;
59         } else {
60             fXpt[i]  = 0.;
61             fXeta[i] = 0.;
62             fXphi[i] = 0.;
63             fXid[i]  = 0; 
64         }
65     }
66 }
67
68  AliGenReaderEcalJets::AliGenReaderEcalJets(const AliGenReaderEcalJets &reader)
69      :AliGenReader(reader),
70       fNcurrent(0),
71       fNparticle(0),
72       fNev(0),
73       fNpart(0),
74       fNjet(0),     
75       fNsjet(0),
76       fNpjet(0),
77       fTreeNtuple(0) 
78 {
79     for (Int_t i = 0; i < 200; i++) {
80         if (i < 2) {
81             fX[i]     = 0.;    
82             fXtyp[i]  = 0; 
83         } else if (i < 10) {
84             fJet[i]   = 0.;  
85             fJeta[i]  = 0.; 
86             fJphi[i]  = 0.; 
87             fJset[i]  = 0.; 
88             fJseta[i] = 0.;
89             fJsphi[i] = 0.;
90             fJpet[i]  = 0.; 
91             fJpeta[i] = 0.;
92             fJpphi[i] = 0.;
93         } else {
94             fXpt[i]  = 0.;
95             fXeta[i] = 0.;
96             fXphi[i] = 0.;
97             fXid[i]  = 0; 
98         }
99     }
100
101     // Copy Constructor
102     reader.Copy(*this);
103 }
104
105 void AliGenReaderEcalJets::Init() 
106 {
107 //
108 // reset the existing file environment and open a new root file if
109 // the pointer to the Fluka tree is null
110     
111     TFile *pFile=0;
112     if (!pFile) {
113         pFile = new TFile(fFileName);
114         pFile->cd();
115         printf("\n I have opened %s file \n", fFileName);
116     }
117 // get the tree address in the Fluka boundary source file
118     fTreeNtuple = (TTree*)gDirectory->Get("h1");
119     TTree *h1=fTreeNtuple;
120     h1->SetMakeClass(1);
121 //Set branch addresses
122     h1->SetBranchAddress("nev",   &fNev);
123     h1->SetBranchAddress("x",      fX);
124     h1->SetBranchAddress("xtyp",   fXtyp);
125     h1->SetBranchAddress("npart", &fNpart);
126     h1->SetBranchAddress("xpt",    fXpt);
127     h1->SetBranchAddress("xeta",   fXeta);
128     h1->SetBranchAddress("xphi",   fXphi);
129     h1->SetBranchAddress("xid",    fXid);
130     h1->SetBranchAddress("njet",  &fNjet);
131     h1->SetBranchAddress("jet",    fJet);
132     h1->SetBranchAddress("jeta",   fJeta);
133     h1->SetBranchAddress("jphi",   fJphi);
134     h1->SetBranchAddress("nsjet", &fNsjet);
135     h1->SetBranchAddress("jset",   fJset);
136     h1->SetBranchAddress("jseta",  fJseta);
137     h1->SetBranchAddress("jsphi",  fJsphi);
138     h1->SetBranchAddress("npjet", &fNpjet);
139     h1->SetBranchAddress("jpet",   fJpet);
140     h1->SetBranchAddress("jpeta",  fJpeta);
141     h1->SetBranchAddress("jpphi",  fJpphi);
142 }
143
144 Int_t AliGenReaderEcalJets::NextEvent() 
145 {
146 // Read the next event  
147     Int_t nTracks=0, nread=0;
148     
149     TFile* pFile = fTreeNtuple->GetCurrentFile();
150     pFile->cd();
151
152     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
153     if (fNcurrent < nentries) {
154         Int_t nb = (Int_t)fTreeNtuple->GetEvent(fNcurrent);
155         nread += nb;
156         fNcurrent++;
157         printf("\n Next event contains %d tracks! \n", fNpjet);
158         nTracks    = fNpjet;
159         fNparticle = 0;
160         return nTracks;
161     }
162     return 0;
163 }
164
165 TParticle* AliGenReaderEcalJets::NextParticle() 
166 {
167 // Read the next particle
168
169     Float_t p[4];
170     Int_t    ipart  = fXid[fNparticle];
171     Float_t  pt     = fXpt[fNparticle];
172     Float_t  eta    = fXeta[fNparticle];
173     Float_t  phi    = fXphi[fNparticle];
174     Float_t  theta  = 2.*TMath::ATan(TMath::Exp(-eta));
175     Double_t amass  = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
176
177     p[0] = pt*TMath::Sin(phi);
178     p[1] = pt*TMath::Cos(phi);      
179     p[2] = pt/TMath::Cos(theta);
180     p[3] = TMath::Sqrt(pt*pt+p[2]*p[2]+amass*amass);
181     
182
183     TParticle* particle = 
184         new TParticle(ipart, 0, -1, -1, -1, -1, p[0], p[1], p[2], p[3], 
185                       0., 0., 0., 0.);
186     fNparticle++;
187     return particle;
188 }
189
190
191
192 AliGenReaderEcalJets& AliGenReaderEcalJets::operator=(const  AliGenReaderEcalJets& rhs)
193 {
194 // Assignment operator
195     rhs.Copy(*this);
196     return (*this);
197 }
198
199 void AliGenReaderEcalJets::Copy(TObject&) const
200 {
201     //
202     // Copy 
203     //
204     Fatal("Copy","Not implemented!\n");
205 }
206
207
208
209
210
211
212