Transition to NewIO
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderPPprod.cxx
CommitLineData
1b446896 1#include "AliHBTReaderPPprod.h"
2
d0c23b58 3#include <Riostream.h>
4#include <Riostream.h>
1b446896 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 //
88cb7938 91 cout<<"New I/O not implemented\n";
92 return 0;
93
1b446896 94 Int_t i; //iterator and some temprary values
95 Int_t Nevents;
96 if (!particles) //check if an object is instatiated
97 {
98 Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
99 }
100 if (!tracks) //check if an object is instatiated
101 {
102 Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
103 }
104 particles->Reset();//clear runs == delete all old events
105 tracks->Reset();
106
107 if( (i=OpenFiles()) )
108 {
109 Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
110 return i;
111 }
112
113 AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
114 if (!goodTPCTracks)
115 {
116 Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
117 return 1;
118 }
119
120 Nevents = 100;
121
122 fClustersFile->cd();
123 AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
124 if (!TPCParam)
125 {
126 Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
127 return 1;
128 }
129
130 TObjArray *tarray = new TObjArray(5000);
131 tarray->SetOwner();// this causes memory leak, but in some cases deleting is infinite loop
132
133
134 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
135 {
136 cout<<"Reading Event "<<currentEvent<<endl;
137 /**************************************/
138 /**************************************/
139 /**************************************/
140 fTracksFile->cd();
141 Char_t treename[100];
142 sprintf(treename,"TreeT_TPC_%d",currentEvent);
143
144 TTree *tracktree=0;
145
146 tracktree=(TTree*)fTracksFile->Get(treename);
147 if (!tracktree)
148 {
149 Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n");
150 return 1;
151 }
152
153 TBranch *trackbranch=tracktree->GetBranch("tracks");
154 if (!trackbranch)
155 {
156 Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n");
157 return 2;
158 }
159 Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
160 cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
161 //Copy tracks to array
162
163 AliTPCtrack *iotrack=0;
164
165 fClustersFile->cd();
88cb7938 166 AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,"");
1b446896 167 if (!tracker)
168 {
169 Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n");
170 return 3;
171 }
88cb7938 172 tracker->LoadInnerSectors();
173 tracker->LoadOuterSectors();
1b446896 174
175 for (i=0; i<NTPCtracks; i++)
176 {
177 iotrack=new AliTPCtrack;
178 trackbranch->SetAddress(&iotrack);
179 tracktree->GetEvent(i);
180 tracker->CookLabel(iotrack,0.1);
181 tarray->AddLast(iotrack);
182 }
183
184
185 fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
186 fTracksFile->Delete("tracks");
187
188 delete tracker;
189
190 tracker = 0x0;
191 trackbranch = 0x0;
192 tracktree = 0x0;
193
194 Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
195
196 Double_t xk;
197 Double_t par[5];
198 Float_t phi, lam, pt;
199 Int_t label;
200 Bool_t found;
201
202 for (i=0; i<ngood; i++)
203 {
204 const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
205
206 if(Pass(gt.code)) continue;
207
208 label = gt.lab;
209 found = kFALSE; //guard in case we don't find track with such a label
210 for (Int_t j=0;j<NTPCtracks;j++)
211 {
212 iotrack = (AliTPCtrack*)tarray->At(j);
213 if (iotrack->GetLabel() == label)
214 {
215 found = kTRUE;
216 break;
217 }
218 }
219 if(!found)
220 {
221 Warning("Read",
222 "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
223 continue; //put comunicate on the screen and continue loop
224 }
225
226 Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
227 Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
228
88cb7938 229 AliHBTParticle* part = new AliHBTParticle(gt.code,i, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
1b446896 230 if(Pass(part)) continue;
231
232
233 iotrack->PropagateTo(gt.x);
234 iotrack->GetExternalParameters(xk,par);
235 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
236 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
237 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
238 lam=par[3];
239 pt=1.0/TMath::Abs(par[4]);
240
241 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
242 Double_t tpy = pt * TMath::Sin(phi); //track x coordinate of momentum
243 Double_t tpz = pt * lam;
244
245 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
246
88cb7938 247 AliHBTParticle* track = new AliHBTParticle(gt.code,i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
1b446896 248 if(Pass(track)) continue;
249
250 particles->AddParticle(currentEvent,part);
251 tracks->AddParticle(currentEvent,track);
252
253 }
254 tarray->Clear();
255 /**************************************/
256 /**************************************/
257 /**************************************/
258 }
259
260 CloseFiles();
261 delete tarray;
262 delete goodTPCTracks;
263 fIsRead = kTRUE;
264 return 0;
265 }
266
267/********************************************************************/
268Int_t AliHBTReaderPPprod::OpenFiles()
269{
270
271 fTracksFile = 0;
272 fTracksFile=TFile::Open("AliTPCtracks.root");
273 if (!fTracksFile->IsOpen())
274 {
275 Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
276 return 1;
277 }
278
279 fClustersFile = 0;
280
281 fClustersFile=TFile::Open("AliTPCclusters.root");
282 if (!fClustersFile->IsOpen())
283 {
284 Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
285 return 2;
286 }
287
288
289
290 return 0;
291}
292
293
294
295/********************************************************************/
296
297void AliHBTReaderPPprod::CloseFiles()
298{
299 fTracksFile->Close();
300 fClustersFile->Close();
301}
302
303/********************************************************************/
304
305/********************************************************************/
306/********************************************************************/
307/********************************************************************/
308
309
310AliGoodTracksPP::~AliGoodTracksPP()
311{
312 delete [] fGoodInEvent;
313 for (Int_t i = 0;i<fNevents;i++)
314 delete [] fData[i];
315 delete [] fData;
316}
317/********************************************************************/
318AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
319{
320
321 fNevents = 100;
322 ifstream in(infilename.Data());
323
324 if(!in)
325 {
326 cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
327 delete this;
328 return;
329 }
330
331
332 fGoodInEvent = new Int_t[fNevents];
333 fData = new struct GoodTrack* [fNevents];
334
335 Int_t i;
336 for( i = 0;i<fNevents;i++)
337 {
338 fGoodInEvent[i] =0;
339 fData[i] = new struct GoodTrack[500];
340 }
341 Float_t tmp;
342 Int_t evno;
343 while(in>>evno)
344 {
345 if(fGoodInEvent[evno]>=500)
346 {
347 cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
348 continue;
349 }
350 in>>fData[evno][fGoodInEvent[evno]].lab;
351 in>>fData[evno][fGoodInEvent[evno]].code;
352 in>>fData[evno][fGoodInEvent[evno]].px;
353 in>>fData[evno][fGoodInEvent[evno]].py;
354 in>>fData[evno][fGoodInEvent[evno]].pz;
355 in>>fData[evno][fGoodInEvent[evno]].x;
356 in>>fData[evno][fGoodInEvent[evno]].y;
357 in>>fData[evno][fGoodInEvent[evno]].z;
358
359 in>>tmp;
360 in>>tmp;
361 in>>tmp;
362
363 /* cout<<evno<<" ";
364 cout<<fData[evno][fGoodInEvent[evno]].lab;
365 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
366 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
367 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
368 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
369 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
370 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
371 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
372 cout<<"\n";
373 */
374 fGoodInEvent[evno]++;
375 }
376 in.close();
377 cout<<"AliGoodTracksPP::AliGoodTracksPP() .... Done\n";
378}
379
380
381
382const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
383 {
384
385 if( (event>fNevents) || (event<0))
386 {
387 gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
388 }
389 if( (n>fGoodInEvent[event]) || (n<0))
390 {
391 gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
392 }
393 return fData[event][n];
394
395 }