Update concerning changes in pair cut
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
CommitLineData
1b446896 1#include "AliHBTReaderTPC.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(AliHBTReaderTPC)
21//reader for TPC tracking
22//needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc
23//
24//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
25//Piotr.Skowronski@cern.ch
26
27AliHBTReaderTPC::
28 AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
29 const Char_t* goodtracksfilename,const Char_t* galicefilename):
30 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
31 fGAliceFileName(galicefilename),
32 fGoodTPCTracksFileName(goodtracksfilename)
33{
34 //constructor, only file names are set
35 //Defaults:
36 // trackfilename = "AliTPCtracks.root"
37 // clusterfilename = "AliTPCclusters.root"
38 // goodtracksfilename = "good_tracks_tpc"
39 // galicefilename = "" - this means: Do not open gAlice file -
40 // just leave the global pointer untached
41
42 fParticles = new AliHBTRun();
43 fTracks = new AliHBTRun();
44
45 fTracksFile = 0x0; //files are opened during reading only
46 fClustersFile = 0x0;
47
48 fIsRead = kFALSE;
49}
50/********************************************************************/
51
52AliHBTReaderTPC::~AliHBTReaderTPC()
53 {
54 //desctructor
55 delete fParticles;
56 delete fTracks;
57 }
58/********************************************************************/
59
60AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
61 {
62 //returns Nth event with simulated particles
63 if (!fIsRead) Read(fParticles,fTracks);
64 return fParticles->GetEvent(n);
65 }
66/********************************************************************/
67AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
68 {
69 //returns Nth event with reconstructed tracks
70 if (!fIsRead) Read(fParticles,fTracks);
71 return fTracks->GetEvent(n);
72 }
73/********************************************************************/
74
75Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
76 {
77 //returns number of events of particles
78 if (!fIsRead) Read(fParticles,fTracks);
79 return fParticles->GetNumberOfEvents();
80 }
81
82/********************************************************************/
83Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
84 {
85 //returns number of events of tracks
86 if (!fIsRead) Read(fParticles,fTracks);
87 return fTracks->GetNumberOfEvents();
88 }
89/********************************************************************/
90
91
92Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
93 {
94 //reads data and puts put to the particles and tracks objects
95 //reurns 0 if everything is OK
96 //
97 Int_t i; //iterator and some temprary values
98 Int_t Nevents;
99 if (!particles) //check if an object is instatiated
100 {
101 Error("AliHBTReaderTPC::Read"," particles object must instatiated before passing it to the reader");
102 }
103 if (!tracks) //check if an object is instatiated
104 {
105 Error("AliHBTReaderTPC::Read"," tracks object must instatiated before passing it to the reader");
106 }
107 particles->Reset();//clear runs == delete all old events
108 tracks->Reset();
109
110 if( (i=OpenFiles()) )
111 {
112 Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
113 return i;
114 }
115
116 AliGoodTracks *goodTPCTracks = new AliGoodTracks(fGoodTPCTracksFileName);
117 if (!goodTPCTracks)
118 {
119 Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
120 return 1;
121 }
122
123
124 if (gAlice->TreeE())//check if tree E exists
125 {
126 Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
127 cout<<"Found "<<Nevents<<endl;
128 }
129 else
130 {//if not return an error
131 Error("AliHBTReaderPPprod::Read","Can not find Header tree (TreeE) in gAlice");
132 return 1;
133 }
134
135 fClustersFile->cd();//set cluster file active
136 AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
137 if (!TPCParam)
138 {
139 Error("AliHBTReaderTPC::Read","TPC parameters have not been found !\n");
140 return 1;
141 }
142
143 TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
144 tarray->SetOwner(); //set the ownership of the objects it contains
145 //when array is is deleted or cleared all objects
146 //that it contains are deleted
147
148 for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
149 {
150 cout<<"Reading Event "<<currentEvent<<endl;
151 /**************************************/
152 /**************************************/
153 /**************************************/
154 fTracksFile->cd();//set track file active
155
156 Char_t treename[100];
157 sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
158
159 TTree *tracktree=0;
160
161 tracktree=(TTree*)fTracksFile->Get(treename);//get the tree
162 if (!tracktree) //check if we got the tree
163 {//if not return with error
164 Error("AliHBTReaderTPC::Read","Can't get a tree with TPC tracks !\n");
165 return 1;
166 }
167
168 TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
169 if (!trackbranch) ////check if we got the branch
170 {//if not return with error
171 Error("AliHBTReaderTPC::Read","Can't get a branch with TPC tracks !\n");
172 return 2;
173 }
174 Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
175 cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
176 //Copy tracks to array
177
178 AliTPCtrack *iotrack=0;
179
180 fClustersFile->cd();//set cluster file active
181 AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
182 if (!tracker) //check if it has created succeffuly
183 {//if not return with error
184 Error("AliHBTReaderTPC::Read","Can't get a tracker !\n");
185 return 3;
186 }
187 tracker->LoadInnerSectors();
188 tracker->LoadOuterSectors();
189
190 for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
191 {
192 iotrack=new AliTPCtrack; //create new tracks
193 trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
194 tracktree->GetEvent(i); //stream track i to the iotrack
195 tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
196 //which is the label of corresponding simulated particle
197 tarray->AddLast(iotrack); //put the track in the array
198 }
199
200 fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
201 fTracksFile->Delete("tracks");//delete branch from memmory
202 delete tracker; //delete tracker
203
204 tracker = 0x0;
205 trackbranch = 0x0;
206 tracktree = 0x0;
207
208 Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent]; //number of good tracks in the current event
209
210 Double_t xk;
211 Double_t par[5];
212 Float_t phi, lam, pt;//angles and transverse momentum
213 Int_t label; //label of the current track
214 Bool_t found; //flag indicated wether we managed to match good_tpc_track with track
215
216 for (i=0; i<ngood; i++) //loop over all good tracks
217 {
218 const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i); //get ith goog track
219
220 if(Pass(gt.code)) continue; //check if we are intersted with particles of this type
221 //if not take next partilce
222
223 label = gt.lab;
224 found = kFALSE; //guard in case we don't find track with such a label
225 for (Int_t j=0;j<NTPCtracks;j++)//lopp over all tpc tracks
226 {
227 iotrack = (AliTPCtrack*)tarray->At(j);
228 if (iotrack->GetLabel() == label) //if the label is the same
229 {
230 found = kTRUE; //we found the track
231 break;
232 }
233 }
234 if(!found) //check if we found the track
235 {
236 Warning("Read",
237 "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
238 continue; //put comunicate on the screen and continue loop
239 }
240
241 Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();//CMS mass of this particle
242 Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass); //particle total energy
243
244 AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
245 if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
246 //if it does not delete it and take next good track
247
248 iotrack->PropagateTo(gt.x);
249 iotrack->GetExternalParameters(xk,par); //get properties of the track
250 phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
251 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
252 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
253 lam=par[3];
254 pt=1.0/TMath::Abs(par[4]);
255
256 Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
257 Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
258 Double_t tpz = pt * lam; //track z coordinate of momentum
259
260 Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
261
262 AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
263 if(Pass(track)) { delete track;continue;}//check if meets all criteria of any of our cuts
264 //if it does not delete it and take next good track
265
266 particles->AddParticle(currentEvent,part);//put track and particle on the run
267 tracks->AddParticle(currentEvent,track);
268
269 }
270 tarray->Clear(); //clear the array
271
272 /**************************************/
273 /**************************************/
274 /**************************************/
275 }
276
277 //save environment (resouces) --
278 //clean your place after the work
279 CloseFiles();
280 delete tarray;
281 delete goodTPCTracks;
282 fIsRead = kTRUE;
283 return 0;
284 }
285
286/********************************************************************/
287Int_t AliHBTReaderTPC::OpenFiles()
288{
289 //opens all the files
290 fTracksFile = 0;
291 fTracksFile=TFile::Open(fTrackFileName.Data());
292 if (!fTracksFile->IsOpen())
293 {
294 Error("AliHBTReaderTPC::OpenFiles","Can't open file with tacks named ",fTrackFileName.Data());
295 return 1;
296 }
297
298 fClustersFile = 0;
299
300 fClustersFile=TFile::Open(fClusterFileName.Data());
301 if (!fClustersFile->IsOpen())
302 {
303 Error("AliHBTReaderTPC::OpenFiles","Can't open file with TPC clusters named ",fClusterFileName.Data());
304 return 2;
305 }
306
307 return 0;
308}
309
310
311
312/********************************************************************/
313
314void AliHBTReaderTPC::CloseFiles()
315{
316 //closes the files
317 fTracksFile->Close();
318 fClustersFile->Close();
319}
320
321/********************************************************************/
322
323/********************************************************************/
324/********************************************************************/
325/********************************************************************/
326
327
328AliGoodTracks::~AliGoodTracks()
329{
330//destructor
331 delete [] fGoodInEvent;
332 for (Int_t i = 0;i<fNevents;i++)
333 delete [] fData[i];
334 delete [] fData;
335}
336/********************************************************************/
337AliGoodTracks::AliGoodTracks(const TString& infilename)
338{
339
340 cout<<"AliGoodTracks::AliGoodTracks() ....\n";
341 if(!gAlice)
342 {
343 cerr<<"There is no gAlice"<<endl;
344 delete this;
345 return;
346 }
347
348 if (!gAlice->TreeE())
349 {
350 cerr<<"Can not find Header tree (TreeE) in gAlice"<<endl;
351 delete this;
352 return;
353 }
354
355 fNevents = (Int_t)gAlice->TreeE()->GetEntries();
356 //fNevents = 100;
357 cout<<fNevents<<" FOUND"<<endl;
358 ifstream in(infilename.Data());
359
360 if(!in)
361 {
362 cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
363 delete this;
364 return;
365 }
366
367
368 fGoodInEvent = new Int_t[fNevents];
369 fData = new struct GoodTrack* [fNevents];
370
371 Int_t i;
372 for( i = 0;i<fNevents;i++)
373 {
374 fGoodInEvent[i] =0;
375 fData[i] = new struct GoodTrack[50000];
376 }
377
378 Int_t evno;
379 while(in>>evno)
380 {
381 if(fGoodInEvent[evno]>=50000)
382 {
383 cerr<<"AliGoodTracks::AliGoodTracks() : Not enough place in the array\n";
384 continue;
385 }
386 in>>fData[evno][fGoodInEvent[evno]].lab;
387 in>>fData[evno][fGoodInEvent[evno]].code;
388 in>>fData[evno][fGoodInEvent[evno]].px;
389 in>>fData[evno][fGoodInEvent[evno]].py;
390 in>>fData[evno][fGoodInEvent[evno]].pz;
391 in>>fData[evno][fGoodInEvent[evno]].x;
392 in>>fData[evno][fGoodInEvent[evno]].y;
393 in>>fData[evno][fGoodInEvent[evno]].z;
394
395 /* cout<<evno<<" ";
396 cout<<fData[evno][fGoodInEvent[evno]].lab;
397 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
398 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
399 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
400 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
401 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
402 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
403 cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
404 cout<<"\n";
405 */
406 fGoodInEvent[evno]++;
407 }
408 in.close();
409 cout<<"AliGoodTracks::AliGoodTracks() .... Done\n";
410}
411
412
413
414const GoodTrack& AliGoodTracks::GetTrack(Int_t event, Int_t n) const
415 {
416
417 if( (event>fNevents) || (event<0))
418 {
419 gAlice->Fatal("AliGoodTracks::GetTrack","No such Event %d",event);
420 }
421 if( (n>fGoodInEvent[event]) || (n<0))
422 {
423 gAlice->Fatal("AliGoodTracks::GetTrack","No such Good TPC Track %d",n);
424 }
425 return fData[event][n];
426
427 }