Manager class for merging/digitization
[u/mrichter/AliRoot.git] / STEER / AliRunDigitizer.cxx
CommitLineData
9ce40367 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/*
17$Log$
18*/
19
20////////////////////////////////////////////////////////////////////////
21//
22// AliRunDigitizer.cxx
23//
24// Manager object for merging/digitization
25//
26// Instance of this class manages the digitization and/or merging of
27// Sdigits into Digits.
28//
29// Only one instance of this class is created in the macro:
30// AliRunDigitizer * manager = new AliRunDigitizer();
31// Then instances of specific detector digitizers are created:
32// AliMUONDigitizer *dMUON = new AliMUONDigitizer(manager)
33// and the manager is configured (you have to specify input files
34// and an output file). The manager generates a combination of input
35// events according to an option set by SetCombinationType. Then it
36// connects appropriate trees from the input files, creates TreeD
37// in the output and runs once per event Digitize method of all existing
38// AliDetDigitizers (without any option). AliDetDigitizers ask manager
39// for a TTree with input (manager->GetNextTreeH(TTree *)
40// or manager->GetNextTreeS(TTree *),
41// merge all inputs, digitize it, and save it in the TreeD
42// obtained by manager->GetTreeD(). Output events are stored with
43// numbers from 0, this default can be changed by
44// manager->SetFirstOutputEventNr(Int_t) method. The particle numbers
45// in the output are shifted by MASK, which is taken from manager.
46//
47// Single input file is permitted. Maximum MAXFILESTOMERGE can be merged.
48// Input from the memory (on-the-fly merging) is not yet
49// supported, as well as access to the input data by invoking methods
50// on the output data.
51//
52// Access to the geometrical data is via gAlice for now (supposing the
53// same geometry in all input files), gAlice is taken from the defined
54// input file.
55//
56// Example with MUON digitizer:
57//
58// AliRunDigitizer * manager = new AliRunDigitizer();
59// manager->SetInput("1track_10events_phi45_60.root");
60// manager->SetInput("1track_10events_phi120_135.root");
61// manager->SetOutputDir("/home/z2/jchudoba/batch/jobtmp");
62// manager->SetOutputFile("digits.root");
63// AliMUONDigitizer *dMUON = new AliMUONDigitizer(manager);
64// manager->SetNrOfEventsToWrite(3);
65// manager->SetCopyTreesFromInput(1);
66// manager->Digitize();
67//
68////////////////////////////////////////////////////////////////////////
69
70// system includes
71
72#include <iostream.h>
73
74// ROOT includes
75
76#include "TFile.h"
77#include "TTree.h"
78
79// AliROOT includes
80
81#include "AliRunDigitizer.h"
82#include "AliDigitizer.h"
83#include "AliRun.h"
84#include "AliHeader.h"
85#include "TParticle.h"
86
87ClassImp(AliRunDigitizer)
88
89////////////////////////////////////////////////////////////////////////
90
91 AliRunDigitizer::AliRunDigitizer() : TNamed("AliRunDigitizer","")
92{
93// default ctor
94
95 for (Int_t i=0;i<MAXDETECTORS;i++) fDigitizers[i]=0;
96 fNDigitizers = 0;
97 fNinputs = 0;
98 fOutputFileName = "digits.root";
99 fOutputDirName = "/tmp/";
100 fCombination.Set(MAXFILESTOMERGE);
101 for (Int_t i=0;i<MAXFILESTOMERGE;i++) {
102 fArrayTreeS[i]=fArrayTreeH[i]=fArrayTreeTPCS[i]=NULL;
103 fCombination[i]=-1;
104 }
105 fkMASKSTEP = 10000000;
106 fkMASK[0] = 0;
107 for (Int_t i=0;i<MAXFILESTOMERGE;i++) {
108 fkMASK[i] = fkMASK[i-1] + fkMASKSTEP;
109 }
110 fInputFileNames = new TClonesArray("TObjString",1);
111 fInputFiles = new TClonesArray("TFile",1);
112 fMinNEvents = 99999999;
113 fCombinationType=1;
114 fCombinationFileName = "combination.txt";
115 fOutput = 0;
116 fEvent = 0;
117 fNrOfEventsToWrite = 0;
118 fNrOfEventsWritten = 0;
119 fCopyTreesFromInput = -1;
120 fDebug = 3;
121 if (GetDebug()>2)
122 cerr<<"AliRunDigitizer::AliRunDigitizer() called"<<endl;
123}
124
125////////////////////////////////////////////////////////////////////////
126
127AliRunDigitizer::~AliRunDigitizer() {
128// dtor
129
130 if (fInputFiles) delete fInputFiles;
131 if (fInputFileNames) delete fInputFileNames;
132}
133
134////////////////////////////////////////////////////////////////////////
135void AliRunDigitizer::AddDigitizer(AliDigitizer *digitizer)
136{
137// add digitizer to the list of active digitizers
138
139 if (fNDigitizers >= MAXDETECTORS) {
140 cerr<<"Too many detectors to digitize. Increase value of MAXDETECTORS"
141 <<" constant in AliRunDigitizer.h and recompile or decrease the"
142 <<" the number of detectors"<<endl;
143 } else {
144 fDigitizers[fNDigitizers++] = digitizer;
145 }
146}
147
148////////////////////////////////////////////////////////////////////////
149
150Bool_t AliRunDigitizer::SetInput(char *inputFile)
151{
152// receives the name of the input file. Opens it and stores pointer
153// to it, returns kFALSE if fails
154// if inputFile = MEMORY, uses pointers from gAlice - not yet implemented
155
156// input cannot be output - open READ only
157 TFile *file = new((*fInputFiles)[fNinputs]) TFile(inputFile,"READ");
158
159 if (GetDebug()>2) cerr<<"AliRunDigitizer::SetInput: file = "<<file<<endl;
160 if (!file->IsOpen()) {
161 cerr<<"ERROR: AliRunDigitizer::SetInput: cannot open file "
162 <<inputFile<<endl;
163 return kFALSE;
164 }
165
166// find Header and get nr of events there
167 TTree * te = (TTree *) file->Get("TE") ;
168 if (!te) {
169 cerr<<"ERROR: AliRunDigitizer::SetInput: input file does "
170 <<"not contain TE"<<endl;
171 return kFALSE;
172 }
173 Int_t nEntries = (Int_t) te->GetEntries();
174
175 if (GetDebug()>2) cerr<<"AliRunDigitizer::SetInput: nEntries = "
176 <<nEntries<<endl;
177 if (nEntries < 1) {
178 cerr<<"ERROR: AliRunDigitizer::SetInput: input file does "
179 <<"not contain any event"<<endl;
180 return kFALSE;
181 }
182
183 if (nEntries < fMinNEvents) fNrOfEventsToWrite = fMinNEvents = nEntries;
184
185// find gAlice object if it is a first input file
186// this is only temporary solution, we need gAlice to access detector
187// geometry parameters. Unfortunately I have to include AliRun header file.
188 if (fNinputs == 0) {
189 gAlice = (AliRun*)file->Get("gAlice");
190 if (GetDebug() > 2) cerr<<"gAlice taken from the first input: "
191 <<gAlice<<endl;
192 if (!gAlice) {
193 cerr<<"ERROR: AliRunDigitizer::SetInput: first input file "
194 <<"does not contain gAlice object"<<endl;
195 return kFALSE;
196 }
197 }
198
199// store this file name if it is OK
200 new((*fInputFileNames)[fNinputs]) TObjString(inputFile);
201 fNinputs++;
202 if (GetDebug() > 2) cerr<<"fNinputs = "<<fNinputs<<endl;
203 return kTRUE;
204}
205
206////////////////////////////////////////////////////////////////////////
207Bool_t AliRunDigitizer::MakeCombination()
208{
209// make a new combination of events from different files
210
211 Int_t type = fCombinationType;
212 if (fNrOfEventsWritten >= fNrOfEventsToWrite) return kFALSE;
213
214 switch (type) {
215
216 case 1:
217// type = 1: 1-1-1 - take the same event number from each file
218 if (fCombination[0]<fMinNEvents-1) {
219 for (Int_t i=0;i<fNinputs;i++) fCombination[i]++;
220 return kTRUE;
221 }
222 if (GetDebug()>2)
223 cerr<<"AliRunDigitizer::MakeCombination: Warning: "
224 <<"maximum number of Events in an input file "
225 <<"was reached."<<endl;
226 break;
227
228 case 2:
229// type = 2: read from the file combinations.ascii
230// not yet implemented 100% correctly - requires 4 entries in the row
231 static FILE *fp ;
232 static Int_t linesRead;
233 if (!fp) {
234 fp = fopen(fCombinationFileName.Data(),"r");
235 linesRead = 0;
236 }
237 if (!fp) {
238 cerr<<"AliRunDigitizer::MakeCombination ERROR: "
239 <<"Cannot open input file with combinations."<<endl;
240 return kFALSE;
241 }
242 Int_t var[4], nInputs;
243 char line[80];
244// since I do not close or rewind the file, the position should be correct
245 if (fgets(line,80,fp)) {
246 nInputs = sscanf(&line[0],"%d%d%d%d",&var[0],&var[1],&var[2],&var[3]);
247 if (nInputs != fNinputs) {
248 cerr<<"AliRunDigitizer::MakeCombination ERROR: "
249 <<"Nr. of input files is different from nr "
250 <<"integers defining the combination"<<endl;
251 return kFALSE;
252 }
253 while(nInputs--) {
254 fCombination[nInputs] = var[nInputs];
255 }
256 return kTRUE;
257 } else {
258 cerr<<"AliRunDigitizer::MakeCombination ERROR: "
259 <<"no more input in the file with combinations"<<endl;
260 return kFALSE;
261 }
262
263 default:
264 cerr<<"AliRunDigitizer::MakeCombination: ERROR: "
265 <<"wrong type of required combination type: "<<type<<endl;
266 }
267 return kFALSE;
268}
269
270////////////////////////////////////////////////////////////////////////
271void AliRunDigitizer::PrintCombination()
272{
273// debug method to print current combination
274
275 cerr<<"AliRunDigitizer::PrintCombination: Current events combination:"<<endl;
276 for (Int_t i=0;i<fNinputs;i++) {
277 cerr<<"File: "<<((TObjString *)fInputFileNames->At(i))->GetString()<<"\tEvent: "<<fCombination[i]<<endl;
278 }
279}
280
281////////////////////////////////////////////////////////////////////////
282void AliRunDigitizer::Digitize()
283{
284// get a new combination of inputs, connect input trees and loop
285// over all digitizers
286
287 if (!InitGlobal()) {
288 cerr<<"False from InitGlobal"<<endl;
289 return;
290 }
291// take gAlice from the first input file. It is needed to access
292// geometry data
293 while (MakeCombination()) {
294 if (GetDebug()>2) PrintCombination();
295 ConnectInputTrees();
296 InitPerEvent();
297// loop over all registered digitizers and let them do the work
298 for (Int_t i=0;i<fNDigitizers; i++) {
299 fDigitizers[i]->Digitize();
300 }
301 FinishPerEvent();
302 }
303 FinishGlobal();
304}
305
306////////////////////////////////////////////////////////////////////////
307void AliRunDigitizer::ConnectInputTrees()
308{
309// fill arrays fArrayTreeS, fArrayTreeH and fArrayTreeTPCS with
310// pointers to the correct events according fCombination values
311// null pointers can be in the output, AliDigitizer has to check it
312
313 TFile *file;
314 TTree *tree;
315 char treeName[20];
316 for (Int_t i=0;i<fNinputs;i++) {
317 file = (TFile*)(*fInputFiles)[i];
318 sprintf(treeName,"TreeS%d",fCombination[i]);
319 tree = (TTree *) file->Get(treeName);
320 fArrayTreeS[i] = tree;
321 sprintf(treeName,"TreeH%d",fCombination[i]);
322 tree = (TTree *) file->Get(treeName);
323 fArrayTreeH[i] = tree;
324 sprintf(treeName,"TreeD_75x40_100x60_%d",fCombination[i]);
325 tree = (TTree *) file->Get(treeName);
326 fArrayTreeTPCS[i] = tree;
327 }
328}
329
330////////////////////////////////////////////////////////////////////////
331Bool_t AliRunDigitizer::InitGlobal()
332{
333// called once before Digitize() is called, initialize digitizers and output
334
335 if (!InitOutputGlobal()) return kFALSE;
336 for (Int_t i=0;i<fNDigitizers; i++) {
337 if (!fDigitizers[i]->Init()) return kFALSE;
338 }
339 return kTRUE;
340}
341
342////////////////////////////////////////////////////////////////////////
343Bool_t AliRunDigitizer::InitOutputGlobal()
344{
345// Creates the output file, called once by InitGlobal()
346
347 TString fn;
348 fn = fOutputDirName + '/' + fOutputFileName;
349 fOutput = new TFile(fn,"recreate");
350 if (GetDebug()>2) {
351 cerr<<"file "<<fn.Data()<<" was opened"<<endl;
352 cerr<<"fOutput = "<<fOutput<<endl;
353 }
354 if (fOutput) return kTRUE;
355 Error("InitOutputGlobal","Could not create output file.");
356 return kFALSE;
357}
358
359
360////////////////////////////////////////////////////////////////////////
361void AliRunDigitizer::InitPerEvent()
362{
363// Creates TreeDxx in the output file, called from Digitize() once for
364// each event. xx = fEvent
365
366 if (GetDebug()>2)
367 cerr<<"AliRunDigitizer::InitPerEvent: fEvent = "<<fEvent<<endl;
368 fOutput->cd();
369 char hname[30];
370 sprintf(hname,"TreeD%d",fEvent);
371 fTreeD = new TTree(hname,"Digits");
372// fTreeD->Write(); // Do I have to write it here???
373}
374
375
376////////////////////////////////////////////////////////////////////////
377TTree* AliRunDigitizer::GetNextTreeH(TTree *current) const
378{
379// returns next one after the current one
380// returns the first if the current is NULL
381// returns NULL if the current is the last one
382
383 if (fNinputs <= 0) return 0;
384 if (current == 0) return fArrayTreeH[0];
385 for (Int_t i=0; i<fNinputs-1; i++) {
386 if (current == fArrayTreeH[i]) return fArrayTreeH[i+1];
387 }
388 return 0;
389}
390////////////////////////////////////////////////////////////////////////
391TTree* AliRunDigitizer::GetNextTreeS(TTree *current) const
392{
393// returns next one after the current one
394// returns the first if the current is NULL
395// returns NULL if the current is the last one
396
397 if (fNinputs <= 0) return 0;
398 if (current == 0) return fArrayTreeS[0];
399 for (Int_t i=0; i<fNinputs-1; i++) {
400 if (current == fArrayTreeS[i]) return fArrayTreeS[i+1];
401 }
402 return 0;
403}
404////////////////////////////////////////////////////////////////////////
405TTree* AliRunDigitizer::GetNextTreeTPCS(TTree *current) const
406{
407// returns next one after the current one
408// returns the first if the current is NULL
409// returns NULL if the current is the last one
410
411 if (fNinputs <= 0) return 0;
412 if (current == 0) return fArrayTreeTPCS[0];
413 for (Int_t i=0; i<fNinputs-1; i++) {
414 if (current == fArrayTreeTPCS[i]) return fArrayTreeTPCS[i+1];
415 }
416 return 0;
417}
418
419////////////////////////////////////////////////////////////////////////
420Int_t AliRunDigitizer::GetNextMask(Int_t current) const
421{
422// returns next one after the current one
423// returns the first if the current is negative
424// returns negative if the current is the last one
425
426 if (fNinputs <= 0) return -1;
427 if (current < 0) return fkMASK[0];
428 for (Int_t i=0; i<fNinputs-1; i++) {
429 if (current == fkMASK[i]) return fkMASK[i+1];
430 }
431 return -1;
432}
433////////////////////////////////////////////////////////////////////////
434void AliRunDigitizer::FinishPerEvent()
435{
436// called at the end of loop over digitizers
437
438 fOutput->cd();
439 if (fCopyTreesFromInput > -1) {
440 char treeName[20];
441 Int_t i = fCopyTreesFromInput;
442 sprintf(treeName,"TreeK%d",fCombination[i]);
443 ((TFile*)fInputFiles->At(i))->Get(treeName)->Clone()->Write();
444 sprintf(treeName,"TreeH%d",fCombination[i]);
445 ((TFile*)fInputFiles->At(i))->Get(treeName)->Clone()->Write();
446 }
447 fEvent++;
448 fNrOfEventsWritten++;
449 if (fTreeD) {
450 delete fTreeD;
451 fTreeD = 0;
452 }
453}
454////////////////////////////////////////////////////////////////////////
455void AliRunDigitizer::FinishGlobal()
456{
457// called at the end of Exec
458// save unique objects to the output file
459
460 fOutput->cd();
461 this->Write();
462 if (fCopyTreesFromInput > -1) {
463 ((TFile*)fInputFiles->At(fCopyTreesFromInput))->Get("TE")
464 ->Clone()->Write();
465 gAlice->Write();
466 }
467 fOutput->Close();
468}
469
470
471////////////////////////////////////////////////////////////////////////
472Int_t AliRunDigitizer::GetNParticles(Int_t event)
473{
474// return number of particles in all input files for a given
475// event (as numbered in the output file)
476// return -1 if some file cannot be accessed
477
478 Int_t sum = 0;
479 Int_t sumI;
480 for (Int_t i = 0; i < fNinputs; i++) {
481 sumI = GetNParticles(GetInputEventNumber(event,i), i);
482 if (sumI < 0) return -1;
483 sum += sumI;
484 }
485 return sum;
486}
487
488////////////////////////////////////////////////////////////////////////
489Int_t AliRunDigitizer::GetNParticles(Int_t event, Int_t input)
490{
491// return number of particles in input file input for a given
492// event (as numbered in this input file)
493// return -1 if some error
494
495 TFile *file = ConnectInputFile(input);
496 if (!file) {
497 Error("GetNParticles","Cannot open input file");
498 return -1;
499 }
500
501// find the header and get Nprimaries and Nsecondaries
502 TTree* tE = (TTree *)file->Get("TE") ;
503 if (!tE) {
504 Error("GetNParticles","input file does not contain TE");
505 return -1;
506 }
507 AliHeader* header;
508 header = 0;
509 tE->SetBranchAddress("Header", &header);
510 if (!tE->GetEntry(event)) {
511 Error("GetNParticles","event %d not found",event);
512 return -1;
513 }
514 if (GetDebug()>2) {
515 cerr<<"Nprimary: "<< header->GetNprimary()<<endl;
516 cerr<<"Nsecondary: "<<header->GetNsecondary()<<endl;
517 }
518 return header->GetNprimary() + header->GetNsecondary();
519}
520
521////////////////////////////////////////////////////////////////////////
522Int_t* AliRunDigitizer::GetInputEventNumbers(Int_t event)
523{
524// return pointer to an int array with input event numbers which were
525// merged in the output event event
526
527// simplified for now, implement later
528 Int_t a[MAXFILESTOMERGE];
529 for (Int_t i = 0; i < fNinputs; i++) {
530 a[i] = event;
531 }
532 return a;
533}
534////////////////////////////////////////////////////////////////////////
535Int_t AliRunDigitizer::GetInputEventNumber(Int_t event, Int_t input)
536{
537// return an event number of an eventInput from input file input
538// which was merged to create output event event
539
540// simplified for now, implement later
541 return event;
542}
543////////////////////////////////////////////////////////////////////////
544TFile* AliRunDigitizer::ConnectInputFile(Int_t input)
545{
546// open input file with index input
547// 1st file has index 0
548// return 0x0 if fails
549
550// check if file with index input is already open
551 TFile *file = 0;
552 if (fInputFiles->GetEntriesFast() > input)
553 file = static_cast<TFile *>(fInputFiles->At(input));
554
555 if (!file) {
556 // find the file name and open it
557 TObjString *fn = static_cast<TObjString *>(fInputFileNames->At(input));
558 file = new((*fInputFiles)[input]) TFile((fn->GetString()).Data(),"READ");
559 if (!file) {
560 Error("ConnectInputFile","Cannot open input file");
561 return 0;
562 }
563 if (!file->IsOpen()) {
564 Error("ConnectInputFile","Cannot open input file");
565 return 0;
566 }
567 }
568 return file;
569}
570
571////////////////////////////////////////////////////////////////////////
572TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t event)
573{
574// return pointer to particle with index i (index with mask)
575
576// decode the MASK
577 Int_t input = i/fkMASKSTEP;
578 return GetParticle(i,input,GetInputEventNumber(event,input));
579}
580
581////////////////////////////////////////////////////////////////////////
582TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t input, Int_t event)
583{
584// return pointer to particle with index i in the input file input
585// (index without mask)
586// event is the event number in the file input
587// return 0 i fit does not exist
588
589 TFile *file = ConnectInputFile(input);
590 if (!file) {
591 Error("GetParticle","Cannot open input file");
592 return 0;
593 }
594
595// find the header and get Nprimaries and Nsecondaries
596 TTree* tE = (TTree *)file->Get("TE") ;
597 if (!tE) {
598 Error("GetParticle","input file does not contain TE");
599 return 0;
600 }
601 AliHeader* header;
602 header = 0;
603 tE->SetBranchAddress("Header", &header);
604 if (!tE->GetEntry(event)) {
605 Error("GetParticle","event %d not found",event);
606 return 0;
607 }
608
609// connect TreeK
610 char treeName[30];
611 sprintf(treeName,"TreeK%d",event);
612 TTree* tK = static_cast<TTree*>(file->Get(treeName));
613 if (!tK) {
614 Error("GetParticle","input file does not contain TreeK%d",event);
615 return 0;
616 }
617 TParticle *particleBuffer;
618 particleBuffer = 0;
619 tK->SetBranchAddress("Particles", &particleBuffer);
620
621
622// algorithmic way of getting entry index
623// (primary particles are filled after secondaries)
624 Int_t entry;
625 if (i<header->GetNprimary())
626 entry = i+header->GetNsecondary();
627 else
628 entry = i-header->GetNprimary();
629// tK->GetEntry(0); // do I need this???
630 Int_t bytesRead = tK->GetEntry(entry);
631// new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
632 if (bytesRead)
633 return particleBuffer;
634 return 0;
635}
636////////////////////////////////////////////////////////////////////////
637