Workaround to compute conformal track parameters at vertex.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderKine.cxx
CommitLineData
d7c6ab14 1// $Id$
2
3//_______________________________________________________________________
4/////////////////////////////////////////////////////////////////////////
5//
6// class AliJetParticlesReaderKine
7//
8// Reader for Kinematics
9//
10// loizides@ikf.uni-frankfurt.de
11//
12/////////////////////////////////////////////////////////////////////////
13
14#include <Riostream.h>
bad753b1 15#include <TFile.h>
d7c6ab14 16#include <TString.h>
17#include <TParticle.h>
d7073a83 18#include <TLorentzVector.h>
d7c6ab14 19#include <AliRunLoader.h>
20#include <AliStack.h>
d7073a83 21#include <AliHeader.h>
22#include <AliGenEventHeader.h>
23#include <AliGenPythiaEventHeader.h>
24#include <AliGenHijingEventHeader.h>
d7c6ab14 25#include "AliJetParticle.h"
26#include "AliJetEventParticles.h"
27#include "AliJetParticlesReaderKine.h"
28
29ClassImp(AliJetParticlesReaderKine)
30
d7c6ab14 31AliJetParticlesReaderKine::AliJetParticlesReaderKine() :
32 AliJetParticlesReader(),
33 fFileName("galice.root"),
34 fRunLoader(0),
35 fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
36 fUseTracks(kFALSE)
37{
38 //constructor
39}
40
d7c6ab14 41AliJetParticlesReaderKine::AliJetParticlesReaderKine(TString& fname) :
42 AliJetParticlesReader(),
43 fFileName(fname),
44 fRunLoader(0),
45 fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
46 fUseTracks(kFALSE)
47{
48 //constructor
49}
50
d7c6ab14 51AliJetParticlesReaderKine::AliJetParticlesReaderKine(TObjArray* dirs, const Char_t *filename):
52 AliJetParticlesReader(dirs),
53 fFileName(filename),
54 fRunLoader(0),
55 fNeutral(kTRUE),fCharged(kTRUE),fEM(kTRUE),
56 fUseTracks(kFALSE)
57{
58 //constructor
59}
60
d7c6ab14 61AliJetParticlesReaderKine::~AliJetParticlesReaderKine()
62{
63 //destructor
64 if(fRunLoader) delete fRunLoader;
65}
66
d7c6ab14 67void AliJetParticlesReaderKine::Rewind()
68{
69 //Rewinds to the beginning
70 if(fRunLoader) delete fRunLoader;
04a02430 71 fRunLoader = 0;
72 fCurrentDir = 0;
73 fNEventsRead = 0;
d7c6ab14 74}
75
d7c6ab14 76Int_t AliJetParticlesReaderKine::ReadNext()
77{
78 //Reads Kinematics Tree
79
80 if((!fOwner) || (fEventParticles == 0))
81 fEventParticles = new AliJetEventParticles();
82
b2760c9e 83 while(fCurrentDir < GetNumberOfDirs())
d7c6ab14 84 {
04a02430 85
86 if (!OpenFile(fCurrentDir))
d7c6ab14 87 {
04a02430 88 delete fRunLoader; //close current session
89 fRunLoader = 0; //assure pointer is null
d7c6ab14 90 fCurrentDir++;
91 continue;
92 }
93
94 if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
95 {
96 //read next directory
97 delete fRunLoader; //close current session
98 fRunLoader = 0; //assure pointer is null
99 fCurrentDir++; //go to next dir
100 continue;
101 }
102
a86edc4e 103 Info("ReadNext","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
d7c6ab14 104 fRunLoader->GetEvent(fCurrentEvent);
105 AliStack* stack = fRunLoader->Stack();
106 if (!stack)
107 {
04a02430 108 Error("ReadNext","Can't get stack for event %d",fCurrentEvent);
d7c6ab14 109 continue;
110 }
111
04a02430 112 //clear old event
113 Int_t nprim = stack->GetNprimary();
114 Int_t npart = nprim;
115 if(fUseTracks)
116 npart = stack->GetNtrack();
117 fEventParticles->Reset(npart);
118
119 TString headdesc="";
d7073a83 120 Int_t gentype=-1;
04a02430 121 AliHeader *header=fRunLoader->GetHeader();
122 if(!header) {
123 Warning("ReadNext","Header not found in event %d",fCurrentEvent);
124 } else {
d7073a83 125 AliGenEventHeader *genheader=header->GenEventHeader();
126 TString hname=genheader->GetName();
127 if(hname.CompareTo("HIJINGparam")==0) {
128 gentype=1; ;
129 headdesc+="HIJINGparam";
130 } else if(hname.CompareTo("Hijing")==0) {
131 gentype=2;
132 AliGenHijingEventHeader *hheader=(AliGenHijingEventHeader*)header->GenEventHeader();
133 headdesc+="Hijing";
134 if(!hheader) {
135 Warning("ReadNext","Hijing-Header not found in event %d",fCurrentEvent);
136 } else {
137 Int_t ntrials=hheader->Trials();
138 fEventParticles->SetTrials(ntrials);
139 fEventParticles->SetImpact(hheader->ImpactParameter());
140 fEventParticles->SetNhard(hheader->HardScatters());
141 fEventParticles->SetNpart(hheader->NwNw());
142
143 TLorentzVector jet1;
144 TLorentzVector jet2;
145 TLorentzVector jetFsr1;
146 TLorentzVector jetFsr2;
147 hheader->GetJets(jet1,jet2,jetFsr1,jetFsr2);
148 fEventParticles->AddHard(0,jet1.X(),jet1.Y(),jet1.Z(),jet1.T(),0);
149 fEventParticles->AddHard(1,jet2.X(),jet2.Y(),jet2.Z(),jet2.T(),0);
150 }
151 } else if(hname.CompareTo("Pythia")==0){
152 gentype=3;
153 headdesc+="Pythia";
154 AliGenPythiaEventHeader *pheader=(AliGenPythiaEventHeader*)header->GenEventHeader();
155 if(!pheader) {
156 Warning("ReadNext","Pythia-Header not found in event %d",fCurrentEvent);
157 } else {
158 Int_t ntruq=0;
042d3031 159#ifndef NOUQHEADERINFO
d7073a83 160 ntruq=pheader->NUQTriggerJets();
161 if(ntruq){
162 Double_t x0=pheader->GetXJet();
163 Double_t y0=pheader->GetYJet();
164 if(x0==y0==-1){
165 x0=y0=0.;
166 }
167 Double_t zquench[4];
168 pheader->GetZQuench(zquench);
169 fEventParticles->SetXYJet(x0,y0);
170 fEventParticles->SetZQuench(zquench);
171 for(Int_t j=0;j<ntruq;j++){
172 Float_t pjet[4];
173 pheader->UQJet(j,pjet);
174 fEventParticles->AddUQJet(pjet);
175 }
04a02430 176 }
042d3031 177#endif
d7073a83 178 //Int_t ptyp=pheader->ProcessType();
179 Int_t ntrials=pheader->Trials();
180 fEventParticles->SetTrials(ntrials);
181
182 Int_t ntr=pheader->NTriggerJets();
183 if(ntr){
184 for(Int_t j=0;j<ntr;j++){
185 Float_t pjet[4];
186 pheader->TriggerJet(j,pjet);
187 fEventParticles->AddJet(pjet);
188 if(!ntruq) fEventParticles->AddUQJet(pjet);
189 }
190 }
191 for(Int_t i=6;i<=7;i++){
192 TParticle *MP = stack->Particle(i);
193 if(!MP) break;
194 Int_t type=0;
195 if(MP->GetPdgCode()==21) type=2;
196 else type=1;
197 fEventParticles->AddHard(i-6,MP->Px(),MP->Py(),MP->Pz(),MP->Energy(),type);
04a02430 198 }
199 }
200 }
301a24f1 201 headdesc+=" Run ";
d7073a83 202 headdesc+=header->GetRun();
203 headdesc+=": Ev ";
204 headdesc+=header->GetEventNrInRun();
04a02430 205 }
206 fEventParticles->SetHeader(headdesc);
207
d7c6ab14 208 //get vertex
209 const TParticle *kv = stack->Particle(0);
210 if(kv) {
211 fEventParticles->SetVertex(kv->Vx(),kv->Vy(),kv->Vz());
212 }
213
04a02430 214 //loop over particles
d7c6ab14 215 for (Int_t i = 0;i<npart; i++)
216 {
217 TParticle *p = stack->Particle(i);
218 if(!p) continue;
b2760c9e 219 Int_t child1 = p->GetFirstDaughter();
220 //Int_t child2 = p->GetLastDaughter();
221 //Int_t mother = p->GetFirstMother();
b2760c9e 222 //cout << child1 << " " << child2 << " " << mother << endl;
04a02430 223 if((child1>=0) && (child1<nprim)) continue;
d7073a83 224
225 //check status code depending on gentype
226 if(gentype==1){ // Hijing Param
227 if(p->GetStatusCode()!=0) continue;
228 } else if(gentype==2){ //Hijing
229 //if( p->GetStatusCode()!=0){
230 // cout <<p->GetStatusCode() << " ";p->Print();
231 //}
232 } else if(gentype==3){ //Pythia
233 if(p->GetStatusCode()!=1) continue;
301a24f1 234 p->SetWeight(-123); //mark particle
04a02430 235 }
d7073a83 236
237 //kinematic cuts
238 if(IsAcceptedParticle(p)){ //put particle in event
239 //if(p->Pt()>20){cout <<p->GetStatusCode() << " ";p->Print(); }
75a0c43e 240 fEventParticles->AddParticle(p,i,i);//what is the label in case
241 //events have been merged
d7073a83 242 }
d7c6ab14 243 }
d7c6ab14 244 fCurrentEvent++;
245 fNEventsRead++;
246 return kTRUE;
b2760c9e 247 }
04a02430 248
249 //end of loop over directories specified in fDirs ObjArray
d7c6ab14 250 return kFALSE;
251}
252
d7c6ab14 253Int_t AliJetParticlesReaderKine::OpenFile(Int_t n)
254{
255 //opens file with kine tree
256
04a02430 257 if(fRunLoader){
258 if(fCurrentEvent < fRunLoader->GetNumberOfEvents()) return kTRUE;
259 else return kFALSE;
260 }
261
d7c6ab14 262 const TString& dirname = GetDirName(n);
263 if (dirname == "")
264 {
04a02430 265 Error("OpenNextFile","Can't get directory name with index %d",n);
d7c6ab14 266 return kFALSE;
267 }
268
269 TString filename = dirname +"/"+ fFileName;
bad753b1 270 TFile file(filename);
271 if ( file.IsOpen() == 0)
d7c6ab14 272 {
273 Error("OpenNextFile","Can't open session from file %s",filename.Data());
274 return kFALSE;
275 }
bad753b1 276 file.Close();
277 fRunLoader = AliRunLoader::Open(filename.Data());
d7c6ab14 278
279 if (fRunLoader->GetNumberOfEvents() <= 0)
280 {
281 Error("OpenNextFile","There is no events in this directory.");
282 delete fRunLoader;
283 fRunLoader = 0;
284 return kFALSE;
285 }
04a02430 286
d7c6ab14 287 if (fRunLoader->LoadKinematics())
288 {
289 Error("OpenNextFile","Error occured while loading kinematics.");
290 return kFALSE;
291 }
04a02430 292
d7c6ab14 293 fCurrentEvent = 0;
294 return kTRUE;
295}