c2aa96289b71d0e588d68b76b1c99096c581ff92
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
1 #include "AliHBTReaderTPC.h"
2
3 #include <TTree.h>
4 #include <TFile.h>
5 #include <TParticle.h>
6
7 #include <Varargs.h>
8
9 #include <AliRun.h>
10 #include <AliLoader.h>
11 #include <AliStack.h>
12 #include <AliMagF.h>
13 #include <AliTPCtrack.h>
14 #include <AliTPCParam.h>
15 #include <AliTPCtracker.h>
16 #include <AliTPCLoader.h>
17
18 #include "AliHBTRun.h"
19 #include "AliHBTEvent.h"
20 #include "AliHBTParticle.h"
21 #include "AliHBTParticleCut.h"
22
23 ClassImp(AliHBTReaderTPC)
24 //______________________________________________
25 //
26 // class AliHBTReaderTPC
27 //
28 //reader for TPC tracking
29 //needs galice.root, AliTPCtracks.root, AliTPCclusters.root
30 //
31 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
32 //Piotr.Skowronski@cern.ch
33 AliHBTReaderTPC::AliHBTReaderTPC():
34  fFileName("galice.root"),
35  fRunLoader(0x0),
36  fTPCLoader(0x0),
37  fMagneticField(0.0),
38  fUseMagFFromRun(kTRUE)
39 {
40   //constructor, 
41   //Defaults:
42   //  galicefilename = ""  - this means: Do not open gAlice file - 
43   //                         just leave the global pointer untouched
44   
45 }
46
47 AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
48  fFileName(galicefilename),
49  fRunLoader(0x0),
50  fTPCLoader(0x0),
51  fMagneticField(0.0),
52  fUseMagFFromRun(kTRUE)
53 {
54   //constructor, 
55   //Defaults:
56   //  galicefilename = ""  - this means: Do not open gAlice file - 
57   //                         just leave the global pointer untouched
58   
59 }
60 /********************************************************************/
61 AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
62  AliHBTReader(dirs), 
63  fFileName(galicefilename),
64  fRunLoader(0x0),
65  fTPCLoader(0x0),
66  fMagneticField(0.0),
67  fUseMagFFromRun(kTRUE)
68 {
69   //constructor, 
70   //Defaults:
71   //  galicefilename = ""  - this means: Do not open gAlice file - 
72   //                         just leave the global pointer untached
73   
74 }
75 /********************************************************************/
76
77 AliHBTReaderTPC::~AliHBTReaderTPC()
78 {
79  //desctructor
80    delete fRunLoader;
81 }
82 /********************************************************************/
83  
84 void AliHBTReaderTPC::Rewind()
85 {
86   delete fRunLoader;
87   fRunLoader = 0x0;
88   fCurrentDir = 0;
89   fNEventsRead= 0;
90 }
91 /********************************************************************/
92
93 Int_t AliHBTReaderTPC::ReadNext()
94  {
95  //reads data and puts put to the particles and tracks objects
96  //reurns 0 if everything is OK
97  //
98   Info("Read","");
99
100  
101   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
102   tarray->SetOwner(); //set the ownership of the objects it contains
103                       //when array is is deleted or cleared all objects 
104                       //that it contains are deleted
105
106   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
107   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
108
109   fParticlesEvent->Reset();
110   fTracksEvent->Reset();
111
112   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
113    {
114       
115     if (fRunLoader == 0x0) 
116       if (OpenNextSession()) continue;//directory counter is increased inside in case of error
117
118     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
119      {
120        //read next directory
121        delete fRunLoader;//close current session
122        fRunLoader = 0x0;//assure pointer is null
123        fCurrentDir++;//go to next dir
124        continue;//directory counter is increased inside in case of error
125      }
126      
127     Info("ReadNext","Reading Event %d",fCurrentEvent);
128     
129     fRunLoader->GetEvent(fCurrentEvent);
130     TTree *tracktree = fTPCLoader->TreeT();//get the tree 
131     if (!tracktree) //check if we got the tree
132       {//if not return with error
133          Error("ReadNext","Can't get a tree with TPC tracks !\n"); 
134          continue;
135       }
136    
137     TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
138     if (!trackbranch) ////check if we got the branch
139       {//if not return with error
140         Error("ReadNext","Can't get a branch with TPC tracks !\n"); 
141         continue;
142       }
143     Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
144     Info("ReadNext","Found %d TPC tracks.",ntpctracks);
145     //Copy tracks to array
146
147     AliTPCtrack *iotrack=0;
148
149     for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
150      {
151        iotrack=new AliTPCtrack;   //create new tracks
152        trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
153        tracktree->GetEvent(i); //stream track i to the iotrack
154        tarray->AddLast(iotrack); //put the track in the array
155      }
156
157     Double_t xk;
158     Double_t par[5];
159     Float_t phi, lam, pt;//angles and transverse momentum
160     Int_t label; //label of the current track
161
162     AliStack* stack = fRunLoader->Stack();
163     if (stack == 0x0)
164      {
165        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
166        fCurrentEvent++;
167        continue;
168      }
169     stack->Particles();
170
171     for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
172      { 
173        iotrack = (AliTPCtrack*)tarray->At(i);
174        label = iotrack->GetLabel();
175
176        if (label < 0) continue;
177
178        TParticle *p = (TParticle*)stack->Particle(label);
179
180        if(p == 0x0) continue; //if returned pointer is NULL
181        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
182
183        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
184                                    //if not take next partilce
185
186        AliHBTParticle* part = new AliHBTParticle(*p,i);
187        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
188                                                //if it does not delete it and take next good track
189
190 //       iotrack->PropagateTo(3.,0.0028,65.19);
191 //       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
192        iotrack->PropagateToVertex(50.,0.0353);
193        
194        iotrack->GetExternalParameters(xk,par);     //get properties of the track
195        phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
196        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
197        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
198        lam=par[3]; 
199        pt=1.0/TMath::Abs(par[4]);
200
201        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
202        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
203        Double_t tpz = pt * lam; //track z coordinate of momentum
204
205        Double_t mass = p->GetMass();
206        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
207
208        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
209        if(Pass(track))//check if meets all criteria of any of our cuts
210                     //if it does not delete it and take next good track
211         { 
212           delete track;
213           delete part;
214           continue;
215         }
216        fParticlesEvent->AddParticle(part);
217        fTracksEvent->AddParticle(track);
218       }
219       
220     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
221             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
222             fNEventsRead,fCurrentEvent,fCurrentDir);
223     
224     fCurrentEvent++;
225     fNEventsRead++;
226     delete tarray;
227     return 0;
228    }while(fCurrentDir < GetNumberOfDirs());
229
230   delete tarray;
231   return 1;
232  }
233 /********************************************************************/
234
235 Int_t AliHBTReaderTPC::OpenNextSession()
236 {
237   TString filename = GetDirName(fCurrentDir);
238   if (filename.IsNull())
239    {
240      DoOpenError("Can not get directory name");
241      return 1;
242    }
243   filename = filename +"/"+ fFileName;
244   fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
245   if( fRunLoader == 0x0)
246    {
247      DoOpenError("Can not open session.");
248      return 1;
249    }
250
251   fRunLoader->LoadHeader();
252   if (fRunLoader->GetNumberOfEvents() <= 0)
253    {
254      DoOpenError("There is no events in this directory.");
255      return 1;
256    }
257
258   if (fRunLoader->LoadKinematics())
259    {
260      DoOpenError("Error occured while loading kinematics.");
261      return 1;
262    }
263   fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
264   if ( fTPCLoader == 0x0)
265    {
266      DoOpenError("Exiting due to problems with opening files.");
267      return 1;
268    }
269   Info("OpenNextSession","________________________________________________________");
270   Info("OpenNextSession","Found %d event(s) in directory %s",
271         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
272   Float_t mf;
273   if (fUseMagFFromRun)
274    {
275      fRunLoader->LoadgAlice();
276      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
277      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
278      fRunLoader->UnloadgAlice();
279    }
280   else
281    {
282      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
283      if (fMagneticField == 0x0)
284       {
285         Fatal("OpenNextSession","Magnetic field can not be 0.");
286         return 1;//pro forma
287       }
288      mf = fMagneticField*10.;
289    }
290   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
291
292
293   fRunLoader->CdGAFile();
294   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
295   if (!TPCParam) 
296    {
297     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
298     if (!TPCParam) 
299      { 
300        DoOpenError("TPC parameters have not been found !\n");
301        return 1;
302      }
303    }
304
305   if (fTPCLoader->LoadTracks())
306    {
307      DoOpenError("Error occured while loading TPC tracks.");
308      return 1;
309    }
310   
311   fCurrentEvent = 0;
312   return 0;
313 }
314 /********************************************************************/
315
316 void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
317 {
318   // Does error display and clean-up in case error caught on Open Next Session
319
320    va_list ap;
321    va_start(ap,va_(fmt));
322    Error("OpenNextSession", va_(fmt), ap);
323    va_end(ap);
324    
325    delete fRunLoader;
326    fRunLoader = 0x0;
327    fTPCLoader = 0x0;
328    fCurrentDir++;
329 }
330
331 /********************************************************************/
332 /********************************************************************/
333 /********************************************************************/
334