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
 101 #include <TObjArray.h>
 102 #include <TAxis.h>
 103 #include <TTreeFormula.h>
 104 #include <TChain.h>
 105 #include <TFile.h>
 106 #include <Riostream.h>
 107 #include <cstring>
 108
 109 ClassImp(AliEventPoolSparse)
 110
 111 // _________________________________________________________________________
 112 AliEventPoolSparse::AliEventPoolSparse() :
 113   fHnSparseI(),
 114   fChunkSize(1024 * 16),
 115   fN(0),
 116   fPool(0x0),
 117   fCurrentBin(-1),
 118   fTagChain(0x0),
 119   fVars(0x0),
 120   fRunCut(0x0),
 121   fLHCCut(0x0),
 122   fDetCut(0x0),
 123   fEvCut(0x0),
 124   fBinNumber(0)
 125 {
 126   // Default constructor. Initializes the THnSparseI,
 127   // the initial size of the array and the array itself
 128   fN = fChunkSize;
 129   fPool = new TEntryList * [fN];
 130   memset(fPool,0x0,fN*sizeof(TEntryList*));
 131 }
 132
 133 // _________________________________________________________________________
 134 AliEventPoolSparse::AliEventPoolSparse(const char* name, const char* title, TChain * tagchain, Int_t dim,
 135                                        const char ** vars, const Int_t* nbins, const Double_t* xmin,
 136                                        const Double_t* xmax, Int_t chunksize):
 137   fHnSparseI(name, title, dim, nbins, xmin, xmax, chunksize),
 138   fChunkSize(chunksize),
 139   fN(0),
 140   fPool(0x0),
 141   fCurrentBin(-1),
 142   fTagChain(tagchain),
 143   fVars(0x0),
 144   fRunCut(0x0),
 145   fLHCCut(0x0),
 146   fDetCut(0x0),
 147   fEvCut(0x0),
 148   fBinNumber(0){
 149   // Constructor. Initializes the THnSparseI,
 150   // the initial size of the pool array and the array itself
 151   // It uses the provided array of variables to create TTreeFormulas
 152   // that are used when the pools are filled. This is the reason to require the input
 153   // tag chain in the constructor.
 154
 155   fN = fChunkSize;
 156   fPool = new TEntryList * [fN];
 157   memset(fPool,0x0,fN*sizeof(TArrayI*));
 158
 159   // Pool variables
 160   fVars = new TTreeFormula*[dim];
 161
 162   for (Int_t ivar=0; ivar<dim; ++ivar) {
 163     fVars[ivar] = new TTreeFormula(vars[ivar],vars[ivar],fTagChain);
 164   }
 165
 166
 167 }
 168
 169 // _________________________________________________________________________
 170 AliEventPoolSparse::~AliEventPoolSparse() {
 171   // Destructor. Delete the pool, the array of TTreeFormula
 172   // and the pointers to cuts
 173   for (Int_t i=0; i<fN; ++i) delete fPool[i];
 174   if (fN>0) delete [] fPool;
 175
 176   Int_t ndim = fHnSparseI.GetNdimensions();
 177   for (Int_t i=0; i<ndim; ++i) delete fVars[i];
 178   delete [] fVars;
 179
 180   delete fRunCut;
 181   delete fLHCCut;
 182   delete fDetCut;
 183   delete fEvCut;
 184 }
 185
 186
 187 // Implementation of the interface functions
 188 // _________________________________________________________________________
 189 TChain* AliEventPoolSparse::GetNextChain(){
 190   // Return the chains one by one. The output is 0x0 if the pool is not initialized
 191   // or the last chain is already reached
 192   if (fCurrentBin<0) {
 193     AliError("The event pool is not initialized");
 194     return 0x0;
 195   }
 196
 197   if (fCurrentBin>=fHnSparseI.GetNbins()) { // Check if >= or >
 198     AliInfo("No more chains");
 199     return 0x0;
 200   }
 201
 202   fBinNumber++;
 203
 204   fChain->SetEntryList(fPool[fCurrentBin++],"ne");
 205   return fChain;
 206 }
 207
 208 // _________________________________________________________________________
 209 void  AliEventPoolSparse::GetCurrentBin(Float_t* xbin) {
 210   // This method fills the center of the current bin in xbin
 211
 212   if (fCurrentBin<0) {
 213     AliError("The event pool is not initialized");
 214     return;
 215   }
 216
 217   Int_t ndim = fHnSparseI.GetNdimensions();
 218   Int_t * coord = new Int_t[ndim];
 219   fHnSparseI.GetBinContent(fCurrentBin,coord);
 220
 221   TObjArray * axes = fHnSparseI.GetListOfAxes();
 222   for (Int_t i=0; i<ndim; ++i)
 223     xbin[i]=((TAxis*)axes->At(i+1))->GetBinCenter(coord[i]);
 224
 225   delete [] coord;
 226 }
 227
 228 // _________________________________________________________________________
 229 void  AliEventPoolSparse::Init(){
 230   // Loop on the tag chain and select the events according
 231   // to the Run, LHC, detector, and event cuts.
 232   // Fill the THnSparse bin and add the event to the corresponding pool
 233   // Taken and modified from AliAnalysisTag
 234
 235   if (!fTagChain) {
 236     AliError("Please provide a tag chain!");
 237     return;
 238   }
 239
 240   Int_t ndim = fHnSparseI.GetNdimensions();
 241   if (ndim<=0) return;
 242
 243   Double_t * x = new Double_t[ndim];
 244
 245   // Tag objects.
 246   AliRunTag *tag = new AliRunTag;
 247   AliEventTag *evTag = new AliEventTag;
 248   fTagChain->SetBranchAddress("AliTAG",&tag);
 249
 250   TString guid("");
 251   TString turl("");
 252   TString path("");
 253
 254   Int_t current = -1;    // Current tree number
 255   for(Int_t iTagFiles = 0; iTagFiles < fTagChain->GetEntries(); iTagFiles++) {
 256     fTagChain->GetEntry(iTagFiles);
 257
 258     if (current != fTagChain->GetTreeNumber()) {
 259       // Update the formula leaves if a new file is processed by the chain
 260       if (fRunCut) fRunCut->UpdateFormulaLeaves();
 261       if (fLHCCut) fLHCCut->UpdateFormulaLeaves();
 262       if (fDetCut) fDetCut->UpdateFormulaLeaves();
 263       if (fEvCut)  fEvCut->UpdateFormulaLeaves();
 264
 265       for (Int_t ivar=0; ivar<fHnSparseI.GetNdimensions(); ++ivar)
 266         if (fVars[ivar]) fVars[ivar]->UpdateFormulaLeaves();
 267
 268       // Create the ESD/AOD chain if not done
 269       if (!fChain) {
 270         // Decide if we have ESD or AOD
 271         TFile * tagfile = fTagChain->GetFile();
 272         if (strstr(tagfile->GetName(),"ESD")) fChain = new TChain("esdTree");
 273         else if (strstr(tagfile->GetName(),"AOD")) fChain = new TChain("aodTree");
 274         else {
 275           AliError("Only ESD and AOD type is implemented!!!");
 276           delete [] x;
 277           return;
 278         }
 279       }
 280
 281       // Update the tree number
 282       current = fTagChain->GetTreeNumber();
 283     }
 284
 285     // Apply Run, LHC, and detector cuts if they exist
 286     if(!fRunCut || fRunCut->EvalInstance(iTagFiles) == 1) {
 287       if(!fLHCCut || fLHCCut->EvalInstance(iTagFiles) == 1) {
 288         if(!fDetCut || fDetCut->EvalInstance(iTagFiles) == 1) {
 289
 290
 291           // Get access to the event data in the TTreeFormula
 292           if (fEvCut) fEvCut->GetNdata();
 293           for (Int_t ivar=0; ivar<fHnSparseI.GetNdimensions(); ++ivar)
 294             if (fVars[ivar]) fVars[ivar]->GetNdata();
 295
 296           // Loop on events
 297           const TClonesArray *tagList = tag->GetEventTags();
 298           Int_t iEvents = tagList->GetEntries();
 299           for(Int_t i = 0; i < iEvents; i++) {
 300             evTag = (AliEventTag *) tagList->At(i);
 301
 302             guid = evTag->GetGUID();
 303             turl = evTag->GetTURL();
 304             path = evTag->GetPath();
 305
 306
 307             if(!fEvCut || fEvCut->EvalInstance(i) == 1) {
 308               TEntryList *fLocalList = new TEntryList();
 309               fLocalList->SetTreeName(fChain->GetName());
 310               fLocalList->SetFileName(turl.Data());
 311               fLocalList->Enter(i);
 312
 313
 314               // Add this event to the corresponding pool
 315               {
 316                 // Increment the bin content corrresponding to the vector "x" by "w",
 317                 // and store the event index iev to the array associated with the bin,
 318                 // then return the bin index.
 319
 320                 for (Int_t ivar=0; ivar<ndim; ++ivar) x[ivar] = fVars[ivar]->EvalInstance(i);
 321
 322                 Int_t bin =  fHnSparseI.Fill(x);
 323                 // Check if we have to enlarge the array of pointers
 324                 if (bin>=fN) Set(bin+fChunkSize);
 325                 // Allocate the TEntryList if this is the first use of it
 326                 if (!fPool[bin]) fPool[bin] = new TEntryList();
 327                 // Add the event iev to the corresponding bin
 328                 fPool[bin]->Add(fLocalList);
 329               }
 330             }
 331           }//event loop
 332
 333           for (Int_t ipool=0; ipool<fHnSparseI.GetNbins(); ++ipool)
 334             fPool[ipool]->OptimizeStorage();
 335
 336           // Add the current file to the ESD/AOD chain
 337           if(!path.IsNull()) fChain->AddFile(path);
 338           else if(!turl.IsNull()) fChain->AddFile(turl);
 339
 340         }//detector tag cuts
 341       }//lhc tag cuts
 342     }//run tag cut
 343   }//tag file loop
 344
 345   delete [] x;
 346   fCurrentBin = 0; // Initialize the current bin
 347 }
 348
 349 // _________________________________________________________________________
 350 void AliEventPoolSparse::SetRunCut(const char * cut){
 351   // Run selection cuts
 352   if (fRunCut) delete fRunCut;
 353   fRunCut = new TTreeFormula("fRun",cut,fTagChain);
 354 }
 355
 356 // _________________________________________________________________________
 357 void AliEventPoolSparse::SetLHCCut(const char * cut){
 358   // LHC selection cuts
 359   if (fLHCCut) delete fLHCCut;
 360   fLHCCut = new TTreeFormula("fLHC",cut,fTagChain);
 361 }
 362
 363 // _________________________________________________________________________
 364 void AliEventPoolSparse::SetDetCut(const char * cut){
 365   // Detector selection cuts
 366   if (fDetCut) delete fDetCut;
 367   fDetCut = new TTreeFormula("fDet",cut,fTagChain);
 368 }
 369
 370 // _________________________________________________________________________
 371 void AliEventPoolSparse::SetEventCut(const char * cut){
 372   // Event selection cuts
 373   if (fEvCut) delete fEvCut;
 374   fEvCut = new TTreeFormula("fEv",cut,fTagChain);
 375 }
 376
 377 // _________________________________________________________________________
 378 void AliEventPoolSparse::Set(Int_t n){
 379   // Set size of the array of pointers to n.
 380   // A new array is created, the old contents copied to the new array,
 381   // then the old array is deleted.
 382   // This function is taken from TArrayI
 383
 384   if (n < 0) return;
 385   if (n != fN) {
 386     TEntryList **temp = fPool;
 387     if (n != 0) {
 388       fPool = new TEntryList*[n];
 389       if (n < fN) memcpy(fPool,temp, n*sizeof(TEntryList*));
 390       else {
 391         memcpy(fPool,temp,fN*sizeof(TEntryList*));
 392         memset(&fPool[fN],0x0,(n-fN)*sizeof(TEntryList*));
 393       }
 394     } else {
 395       fPool = 0x0;
 396     }
 397     if (fN) delete [] temp;
 398     fN = n;
 399   }
 400 }