Bug correction. (C.Finck)
[u/mrichter/AliRoot.git] / ANALYSIS / AliReaderAOD.cxx
1 #include "AliReaderAOD.h"
2 //______________________________________________________________________________
3 ////////////////////////////////////////////////////////////////////////////////
4 //                                                                            //
5 // class AliReaderAOD                                                         //
6 //                                                                            //
7 // Reader and Writer for AOD format.                                          //
8 // AODs are stored in a tree named by the variable fgkTreeName.               //
9 // There is stored 1 or 2 branches. Each of them stores AOD objects           //
10 // First branch is named by the variable fgkReconstructedDataBranchName       //
11 // ("reconstructed.") and keeps reconstructed data.                           //
12 // Second branch is called by the variable fgkSimulatedDataBranchName         //
13 // ("simulated.") and stores Monte carlo truth. If both branches are present  //
14 // AODs are parallel, i.e. nth particle in one branch corresponds to the nth  //
15 // particle in the other one.                                                 //
16 //                                                                            //
17 // Since we accept different formats of particles that are stored in AODs     //
18 // reader must take care of that fact: clean buffer if the next file contains //
19 // different particle type.                                                   //
20 //                                                                            //
21 // Piotr.Skowronski@cern.ch                                                   //
22 //                                                                            //
23 ////////////////////////////////////////////////////////////////////////////////
24
25 ClassImp(AliReaderAOD)
26
27 #include <TError.h>
28 #include <TFile.h>
29 #include <TTree.h>
30 #include <TH1.h>
31 #include "AliAOD.h"
32
33
34 const TString AliReaderAOD::fgkTreeName("TAOD");
35 const TString AliReaderAOD::fgkReconstructedDataBranchName("reconstructed.");
36 const TString AliReaderAOD::fgkSimulatedDataBranchName("simulated.");
37
38 AliReaderAOD::AliReaderAOD(const Char_t* aodfilename):
39  fFileName(aodfilename),
40  fReadSim(kFALSE),
41  fReadRec(kTRUE),
42  fTree(0x0),
43  fFile(0x0),
44  fSimBuffer(0x0),
45  fRecBuffer(0x0)
46 {
47   //ctor
48 }
49 /********************************************************************/
50
51 AliReaderAOD::~AliReaderAOD()
52 {
53 //dtor
54   if (fEventSim == fSimBuffer )
55    {
56     fEventSim = 0x0;
57     fEventRec = 0x0;
58    }
59   delete fSimBuffer;
60   delete fRecBuffer;
61   
62   delete fTree;
63   delete fFile;
64 }
65 /********************************************************************/
66
67 void AliReaderAOD::Rewind()
68 {
69 //Rewinds reading
70   delete fTree;
71   fTree = 0x0;
72   delete fFile;
73   fFile = 0x0;
74   fCurrentDir = 0;
75   fNEventsRead= 0;
76 }
77 /********************************************************************/
78 Int_t AliReaderAOD::ReadNext()
79 {
80 //Reads next event
81   
82   Info("ReadNext","Entered");
83   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
84     {
85       if (fFile == 0x0)
86        {
87          Int_t openfailed = OpenFile(fCurrentDir);//rl is opened here
88          if (openfailed)
89            {
90              //Error("ReadNext","Error Occured while opening directory number %d",fCurrentDir);
91              fCurrentDir++;
92              continue;
93            }
94          fCurrentEvent = 0;
95        }
96       //Tree must exist because OpenFile would reuturn error in the other case
97       if ( fCurrentEvent >= fTree->GetEntries() )
98        {
99          delete fTree;
100          fTree = 0x0;
101          delete fFile;
102          fFile = 0x0;
103          fSimBuffer = 0x0;
104          fRecBuffer = 0x0;
105          fCurrentDir++;
106          continue;
107        }
108        
109       Info("ReadNext","Getting event %d",fCurrentEvent);
110       fTree->GetEvent(fCurrentEvent);
111       Info("ReadNext","Getting event %d Done",fCurrentEvent);
112       
113       Int_t retval = 0;
114       if (fReadRec && fReadSim)
115        {
116          retval = ReadRecAndSim();
117        }
118       else
119        {
120          if (fReadRec) retval = ReadRec();
121          if (fReadSim) retval = ReadSim();
122        } 
123
124       fCurrentEvent++;
125       if (retval == 0) fNEventsRead++;
126
127       return retval;//success -> read one event
128       
129     }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
130    
131   return 1; //no more directories to read
132   
133   
134 }
135 /********************************************************************/
136
137 Int_t AliReaderAOD::ReadRecAndSim()
138 {
139 //Reads raconstructed and simulated data 
140
141  Info("ReadRecAndSim","Found %d reconstructed tracks and %d simulated particles",
142        fRecBuffer->GetNumberOfParticles(),fSimBuffer->GetNumberOfParticles());
143
144  if (fCuts->GetEntriesFast() == 0x0)
145   {//if there is no cuts we return pointer to the buffer
146     if (fEventRec != fRecBuffer)
147      {
148        delete fEventRec;
149        delete fEventSim;
150      }
151     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
152     fEventSim = fSimBuffer;
153   }
154  else
155   {//if there are cuts specified
156     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
157      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
158        fEventRec = new AliAOD();
159        fEventSim = new AliAOD();
160
161        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
162        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
163      } 
164     else
165      {//or simply reset it in case it already exists
166        fEventRec->Reset();
167        fEventSim->Reset();
168      }
169
170     Int_t npart = fRecBuffer->GetNumberOfParticles();
171     for (Int_t i = 0; i < npart; i++)
172      {
173        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
174        if (Rejected(prec)) continue;//we make cuts only on simulated data
175
176        fEventRec->AddParticle(prec);
177        fEventSim->AddParticle( fSimBuffer->GetParticle(i));
178      }
179   }
180
181  Info("ReadRecAndSim","Read %d reconstructed tracks and %d simulated particles",
182        fEventRec->GetNumberOfParticles(),fEventSim->GetNumberOfParticles());
183  
184  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
185  
186  return 0;
187 }
188 /********************************************************************/
189
190 Int_t AliReaderAOD::ReadRec()
191 {
192 //Reads reconstructed data only
193
194  Info("ReadRec","Found %d reconstructed tracks",fRecBuffer->GetNumberOfParticles());
195
196  if (fCuts->GetEntriesFast() == 0x0)
197   {//if there is no cuts we return pointer to the buffer
198     if (fEventRec != fRecBuffer)
199      {
200        delete fEventRec;
201      }
202     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
203   }
204  else
205   {//if there are cuts specified
206     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
207      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
208        fEventRec = new AliAOD();
209
210        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
211      } 
212     else
213      {//or simply reset it in case it already exists
214        fEventRec->Reset();
215      }
216
217     Int_t npart = fRecBuffer->GetNumberOfParticles();
218     for (Int_t i = 0; i < npart; i++)
219      {
220        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
221        if (Rejected(prec)) continue;//we make cuts only on simulated data
222
223        fEventRec->AddParticle(prec);
224      }
225   }
226
227  Info("ReadRec","Read %d reconstructed tracks",fEventRec->GetNumberOfParticles());
228  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
229
230  return 0;
231 }
232 /********************************************************************/
233
234 Int_t AliReaderAOD::ReadSim()
235 {
236 //Reads simulated data only
237
238  Info("ReadSim","Found %d simulated particles",fSimBuffer->GetNumberOfParticles());
239
240  if (fCuts->GetEntriesFast() == 0x0)
241   {//if there is no cuts we return pointer to the buffer
242     if (fEventSim != fSimBuffer)
243      {
244        delete fEventSim;
245      }
246     fEventSim = fSimBuffer;
247   }
248  else
249   {//if there are cuts specified
250     if ( (fEventSim == 0x0) || (fEventSim == fSimBuffer) )
251      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
252        fEventSim = new AliAOD();
253
254        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
255      } 
256     else
257      {//or simply reset it in case it already exists
258        fEventSim->Reset();
259      }
260
261     Int_t npart = fSimBuffer->GetNumberOfParticles();
262     for (Int_t i = 0; i < npart; i++)
263      {
264        AliVAODParticle* prec = fSimBuffer->GetParticle(i);
265        if (Rejected(prec)) continue;//we make cuts only on simulated data
266        fEventSim->AddParticle(prec);
267      }
268   }
269
270  Info("ReadSim","Read %d simulated particles",fEventSim->GetNumberOfParticles());
271  fTrackCounter->Fill(fEventSim->GetNumberOfParticles());
272
273
274  return 0;
275 }
276 /********************************************************************/
277
278 Int_t AliReaderAOD::OpenFile(Int_t n)
279 {
280 //opens fFile with  tree
281
282 // Info("ReadNext","Opening File %d",n);
283  const TString dirname = GetDirName(n);
284  if (dirname == "")
285   {
286     if (AliVAODParticle::GetDebug() > 2 )
287       {
288         Info("OpenFile","Got empty string as a directory name."); 
289       }
290    return 1;
291   }
292  
293  TString filename = dirname +"/"+ fFileName;
294  fFile = TFile::Open(filename.Data()); 
295  if ( fFile == 0x0)
296   {
297     Error("OpenFile","Can't open fFile %s",filename.Data());
298     return 2;
299   }
300  if (!fFile->IsOpen())
301   {
302     Error("OpenFile","Can't open fFile  %s",filename.Data());
303     delete fFile;
304     fFile = 0x0;
305     return 3;
306   }
307
308  Info("ReadNext","File Is Opened, Getting the TREE");
309  
310  fTree = dynamic_cast<TTree*>(fFile->Get(fgkTreeName));
311  if (fTree == 0x0)
312   {
313     if (AliVAODParticle::GetDebug() > 2 )
314       {
315         Info("ReadNext","Can not find TTree object named %s",fgkTreeName.Data());
316       }
317     delete fFile;
318     fFile = 0x0;
319     return 4;
320   }
321
322 //  Info("ReadNext","Got TREE, Setting branch addresses");
323
324   if (fReadRec)
325    {
326      TBranch* branch = fTree->GetBranch(fgkReconstructedDataBranchName);
327      if (branch == 0x0)
328       {
329         Error("OpenFile","Can not find branch %s in file %s",
330               fgkReconstructedDataBranchName.Data(),filename.Data());
331         
332         delete fTree;
333         fTree = 0x0;
334         delete fFile;
335         fFile = 0x0;
336         return 5;
337       }
338      fTree->SetBranchAddress(fgkReconstructedDataBranchName,&fRecBuffer);
339    }
340   
341
342   if (fReadSim)
343    {
344      TBranch* branch = fTree->GetBranch(fgkSimulatedDataBranchName);
345      if (branch == 0x0)
346       {
347         Error("OpenFile","Can not find branch %s in file %s",
348               fgkSimulatedDataBranchName.Data(),filename.Data());
349         
350         delete fTree;
351         fTree = 0x0;
352         delete fFile;
353         fFile = 0x0;
354         return 6;
355       }
356      fTree->SetBranchAddress(fgkSimulatedDataBranchName,&fSimBuffer);
357    }
358 //  Info("ReadNext","Got TREE, Addresses are set.");
359 //  Info("ReadNext","Quitting the method.");
360   
361   return 0;
362  
363 }
364 /********************************************************************/
365
366 Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, const char* pclassname,  Bool_t /*multcheck*/)
367 {
368 //reads tracks from runs and writes them to file
369   ::Info("AliReaderAOD::Write","________________________________________________________");
370   ::Info("AliReaderAOD::Write","________________________________________________________");
371   ::Info("AliReaderAOD::Write","________________________________________________________");
372   
373   if (reader == 0x0)
374    {
375      ::Error("AliReaderAOD::Write","Input Reader is NULL");
376      return -1;
377    }
378   TFile *outfile = TFile::Open(outfilename,"recreate");
379   if (outfile == 0x0)
380    {
381      ::Error("AliReaderAOD::Write","Can not open output file %s",outfilename);
382      return -1;
383    }
384
385   TTree *tree = new TTree(fgkTreeName,"Tree with tracks");
386   
387   TBranch *recbranch = 0x0, *simbranch = 0x0;
388   
389   AliAOD* eventrec = new AliAOD();//must be created before Branch is called. Otherwise clones array is not splitted
390   AliAOD* eventsim = new AliAOD();//AOD together with fParticles clones array knowing exact type of particles
391   
392   eventrec->SetParticleClassName(pclassname);
393   eventsim->SetParticleClassName(pclassname);
394  
395   AliAOD* recbuffer = eventrec;
396   AliAOD* simbuffer = eventsim;
397   
398   if (reader->ReadsRec()) recbranch = tree->Branch(fgkReconstructedDataBranchName,"AliAOD",&recbuffer,32000,99);
399   if (reader->ReadsSim()) simbranch = tree->Branch(fgkSimulatedDataBranchName,"AliAOD",&simbuffer,32000,99);
400
401   reader->Rewind();
402   while (reader->Next() == kFALSE)
403    {
404      
405      if (reader->ReadsRec())
406       {//here we can get AOD that has different particle type
407         AliAOD* event = reader->GetEventRec();
408         if ( eventrec->GetParticleClass() != event->GetParticleClass() )
409          {//if class type is not what what we whant we copy particles
410            eventrec->CopyData(event);
411            recbuffer = eventrec;
412          }
413         else
414          {//else just pointer to event from input reader is passed
415            recbuffer = event;
416          } 
417         recbuffer->GetParticle(0)->Print();
418       }
419
420      if (reader->ReadsSim())
421       {
422         AliAOD* event = reader->GetEventSim();
423         if ( eventsim->GetParticleClass() != event->GetParticleClass() )
424          {//if class type is not what what we whant we copy particles
425            eventsim->CopyData(event);
426            simbuffer = eventrec;
427          }
428         else
429          {//else just pointer to event from input reader is passed
430            simbuffer = event;
431          } 
432         simbuffer->GetParticle(0)->Print();
433       }
434      tree->Fill();
435    }
436   
437   ::Info("AliReaderAOD::Write","Written %d events",tree->GetEntries());
438   outfile->cd();
439   tree->Write();
440
441   delete eventsim;
442   delete eventrec;
443   
444   delete tree;
445   delete outfile;
446   return 0; 
447 }
448