New VZERO tender supply which can be used in order to correct for the LHC clock phase...
[u/mrichter/AliRoot.git] / ANALYSIS / AliEventPoolSparse.cxx
CommitLineData
2af56f28 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
39void 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
109ClassImp(AliEventPoolSparse)
110
111// _________________________________________________________________________
112AliEventPoolSparse::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),
a065d8ed 123 fEvCut(0x0),
124 fBinNumber(0)
125{
2af56f28 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// _________________________________________________________________________
134AliEventPoolSparse::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),
a065d8ed 147 fEvCut(0x0),
148 fBinNumber(0){
2af56f28 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// _________________________________________________________________________
170AliEventPoolSparse::~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// _________________________________________________________________________
189TChain* 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
a065d8ed 202 fBinNumber++;
203
2af56f28 204 fChain->SetEntryList(fPool[fCurrentBin++],"ne");
205 return fChain;
206}
207
208// _________________________________________________________________________
209void 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// _________________________________________________________________________
229void 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// _________________________________________________________________________
350void AliEventPoolSparse::SetRunCut(const char * cut){
351 // Run selection cuts
352 if (fRunCut) delete fRunCut;
353 fRunCut = new TTreeFormula("fRun",cut,fTagChain);
354}
355
356// _________________________________________________________________________
357void AliEventPoolSparse::SetLHCCut(const char * cut){
358 // LHC selection cuts
359 if (fLHCCut) delete fLHCCut;
360 fLHCCut = new TTreeFormula("fLHC",cut,fTagChain);
361}
362
363// _________________________________________________________________________
364void AliEventPoolSparse::SetDetCut(const char * cut){
365 // Detector selection cuts
366 if (fDetCut) delete fDetCut;
367 fDetCut = new TTreeFormula("fDet",cut,fTagChain);
368}
369
370// _________________________________________________________________________
371void AliEventPoolSparse::SetEventCut(const char * cut){
372 // Event selection cuts
373 if (fEvCut) delete fEvCut;
374 fEvCut = new TTreeFormula("fEv",cut,fTagChain);
375}
376
377// _________________________________________________________________________
378void 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}