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