]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/AliHBTReaderPPprod.cxx
Updated Linkdef and libTOF.pkg
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderPPprod.cxx
CommitLineData
1b446896 1#include "AliHBTReaderPPprod.h"
2
3#include <iostream.h>
4#include <fstream.h>
5#include <TString.h>
6#include <TTree.h>
7#include <TFile.h>
8
9#include <AliTPCtrack.h>
10#include <AliTPCParam.h>
11#include <AliTPCtracker.h>
12
13#include "AliRun.h"
14#include "AliHBTRun.h"
15#include "AliHBTEvent.h"
16#include "AliHBTParticle.h"
17#include "AliHBTParticleCut.h"
18
19
20ClassImp(AliHBTReaderPPprod)
21
22AliHBTReaderPPprod::
23 AliHBTReaderPPprod(const Char_t* trackfilename,const Char_t* clusterfilename,
24 const Char_t* goodtracksfilename,const Char_t* galicefilename):
25 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
26 fGAliceFileName(galicefilename),
27 fGoodTPCTracksFileName(goodtracksfilename)
28{
29 //constructor, only file names are set
30 //Defaults:
31 // trackfilename = "AliTPCtracks.root"
32 // clusterfilename = "AliTPCclusters.root"
33 // goodtracksfilename = "good_tracks_tpc"
34 // galicefilename = "" - this means: Do not open gAlice file -
35 // just leave the global pointer untached
36
37 fParticles = new AliHBTRun();
38 fTracks = new AliHBTRun();
39
40 fTracksFile = 0x0; //files are opened during reading only
41 fClustersFile = 0x0;
42
43 fIsRead = kFALSE;
44}
45/********************************************************************/
46
47AliHBTReaderPPprod::~AliHBTReaderPPprod()
48 {
49 delete fParticles;
50 delete fTracks;
51 }
52/********************************************************************/
53
54AliHBTEvent* AliHBTReaderPPprod::GetParticleEvent(Int_t n)
55 {
56 //returns Nth event with simulated particles
57 if (!fIsRead) Read(fParticles,fTracks);
58 return fParticles->GetEvent(n);
59 }
60/********************************************************************/
61AliHBTEvent* AliHBTReaderPPprod::GetTrackEvent(Int_t n)
62 {
63 //returns Nth event with reconstructed tracks
64 if (!fIsRead) Read(fParticles,fTracks);
65 return fTracks->GetEvent(n);
66 }
67/********************************************************************/
68
69Int_t AliHBTReaderPPprod::GetNumberOfPartEvents()
70 {
71 //returns number of events of particles
72 if (!fIsRead) Read(fParticles,fTracks);
73 return fParticles->GetNumberOfEvents();
74 }
75
76/********************************************************************/
77Int_t AliHBTReaderPPprod::GetNumberOfTrackEvents()
78 {
79 //returns number of events of tracks
80 if (!fIsRead) Read(fParticles,fTracks);
81 return fTracks->GetNumberOfEvents();
82 }
83/********************************************************************/
84
85
86Int_t AliHBTReaderPPprod::Read(AliHBTRun* particles, AliHBTRun *tracks)
87 {
88 //reads data and puts put to the particles and tracks objects
89 //reurns 0 if everything is OK
90 //
91 Int_t i; //iterator and some temprary values
92 Int_t Nevents;
93 if (!particles) //check if an object is instatiated
94 {
95 Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
96 }
97 if (!tracks) //check if an object is instatiated
98 {
99 Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
100 }
101 particles->Reset();//clear runs == delete all old events
102 tracks->Reset();
103
104 if( (i=OpenFiles()) )
105 {
106 Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
107 return i;
108 }
109
110 AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
111 if (!goodTPCTracks)
112 {
113 Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
114 return 1;
115 }
116
117 Nevents = 100;
118
119 fClustersFile->cd();
120 AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
121 if (!TPCParam)
122 {
123 Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
124 return 1;
125 }
126
127 TObjArray *tarray = new TObjArray(5000);
128 tarray->SetOwner();// this causes memory leak, but in some cases deleting is infinite loop
129
130
131 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
132 {
133 cout<<"Reading Event "<<currentEvent<<endl;
134 /**************************************/
135 /**************************************/
136 /**************************************/
137 fTracksFile->cd();
138 Char_t treename[100];
139 sprintf(treename,"TreeT_TPC_%d",currentEvent);
140
141 TTree *tracktree=0;
142
143 tracktree=(TTree*)fTracksFile->Get(treename);
144 if (!tracktree)
145 {
146 Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n");
147 return 1;
148 }
149
150 TBranch *trackbranch=tracktree->GetBranch("tracks");
151 if (!trackbranch)
152 {
153 Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n");
154 return 2;
155 }
156 Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
157 cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
158 //Copy tracks to array
159
160 AliTPCtrack *iotrack=0;
161
162 fClustersFile->cd();
163 AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);
164 if (!tracker)
165 {
166 Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n");
167 return 3;
168 }
169 tracker->LoadInnerSectors();
170 tracker->LoadOuterSectors();
171
172 for (i=0; i<NTPCtracks; i++)
173 {
174 iotrack=new AliTPCtrack;
175 trackbranch->SetAddress(&iotrack);
176 tracktree->GetEvent(i);
177 tracker->CookLabel(iotrack,0.1);
178 tarray->AddLast(iotrack);
179 }
180
181
182 fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
183 fTracksFile->Delete("tracks");
184
185 delete tracker;
186
187 tracker = 0x0;
188 trackbranch = 0x0;
189 tracktree = 0x0;
190
191 Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
192
193 Double_t xk;
194 Double_t par[5];
195 Float_t phi, lam, pt;
196 Int_t label;
197 Bool_t found;
198
199 for (i=0; i<ngood; i++)
200 {
201 const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
202
203 if(Pass(gt.code)) continue;
204
205 label = gt.lab;
206 found = kFALSE; //guard in case we don't find track with such a label
207 for (Int_t j=0;j<NTPCtracks;j++)
208 {
209 iotrack = (AliTPCtrack*)tarray->At(j);
210 if (iotrack->GetLabel() == label)
211 {
212 found = kTRUE;
213 break;
214 }
215 }
216 if(!found)
217 {
218 Warning("Read",
219 "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
220 continue; //put comunicate on the screen and continue loop
221 }
222
223 Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
224 Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
225
226 AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
227 if(Pass(part)) continue;
228
229
230 iotrack->PropagateTo(gt.x);
231 iotrack->GetExternalParameters(xk,par);
232 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
233 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
234 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
235 lam=par[3];
236 pt=1.0/TMath::Abs(par[4]);
237
238 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
239 Double_t tpy = pt * TMath::Sin(phi); //track x coordinate of momentum
240 Double_t tpz = pt * lam;
241
242 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
243
244 AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
245 if(Pass(track)) continue;
246
247 particles->AddParticle(currentEvent,part);
248 tracks->AddParticle(currentEvent,track);
249
250 }
251 tarray->Clear();
252 /**************************************/
253 /**************************************/
254 /**************************************/
255 }
256
257 CloseFiles();
258 delete tarray;
259 delete goodTPCTracks;
260 fIsRead = kTRUE;
261 return 0;
262 }
263
264/********************************************************************/
265Int_t AliHBTReaderPPprod::OpenFiles()
266{
267
268 fTracksFile = 0;
269 fTracksFile=TFile::Open("AliTPCtracks.root");
270 if (!fTracksFile->IsOpen())
271 {
272 Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
273 return 1;
274 }
275
276 fClustersFile = 0;
277
278 fClustersFile=TFile::Open("AliTPCclusters.root");
279 if (!fClustersFile->IsOpen())
280 {
281 Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
282 return 2;
283 }
284
285
286
287 return 0;
288}
289
290
291
292/********************************************************************/
293
294void AliHBTReaderPPprod::CloseFiles()
295{
296 fTracksFile->Close();
297 fClustersFile->Close();
298}
299
300/********************************************************************/
301
302/********************************************************************/
303/********************************************************************/
304/********************************************************************/
305
306
307AliGoodTracksPP::~AliGoodTracksPP()
308{
309 delete [] fGoodInEvent;
310 for (Int_t i = 0;i<fNevents;i++)
311 delete [] fData[i];
312 delete [] fData;
313}
314/********************************************************************/
315AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
316{
317
318 fNevents = 100;
319 ifstream in(infilename.Data());
320
321 if(!in)
322 {
323 cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
324 delete this;
325 return;
326 }
327
328
329 fGoodInEvent = new Int_t[fNevents];
330 fData = new struct GoodTrack* [fNevents];
331
332 Int_t i;
333 for( i = 0;i<fNevents;i++)
334 {
335 fGoodInEvent[i] =0;
336 fData[i] = new struct GoodTrack[500];
337 }
338 Float_t tmp;
339 Int_t evno;
340 while(in>>evno)
341 {
342 if(fGoodInEvent[evno]>=500)
343 {
344 cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
345 continue;
346 }
347 in>>fData[evno][fGoodInEvent[evno]].lab;
348 in>>fData[evno][fGoodInEvent[evno]].code;
349 in>>fData[evno][fGoodInEvent[evno]].px;
350 in>>fData[evno][fGoodInEvent[evno]].py;
351 in>>fData[evno][fGoodInEvent[evno]].pz;
352 in>>fData[evno][fGoodInEvent[evno]].x;
353 in>>fData[evno][fGoodInEvent[evno]].y;
354 in>>fData[evno][fGoodInEvent[evno]].z;
355
356 in>>tmp;
357 in>>tmp;
358 in>>tmp;
359
360 /* cout<<evno<<" ";
361 cout<<fData[evno][fGoodInEvent[evno]].lab;
362 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
363 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
364 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
365 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
366 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
367 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
368 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
369 cout<<"\n";
370 */
371 fGoodInEvent[evno]++;
372 }
373 in.close();
374 cout<<"AliGoodTracksPP::AliGoodTracksPP() .... Done\n";
375}
376
377
378
379const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
380 {
381
382 if( (event>fNevents) || (event<0))
383 {
384 gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
385 }
386 if( (n>fGoodInEvent[event]) || (n<0))
387 {
388 gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
389 }
390 return fData[event][n];
391
392 }