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