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 }