Serious memory leak corrected
[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          
104          delete fSimBuffer;
105          delete fRecBuffer;
106          
107          fSimBuffer = 0x0;
108          fRecBuffer = 0x0;
109          fCurrentDir++;
110          continue;
111        }
112       
113       Info("ReadNext","Getting event %d",fCurrentEvent);
114       fTree->GetEvent(fCurrentEvent);
115       Info("ReadNext","Getting event %d Done",fCurrentEvent);
116       
117       Int_t retval = 0;
118       if (fReadRec && fReadSim)
119        {
120          retval = ReadRecAndSim();
121        }
122       else
123        {
124          if (fReadRec) retval = ReadRec();
125          if (fReadSim) retval = ReadSim();
126        } 
127
128       fCurrentEvent++;
129       if (retval != 0) 
130         {
131           //something wrong has happend during reading this event, take next
132           continue;
133         }
134
135       fNEventsRead++;
136       return retval;//success -> read one event
137       
138     }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
139    
140   return 1; //no more directories to read
141   
142   
143 }
144 /********************************************************************/
145
146 Int_t AliReaderAOD::ReadRecAndSim()
147 {
148 //Reads raconstructed and simulated data 
149
150  Info("ReadRecAndSim","Found %d reconstructed tracks and %d simulated particles",
151        fRecBuffer->GetNumberOfParticles(),fSimBuffer->GetNumberOfParticles());
152
153  if (fCuts->GetEntriesFast() == 0x0)
154   {//if there is no cuts we return pointer to the buffer
155     if (fEventRec != fRecBuffer)
156      {
157        delete fEventRec;
158        delete fEventSim;
159      }
160     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
161     fEventSim = fSimBuffer;
162   }
163  else
164   {//if there are cuts specified
165     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
166      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
167        fEventRec = new AliAOD();
168        fEventSim = new AliAOD();
169
170        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
171        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
172      } 
173     else
174      {//or simply reset it in case it already exists
175        fEventRec->Reset();
176        fEventSim->Reset();
177      }
178
179     Int_t npart = fRecBuffer->GetNumberOfParticles();
180     
181     if (npart != fSimBuffer->GetNumberOfParticles())
182      {
183        Error("ReadRecAndSim","There is different number of simulated and reconstructed particles!",
184                               fSimBuffer->GetNumberOfParticles(),npart);
185        return 1;
186      } 
187     for (Int_t i = 0; i < npart; i++)
188      {
189        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
190        AliVAODParticle* psim = fSimBuffer->GetParticle(i);
191        
192        if (prec == 0x0)
193         {
194           Error("ReadRecAndSim","Reconstructed Particle is NULL !!!");
195           continue;
196         }
197        if (psim == 0x0)
198         {
199           Error("ReadRecAndSim","Simulated Particle is NULL !!!");
200           continue;
201         }
202        
203        if (Rejected(prec)) continue;//we make cuts only on reconstructed data
204        
205        fEventRec->AddParticle(prec);
206        fEventSim->AddParticle( fSimBuffer->GetParticle(i));
207      }
208   }
209
210  Info("ReadRecAndSim","Read %d reconstructed tracks and %d simulated particles",
211        fEventRec->GetNumberOfParticles(),fEventSim->GetNumberOfParticles());
212  
213  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
214  
215  return 0;
216 }
217 /********************************************************************/
218
219 Int_t AliReaderAOD::ReadRec()
220 {
221 //Reads reconstructed data only
222
223  Info("ReadRec","Found %d reconstructed tracks",fRecBuffer->GetNumberOfParticles());
224
225  if (fCuts->GetEntriesFast() == 0x0)
226   {//if there is no cuts we return pointer to the buffer
227     if (fEventRec != fRecBuffer)
228      {
229        delete fEventRec;
230      }
231     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
232   }
233  else
234   {//if there are cuts specified
235     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
236      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
237        fEventRec = new AliAOD();
238
239        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
240      } 
241     else
242      {//or simply reset it in case it already exists
243        fEventRec->Reset();
244      }
245
246     Int_t npart = fRecBuffer->GetNumberOfParticles();
247     for (Int_t i = 0; i < npart; i++)
248      {
249        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
250        if (Rejected(prec)) continue;//we make cuts only on simulated data
251
252        fEventRec->AddParticle(prec);
253      }
254   }
255
256  Info("ReadRec","Read %d reconstructed tracks",fEventRec->GetNumberOfParticles());
257  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
258
259  return 0;
260 }
261 /********************************************************************/
262
263 Int_t AliReaderAOD::ReadSim()
264 {
265 //Reads simulated data only
266
267  Info("ReadSim","Found %d simulated particles",fSimBuffer->GetNumberOfParticles());
268
269  if (fCuts->GetEntriesFast() == 0x0)
270   {//if there is no cuts we return pointer to the buffer
271     if (fEventSim != fSimBuffer)
272      {
273        delete fEventSim;
274      }
275     fEventSim = fSimBuffer;
276   }
277  else
278   {//if there are cuts specified
279     if ( (fEventSim == 0x0) || (fEventSim == fSimBuffer) )
280      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
281        fEventSim = new AliAOD();
282
283        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
284      } 
285     else
286      {//or simply reset it in case it already exists
287        fEventSim->Reset();
288      }
289
290     Int_t npart = fSimBuffer->GetNumberOfParticles();
291     for (Int_t i = 0; i < npart; i++)
292      {
293        AliVAODParticle* prec = fSimBuffer->GetParticle(i);
294        if (Rejected(prec)) continue;//we make cuts only on simulated data
295        fEventSim->AddParticle(prec);
296      }
297   }
298
299  Info("ReadSim","Read %d simulated particles",fEventSim->GetNumberOfParticles());
300  fTrackCounter->Fill(fEventSim->GetNumberOfParticles());
301
302
303  return 0;
304 }
305 /********************************************************************/
306
307 Int_t AliReaderAOD::OpenFile(Int_t n)
308 {
309 //opens fFile with  tree
310
311 // Info("ReadNext","Opening File %d",n);
312  const TString dirname = GetDirName(n);
313  if (dirname == "")
314   {
315     if (AliVAODParticle::GetDebug() > 2 )
316       {
317         Info("OpenFile","Got empty string as a directory name."); 
318       }
319    return 1;
320   }
321  
322  TString filename = dirname +"/"+ fFileName;
323  fFile = TFile::Open(filename.Data()); 
324  if ( fFile == 0x0)
325   {
326     Error("OpenFile","Can't open fFile %s",filename.Data());
327     return 2;
328   }
329  if (!fFile->IsOpen())
330   {
331     Error("OpenFile","Can't open fFile  %s",filename.Data());
332     delete fFile;
333     fFile = 0x0;
334     return 3;
335   }
336
337  Info("ReadNext","File %s Is Opened, Getting the TREE",filename.Data());
338  
339  fTree = dynamic_cast<TTree*>(fFile->Get(fgkTreeName));
340  if (fTree == 0x0)
341   {
342     if (AliVAODParticle::GetDebug() > 2 )
343       {
344         Info("ReadNext","Can not find TTree object named %s",fgkTreeName.Data());
345       }
346     delete fFile;
347     fFile = 0x0;
348     return 4;
349   }
350
351 //  Info("ReadNext","Got TREE, Setting branch addresses");
352
353   if (fReadRec)
354    {
355      TBranch* branch = fTree->GetBranch(fgkReconstructedDataBranchName);
356      if (branch == 0x0)
357       {
358         Error("OpenFile","Can not find branch %s in file %s",
359               fgkReconstructedDataBranchName.Data(),filename.Data());
360         
361         delete fTree;
362         fTree = 0x0;
363         delete fFile;
364         fFile = 0x0;
365         return 5;
366       }
367      fTree->SetBranchAddress(fgkReconstructedDataBranchName,&fRecBuffer);
368    }
369   
370
371   if (fReadSim)
372    {
373      TBranch* branch = fTree->GetBranch(fgkSimulatedDataBranchName);
374      if (branch == 0x0)
375       {
376         Error("OpenFile","Can not find branch %s in file %s",
377               fgkSimulatedDataBranchName.Data(),filename.Data());
378         
379         delete fTree;
380         fTree = 0x0;
381         delete fFile;
382         fFile = 0x0;
383         return 6;
384       }
385      fTree->SetBranchAddress(fgkSimulatedDataBranchName,&fSimBuffer);
386    }
387 //  Info("ReadNext","Got TREE, Addresses are set.");
388 //  Info("ReadNext","Quitting the method.");
389   
390   return 0;
391  
392 }
393 /********************************************************************/
394
395 Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, const char* pclassname,  Bool_t /*multcheck*/)
396 {
397 //reads tracks from runs and writes them to file
398   ::Info("AliReaderAOD::Write","________________________________________________________");
399   ::Info("AliReaderAOD::Write","________________________________________________________");
400   ::Info("AliReaderAOD::Write","________________________________________________________");
401   
402   if (reader == 0x0)
403    {
404      ::Error("AliReaderAOD::Write","Input Reader is NULL");
405      return -1;
406    }
407   TFile *outfile = TFile::Open(outfilename,"recreate");
408   if (outfile == 0x0)
409    {
410      ::Error("AliReaderAOD::Write","Can not open output file %s",outfilename);
411      return -1;
412    }
413
414   TTree *tree = new TTree(fgkTreeName,"Tree with tracks");
415   
416   TBranch *recbranch = 0x0, *simbranch = 0x0;
417   
418   AliAOD* eventrec = new AliAOD();//must be created before Branch is called. Otherwise clones array is not splitted
419   AliAOD* eventsim = new AliAOD();//AOD together with fParticles clones array knowing exact type of particles
420   
421   eventrec->SetParticleClassName(pclassname);
422   eventsim->SetParticleClassName(pclassname);
423  
424   AliAOD* recbuffer = eventrec;
425   AliAOD* simbuffer = eventsim;
426   
427   if (reader->ReadsRec()) recbranch = tree->Branch(fgkReconstructedDataBranchName,"AliAOD",&recbuffer,32000,99);
428   if (reader->ReadsSim()) simbranch = tree->Branch(fgkSimulatedDataBranchName,"AliAOD",&simbuffer,32000,99);
429
430   reader->Rewind();
431   while (reader->Next() == kFALSE)
432    {
433      
434      if (reader->ReadsRec())
435       {//here we can get AOD that has different particle type
436         AliAOD* event = reader->GetEventRec();
437         if ( eventrec->GetParticleClass() != event->GetParticleClass() )
438          {//if class type is not what what we whant we copy particles
439            eventrec->CopyData(event);
440            recbuffer = eventrec;
441          }
442         else
443          {//else just pointer to event from input reader is passed
444            recbuffer = event;
445          } 
446       }
447
448      if (reader->ReadsSim())
449       {
450         AliAOD* event = reader->GetEventSim();
451         if ( eventsim->GetParticleClass() != event->GetParticleClass() )
452          {//if class type is not what what we whant we copy particles
453            eventsim->CopyData(event);
454            simbuffer = eventrec;
455          }
456         else
457          {//else just pointer to event from input reader is passed
458            simbuffer = event;
459          } 
460       }
461      tree->Fill();
462    }
463   
464   ::Info("AliReaderAOD::Write","Written %d events",tree->GetEntries());
465   outfile->cd();
466   tree->Write();
467
468   delete eventsim;
469   delete eventrec;
470   
471   delete tree;
472   delete outfile;
473   return 0; 
474 }
475