Merge branch 'master', remote branch 'origin' into TPCdev
[u/mrichter/AliRoot.git] / ANALYSIS / AliEventPoolSparse.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17  
18
19 /*
20   Realisation of an AliVEventPool based on THnSparse
21   Author Peter Hristov
22   Peter.Hristov@cern.ch
23
24   Possible usage: three steps
25   1) Creation of a XML tag collection in aliensh
26
27   aliensh:> find -x charm /alice/sim/PDC_08/LHC08x/180001 tag.root > 180001.xml
28
29   2) Merging the tag files
30
31   TGrid::Connect("alien://");
32   TAlienCollection *collection = TAlienCollection::Open(xmlfile);
33   TGridResult* result = collection->GetGridResult("",0);
34   AliTagCreator *t = new AliTagCreator();
35   t->MergeTags("ESD",result);
36
37   3) Chain the merged tag files and test the event pool
38
39 void testpool(const char * dirname = ".", const char * pattern = "Run180001") {
40
41   gSystem->Load("libANALYSIS.so");
42   gSystem->Load("libANALYSISalice.so");
43
44   // Create a chain
45   TChain * fChain = new TChain("T");
46
47
48   // Chain the tag files in the working directory
49   TString fTagFilename;
50   
51   // Open the working directory
52   void * dirp = gSystem->OpenDirectory(dirname);
53   const char * name = 0x0;
54   // Add all files matching *pattern* to the chain
55   while((name = gSystem->GetDirEntry(dirp))) {
56     cout << name << endl;
57     if (strstr(name,pattern)) { 
58       fTagFilename = dirname;
59       fTagFilename += "/";
60       fTagFilename += name;
61                 
62       fChain->Add(fTagFilename);  
63     }//pattern check
64   }//directory loop
65
66
67   Int_t nruns = fChain->GetEntries();
68
69   cout << nruns << " run(s) found in the tag chain." << endl;
70
71   Int_t dim = 3;
72   const char * vars[] = {"fNumberOfPositiveTracks","fNumberOfNegativeTracks","fPrimaryVertexZ"};
73   Int_t nbins[] = {10,10,10};
74   Double_t xmin[] ={-0.5,-0.5,-20};
75   Double_t xmax[] ={49/5,49.5,20};
76   Int_t chunksize = 100;
77
78   AliEventPoolSparse * pool = 
79     new AliEventPoolSparse("test", "test", fChain, dim, vars, nbins, xmin, xmax, chunksize);
80
81   pool->Init();
82
83   TChain * esdchain = 0x0;
84   Int_t ichain = 0;
85   while (esdchain=pool->GetNextChain()) {
86     cout << "Chain: "<< ichain <<" Events: " << esdchain->GetEntries() << endl;
87     ichain++;
88   }
89
90   delete fChain;
91
92 }
93
94 */
95
96 #include "AliEventPoolSparse.h"
97 #include "AliRunTag.h"
98 #include "AliEventTag.h"
99 #include "AliLog.h"
100 #include "AliRunTagCuts.h"
101 #include "AliLHCTagCuts.h"
102 #include "AliDetectorTagCuts.h"
103 #include "AliEventTagCuts.h"
104
105 #include <TObjArray.h>
106 #include <TAxis.h>
107 #include <TTreeFormula.h>
108 #include <TChain.h>
109 #include <TFile.h>
110 #include <Riostream.h>
111 #include <cstring>
112
113 ClassImp(AliEventPoolSparse)
114
115 // _________________________________________________________________________
116 AliEventPoolSparse::AliEventPoolSparse() :
117   fHnSparseI(),
118   fChunkSize(1024 * 16),
119   fN(0),
120   fPool(0x0),
121   fCurrentBin(-1),
122   fTagChain(0x0),
123   fVars(0x0),
124   fRunCut(0x0),
125   fLHCCut(0x0),
126   fDetCut(0x0),
127   fEvCut(0x0),
128   fRunTagCut(0x0),
129   fEventTagCut(0x0),
130   fDetectorTagCut(0x0),
131   fLHCTagCut(0x0),
132   fBinNumber(0)
133 {
134   // Default constructor. Initializes the THnSparseI,
135   // the initial size of the array and the array itself
136   fN = fChunkSize;
137   fPool = new TEntryList * [fN];
138   memset(fPool,0x0,fN*sizeof(TEntryList*));
139 }
140
141 // _________________________________________________________________________
142 AliEventPoolSparse::AliEventPoolSparse(const char* name, const char* title, TChain * tagchain, Int_t dim,
143                                        const char ** vars, const Int_t* nbins, const Double_t* xmin,
144                                        const Double_t* xmax, Int_t chunksize):
145   fHnSparseI(name, title, dim, nbins, xmin, xmax, chunksize),
146   fChunkSize(chunksize),
147   fN(0),
148   fPool(0x0),
149   fCurrentBin(-1),
150   fTagChain(tagchain),
151   fVars(0x0),
152   fRunCut(0x0),
153   fLHCCut(0x0),
154   fDetCut(0x0),
155   fEvCut(0x0),
156   fRunTagCut(0x0),
157   fEventTagCut(0x0),
158   fDetectorTagCut(0x0),
159   fLHCTagCut(0x0),
160   fBinNumber(0){
161   // Constructor. Initializes the THnSparseI,
162   // the initial size of the pool array and the array itself
163   // It uses the provided array of variables to create TTreeFormulas
164   // that are used when the pools are filled. This is the reason to require the input
165   // tag chain in the constructor.
166
167   fN = fChunkSize;
168   fPool = new TEntryList * [fN];
169   memset(fPool,0x0,fN*sizeof(TArrayI*));
170
171   // Pool variables
172   fVars = new TTreeFormula*[dim];
173
174   for (Int_t ivar=0; ivar<dim; ++ivar) {
175     fVars[ivar] = new TTreeFormula(vars[ivar],vars[ivar],fTagChain);
176   }
177
178
179 }
180
181 // _________________________________________________________________________
182 AliEventPoolSparse::~AliEventPoolSparse() {
183   // Destructor. Delete the pool, the array of TTreeFormula
184   // and the pointers to cuts
185   for (Int_t i=0; i<fN; ++i) delete fPool[i];
186   if (fN>0) delete [] fPool;
187
188   Int_t ndim = fHnSparseI.GetNdimensions();
189   for (Int_t i=0; i<ndim; ++i) delete fVars[i];
190   delete [] fVars;
191
192   delete fRunCut;
193   delete fLHCCut;
194   delete fDetCut;
195   delete fEvCut;
196
197   delete fRunTagCut;
198   delete fEventTagCut;
199   delete fDetectorTagCut;
200   delete fLHCTagCut;
201
202 }
203
204
205 // Implementation of the interface functions
206 // _________________________________________________________________________
207 TChain* AliEventPoolSparse::GetNextChain(){
208   // Return the chains one by one. The output is 0x0 if the pool is not initialized
209   // or the last chain is already reached
210   if (fCurrentBin<0) {
211     AliError("The event pool is not initialized");
212     return 0x0;
213   }
214
215   if (fCurrentBin>=fHnSparseI.GetNbins()) { // Check if >= or >
216     AliInfo("No more chains");
217     return 0x0;
218   }
219
220   fBinNumber++;
221   
222   fChain->SetEntryList(fPool[fCurrentBin++],"ne"); 
223   return fChain;
224 }
225
226 // _________________________________________________________________________
227 void  AliEventPoolSparse::GetCurrentBin(Float_t* xbin) {
228   // This method fills the center of the current bin in xbin
229
230   if (fCurrentBin<0) {
231     AliError("The event pool is not initialized");
232     return;
233   }
234
235   Int_t ndim = fHnSparseI.GetNdimensions();
236   Int_t * coord = new Int_t[ndim]; 
237   fHnSparseI.GetBinContent(fCurrentBin,coord);
238
239   TObjArray * axes = fHnSparseI.GetListOfAxes();
240   for (Int_t i=0; i<ndim; ++i) 
241     xbin[i]=((TAxis*)axes->At(i+1))->GetBinCenter(coord[i]);
242
243   delete [] coord;
244 }
245
246 // _________________________________________________________________________
247 void  AliEventPoolSparse::Init(){
248   // Loop on the tag chain and select the events according
249   // to the Run, LHC, detector, and event cuts.
250   // Fill the THnSparse bin and add the event to the corresponding pool
251   // Taken and modified from AliAnalysisTag
252
253   if (!fTagChain) {
254     AliError("Please provide a tag chain!");
255     return;
256   }
257
258   Int_t ndim = fHnSparseI.GetNdimensions();
259   if (ndim<=0) return;
260
261   Double_t * x = new Double_t[ndim];
262   
263   // Tag objects.
264   AliRunTag *tag = new AliRunTag;        
265   AliEventTag *evTag = 0;  
266   fTagChain->SetBranchAddress("AliTAG",&tag);    
267   
268   TString guid("");      
269   TString turl("");      
270   TString path("");      
271   
272   Int_t current = -1;    // Current tree number
273   for(Int_t iTagFiles = 0; iTagFiles < fTagChain->GetEntries(); iTagFiles++) {   
274     fTagChain->GetEntry(iTagFiles);      
275
276     if (current != fTagChain->GetTreeNumber()) {         
277       // Update the formula leaves if a new file is processed by the chain
278 //       if (fRunCut) fRunCut->UpdateFormulaLeaves();    
279 //       if (fLHCCut) fLHCCut->UpdateFormulaLeaves();    
280 //       if (fDetCut) fDetCut->UpdateFormulaLeaves();    
281 //       if (fEvCut)  fEvCut->UpdateFormulaLeaves();
282
283       if (fEventTagCut) fEventTagCut->InitializeTriggerClasses(tag->GetActiveTriggerClasses());
284
285       for (Int_t ivar=0; ivar<fHnSparseI.GetNdimensions(); ++ivar)
286         if (fVars[ivar]) fVars[ivar]->UpdateFormulaLeaves();
287
288       // Create the ESD/AOD chain if not done
289       if (!fChain) {
290         // Decide if we have ESD or AOD
291         TFile * tagfile = fTagChain->GetFile();
292         if (strstr(tagfile->GetName(),"ESD")) fChain = new TChain("esdTree");
293         else if (strstr(tagfile->GetName(),"AOD")) fChain = new TChain("aodTree");
294         else {
295           AliError("Only ESD and AOD type is implemented!!!");
296           delete [] x;
297           return;
298         }
299       }
300       
301       // Update the tree number
302       current = fTagChain->GetTreeNumber();
303     }
304     
305     // Deprecated use of TTreeFormulas
306 //     // Apply Run, LHC, and detector cuts if they exist
307 //     if(!fRunCut || fRunCut->EvalInstance(iTagFiles) == 1) {   
308 //       if(!fLHCCut || fLHCCut->EvalInstance(iTagFiles) == 1) {         
309 //      if(!fDetCut || fDetCut->EvalInstance(iTagFiles) == 1) {
310          
311
312 //        // Get access to the event data in the TTreeFormula
313 //        if (fEvCut) fEvCut->GetNdata();
314 //        for (Int_t ivar=0; ivar<fHnSparseI.GetNdimensions(); ++ivar)
315 //          if (fVars[ivar]) fVars[ivar]->GetNdata();
316          
317 //        // Loop on events
318 //        //      const TClonesArray *tagList = tag->GetEventTags();
319 //        Int_t iFiles = tag->GetNFiles();
320 //        for (int ifs = 0; ifs<iFiles; ifs++) {
321 //          AliFileTag *eftag = (AliFileTag *) tag->GetFileTag(ifs);
322
323 //          guid = eftag->GetGUID();     
324 //          turl = eftag->GetTURL();     
325 //          path = eftag->GetPath();     
326             
327 //          Int_t iEvents = eftag->GetNEvents();
328 //          for(Int_t i = 0; i < iEvents; i++) {         
329 //            evTag = (AliEventTag *) eftag->GetEventTag(i);     
330               
331               
332 //            if(!fEvCut || fEvCut->EvalInstance(i) == 1) {
333 //              TEntryList *fLocalList = new TEntryList();
334 //              fLocalList->SetTreeName(fChain->GetName());
335 //              fLocalList->SetFileName(turl.Data());
336 //              fLocalList->Enter(i);
337                 
338                 
339 //              // Add this event to the corresponding pool
340 //              {
341 //                // Increment the bin content corrresponding to the vector "x" by "w",
342 //                // and store the event index iev to the array associated with the bin,
343 //                // then return the bin index.
344                   
345 //                for (Int_t ivar=0; ivar<ndim; ++ivar) x[ivar] = fVars[ivar]->EvalInstance(i);
346                   
347 //                Int_t bin =  fHnSparseI.Fill(x);
348 //                // Check if we have to enlarge the array of pointers
349 //                if (bin>=fN) Set(bin+fChunkSize);
350 //                // Allocate the TEntryList if this is the first use of it
351 //                if (!fPool[bin]) fPool[bin] = new TEntryList();
352 //                // Add the event iev to the corresponding bin
353 //                fPool[bin]->Add(fLocalList);
354 //              }
355 //            }
356 //          }//event loop        
357           
358 //          for (Int_t ipool=0; ipool<fHnSparseI.GetNbins(); ++ipool) 
359 //            fPool[ipool]->OptimizeStorage();
360             
361 //          // Add the current file to the ESD/AOD chain
362 //          if(!path.IsNull()) fChain->AddFile(path);    
363 //          else if(!turl.IsNull()) fChain->AddFile(turl);
364 //        }
365 //      }//detector tag cuts
366 //       }//lhc tag cuts
367 //     }//run tag cut    
368
369     // Apply Run, LHC, and detector cuts if they exist
370     if(!fRunTagCut || fRunTagCut->IsAccepted(tag)) {     
371       if(!fLHCTagCut || fLHCTagCut->IsAccepted(tag->GetLHCTag())) {      
372         if(!fDetectorTagCut || fDetectorTagCut->IsAccepted(tag->GetDetectorTags())) {
373          
374 //        // Get access to the event data in the TTreeFormula
375 //        if (fEvCut) fEvCut->GetNdata();
376           for (Int_t ivar=0; ivar<fHnSparseI.GetNdimensions(); ++ivar)
377             if (fVars[ivar]) fVars[ivar]->GetNdata();
378          
379           // Loop on events
380           //      const TClonesArray *tagList = tag->GetEventTags();
381           Int_t iFiles = tag->GetNFiles();
382           for (int ifs = 0; ifs<iFiles; ifs++) {
383             AliFileTag *eftag = (AliFileTag *) tag->GetFileTag(ifs);
384
385             guid = eftag->GetGUID();     
386             turl = eftag->GetTURL();     
387             path = eftag->GetPath();     
388             
389             Int_t iEvents = eftag->GetNEvents();
390             for(Int_t i = 0; i < iEvents; i++) {         
391               evTag = (AliEventTag *) eftag->GetEventTag(i);     
392               
393               
394               if(!fEventTagCut || fEventTagCut->IsAccepted(evTag)) {
395                 TEntryList *fLocalList = new TEntryList();
396                 fLocalList->SetTreeName(fChain->GetName());
397                 fLocalList->SetFileName(turl.Data());
398                 fLocalList->Enter(i);
399                 
400                 
401                 // Add this event to the corresponding pool
402                 {
403                   // Increment the bin content corrresponding to the vector "x" by "w",
404                   // and store the event index iev to the array associated with the bin,
405                   // then return the bin index.
406                   
407                   for (Int_t ivar=0; ivar<ndim; ++ivar) x[ivar] = fVars[ivar]->EvalInstance(i);
408                   
409                   Int_t bin =  fHnSparseI.Fill(x);
410                   // Check if we have to enlarge the array of pointers
411                   if (bin>=fN) Set(bin+fChunkSize);
412                   // Allocate the TEntryList if this is the first use of it
413                   if (!fPool[bin]) fPool[bin] = new TEntryList();
414                   // Add the event iev to the corresponding bin
415                   fPool[bin]->Add(fLocalList);
416                 }
417               }
418             }//event loop        
419           
420             for (Int_t ipool=0; ipool<fHnSparseI.GetNbins(); ++ipool) 
421               fPool[ipool]->OptimizeStorage();
422             
423             // Add the current file to the ESD/AOD chain
424             if(!path.IsNull()) fChain->AddFile(path);    
425             else if(!turl.IsNull()) fChain->AddFile(turl);
426           }
427         }//detector tag cuts
428       }//lhc tag cuts
429     }//run tag cut       
430
431   }//tag file loop       
432
433   delete [] x;
434   fCurrentBin = 0; // Initialize the current bin
435 }
436
437 // _________________________________________________________________________
438 void AliEventPoolSparse::SetRunCut(const char * cut){
439   // Run selection cuts
440   if (fRunCut) delete fRunCut;
441   fRunCut = new TTreeFormula("fRun",cut,fTagChain);
442 }
443
444 // _________________________________________________________________________
445 void AliEventPoolSparse::SetLHCCut(const char * cut){
446   // LHC selection cuts
447   if (fLHCCut) delete fLHCCut;
448   fLHCCut = new TTreeFormula("fLHC",cut,fTagChain);
449 }
450
451 // _________________________________________________________________________
452 void AliEventPoolSparse::SetDetCut(const char * cut){
453   // Detector selection cuts
454   if (fDetCut) delete fDetCut;
455   fDetCut = new TTreeFormula("fDet",cut,fTagChain);
456 }
457
458 // _________________________________________________________________________
459 void AliEventPoolSparse::SetEventCut(const char * cut){
460   // Event selection cuts
461   if (fEvCut) delete fEvCut;
462   fEvCut = new TTreeFormula("fEv",cut,fTagChain);
463 }
464
465 // _________________________________________________________________________
466 void AliEventPoolSparse::SetRunCut(AliRunTagCuts* cut)
467 {
468   fRunTagCut = cut;
469 }
470 // _________________________________________________________________________
471 void AliEventPoolSparse::SetEventCut(AliEventTagCuts* cut)
472 {
473   fEventTagCut = cut;
474 }
475 // _________________________________________________________________________
476 void AliEventPoolSparse::SetDetectorCut(AliDetectorTagCuts* cut)
477 {
478   fDetectorTagCut = cut;
479 }
480 // _________________________________________________________________________
481 void AliEventPoolSparse::SetLHCCut(AliLHCTagCuts* cut)
482 {
483   fLHCTagCut = cut;
484 }
485
486 // _________________________________________________________________________
487 void AliEventPoolSparse::Set(Int_t n){
488   // Set size of the array of pointers to n.
489   // A new array is created, the old contents copied to the new array,
490   // then the old array is deleted.
491   // This function is taken from TArrayI
492
493   if (n < 0) return;
494   if (n != fN) {
495     TEntryList **temp = fPool;
496     if (n != 0) {
497       fPool = new TEntryList*[n];
498       if (n < fN) memcpy(fPool,temp, n*sizeof(TEntryList*));
499       else {
500         memcpy(fPool,temp,fN*sizeof(TEntryList*));
501         memset(&fPool[fN],0x0,(n-fN)*sizeof(TEntryList*));
502       }
503     } else {
504       fPool = 0x0;
505     }
506     if (fN) delete [] temp;
507     fN = n;
508   }
509 }