Add special treatment for TPC and TRD, they use different trees than other detectors
[u/mrichter/AliRoot.git] / STEER / AliRunDigitizer.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 /*
17 $Log$
18 Revision 1.8  2001/10/21 18:38:43  hristov
19 Several pointers were set to zero in the default constructors to avoid memory management problems
20
21 Revision 1.7  2001/10/04 15:56:07  jchudoba
22 TTask inheritance
23
24 Revision 1.4  2001/09/19 06:23:50  jchudoba
25 Move some tasks to AliStream and AliMergeCombi classes
26
27 Revision 1.3  2001/07/30 14:04:18  jchudoba
28 correct bug in the initialization
29
30 Revision 1.2  2001/07/28 10:44:32  hristov
31 Loop variable declared once; typos corrected
32
33 Revision 1.1  2001/07/27 12:59:00  jchudoba
34 Manager class for merging/digitization
35
36 */
37
38 ////////////////////////////////////////////////////////////////////////
39 //
40 // AliRunDigitizer.cxx
41 //
42 // Manager object for merging/digitization
43 //
44 // Instance of this class manages the digitization and/or merging of
45 // Sdigits into Digits. 
46 //
47 // Only one instance of this class is created in the macro:
48 //   AliRunDigitizer * manager = 
49 //      new AliRunDigitizer(nInputStreams,SPERB);
50 // where nInputStreams is number of input streams and SPERB is
51 // signals per background variable, which determines how combinations
52 // of signal and background events are generated.
53 // Then instances of specific detector digitizers are created:
54 //   AliMUONDigitizer *dMUON  = new AliMUONDigitizer(manager)
55 // and the I/O configured (you have to specify input files 
56 // and an output file). The manager connects appropriate trees from 
57 // the input files according a combination returned by AliMergeCombi 
58 // classcreates. It creates TreeD in the output and runs once per 
59 // event Digitize method of all existing AliDetDigitizers 
60 // (without any option). AliDetDigitizers ask manager
61 // for a TTree with input (manager->GetInputTreeS(Int_t i),
62 // merge all inputs, digitize it, and save it in the TreeD 
63 // obtained by manager->GetTreeD(). Output events are stored with 
64 // numbers from 0, this default can be changed by 
65 // manager->SetFirstOutputEventNr(Int_t) method. The particle numbers
66 // in the output are shifted by MASK, which is taken from manager.
67 //
68 // The default output is to the signal file (stream 0). This can be 
69 // changed with the SetOutputFile(TString fn)  method.
70 //
71 // Single input file is permitted. Maximum MAXSTREAMSTOMERGE can be merged.
72 // Input from the memory (on-the-fly merging) is not yet 
73 // supported, as well as access to the input data by invoking methods
74 // on the output data.
75 //
76 // Access to the some data is via gAlice for now (supposing the 
77 // same geometry in all input files), gAlice is taken from the first 
78 // input file on the first stream.
79 //
80 // Example with MUON digitizer:
81 //
82 //  AliRunDigitizer * manager = new AliRunDigitizer(2,1);
83 //  manager->SetInputStream(0,"1track_10events_phi45_60.root");
84 //  manager->SetInputStream(1,"1track_10events_phi120_135.root");
85 //  manager->SetOutputFile("digits.root");
86 //  AliMUONDigitizer *dMUON  = new AliMUONDigitizer(manager);
87 //  manager->SetNrOfEventsToWrite(1);
88 //  manager->Exec("");
89 //
90 //////////////////////////////////////////////////////////////////////// 
91
92 // system includes
93
94 #include <iostream.h>
95
96 // ROOT includes
97
98 #include "TFile.h"
99 #include "TTree.h"
100 #include "TList.h"
101
102 // AliROOT includes
103
104 #include "AliRunDigitizer.h"
105 #include "AliDigitizer.h"
106 #include "AliRun.h"
107 #include "AliHeader.h"
108 #include "TParticle.h"
109 #include "AliStream.h"
110 #include "AliMergeCombi.h"
111
112 ClassImp(AliRunDigitizer)
113
114 ////////////////////////////////////////////////////////////////////////
115 AliRunDigitizer::AliRunDigitizer(Int_t nInputStreams, Int_t sperb) : TTask("AliRunDigitizer","The manager for Merging")
116 {
117 // default ctor
118   if (nInputStreams == 0) {
119     Error("AliRunDigitizer","Specify nr of input streams");
120     return;
121   }
122   Int_t i;
123   fNinputs = nInputStreams;
124   fOutputFileName = "";
125   fOutputDirName = ".";
126   fCombination.Set(MAXSTREAMSTOMERGE);
127   for (i=0;i<MAXSTREAMSTOMERGE;i++) {
128     fArrayTreeS[i]=fArrayTreeH[i]=fArrayTreeTPCS[i]=NULL;
129     fCombination[i]=-1;
130   }
131   fkMASKSTEP = 10000000;
132   fkMASK[0] = 0;
133   for (i=1;i<MAXSTREAMSTOMERGE;i++) {
134     fkMASK[i] = fkMASK[i-1] + fkMASKSTEP;
135   }
136   fInputStreams = new TClonesArray("AliStream",nInputStreams);
137   TClonesArray &lInputStreams = *fInputStreams;
138 // the first Input is open RW to be output as well
139   new(lInputStreams[0]) AliStream("UPDATE");
140   for (i=1;i<nInputStreams;i++) {
141     new(lInputStreams[i]) AliStream();
142   }
143   fOutput = 0;
144   fEvent = 0;
145   fNrOfEventsToWrite = -1;
146   fNrOfEventsWritten = 0;
147   fCopyTreesFromInput = -1;
148   fCombi = new AliMergeCombi(nInputStreams,sperb);
149   fDebug = 0;
150   if (GetDebug()>2) 
151     cerr<<"AliRunDigitizer::AliRunDigitizer() called"<<endl;
152   //PH
153   fTreeD = 0;
154   for (i=0; i<MAXSTREAMSTOMERGE; i++) fInputFiles[i]=0;
155 }
156
157 ////////////////////////////////////////////////////////////////////////
158
159 AliRunDigitizer::~AliRunDigitizer() {
160 // dtor
161
162   if (fInputStreams) {
163     delete fInputStreams;
164     fInputStreams = 0;
165   }
166   if (fCombi) {
167     delete fCombi;
168     fCombi = 0;
169   }
170
171 }
172 ////////////////////////////////////////////////////////////////////////
173 void AliRunDigitizer::AddDigitizer(AliDigitizer *digitizer)
174 {
175 // add digitizer to the list of active digitizers
176   this->Add(digitizer);
177 }
178 ////////////////////////////////////////////////////////////////////////
179 void AliRunDigitizer::SetInputStream(Int_t i, char *inputFile)
180 {
181   if (i > fInputStreams->GetLast()) {
182     Error("SetInputStream","Input stream number too high");
183     return;
184   }
185   static_cast<AliStream*>(fInputStreams->At(i))->AddFile(inputFile);
186 }
187
188 ////////////////////////////////////////////////////////////////////////
189 void AliRunDigitizer::Digitize(Option_t* option)
190 {
191 // get a new combination of inputs, connect input trees and loop 
192 // over all digitizers
193
194   if (!InitGlobal()) {
195     cerr<<"False from InitGlobal"<<endl;
196     return;
197   }
198 // take gAlice from the first input file. It is needed to access
199 //  geometry data
200   if (!static_cast<AliStream*>(fInputStreams->At(0))->ImportgAlice()) {
201     cerr<<"gAlice object not found in the first file of "
202         <<"the 1st stream"<<endl;
203     return;
204   }
205   Int_t eventsCreated = 0;
206 // loop until there is anything on the input in case fNrOfEventsToWrite < 0
207   while ((eventsCreated++ < fNrOfEventsToWrite) || (fNrOfEventsToWrite < 0)) {
208     if (!ConnectInputTrees()) break;
209     InitEvent();
210 // loop over all registered digitizers and let them do the work
211     ExecuteTasks("");
212     CleanTasks();
213     FinishEvent();
214   }
215   FinishGlobal();
216 }
217
218 ////////////////////////////////////////////////////////////////////////
219 Bool_t AliRunDigitizer::ConnectInputTrees()
220 {
221 // fill arrays fArrayTreeS, fArrayTreeH and fArrayTreeTPCS with 
222 // pointers to the correct events according fCombination values
223 // null pointers can be in the output, AliDigitizer has to check it
224
225   TTree *tree;
226   char treeName[50];
227   Int_t serialNr;
228   Int_t eventNr[MAXSTREAMSTOMERGE], delta[MAXSTREAMSTOMERGE];
229   fCombi->Combination(eventNr, delta);
230   for (Int_t i=0;i<fNinputs;i++) {
231     if (delta[i] == 1) {
232       AliStream *iStream = static_cast<AliStream*>(fInputStreams->At(i));
233       if (!iStream->NextEventInStream(serialNr)) return kFALSE;
234       fInputFiles[i]=iStream->CurrentFile();
235       sprintf(treeName,"TreeS%d",serialNr);
236       tree = static_cast<TTree*>(iStream->CurrentFile()->Get(treeName));
237       fArrayTreeS[i] = tree;
238       sprintf(treeName,"TreeH%d",serialNr);
239       tree = static_cast<TTree*>(iStream->CurrentFile()->Get(treeName));
240       fArrayTreeH[i] = tree;
241       sprintf(treeName,"TreeS_75x40_100x60_%d",serialNr);
242       tree = static_cast<TTree*>(iStream->CurrentFile()->Get(treeName));
243       fArrayTreeTPCS[i] = tree;
244       sprintf(treeName,"TreeS%d_TRD",serialNr);
245       tree = static_cast<TTree*>(iStream->CurrentFile()->Get(treeName));
246       fArrayTreeTRDS[i] = tree;
247     } else if (delta[i] != 0) {
248       Error("ConnectInputTrees","Only delta 0 or 1 is implemented");
249       return kFALSE;
250     }
251   }
252   return kTRUE;
253 }
254
255 ////////////////////////////////////////////////////////////////////////
256 Bool_t AliRunDigitizer::InitGlobal()
257 {
258 // called once before Digitize() is called, initialize digitizers and output
259
260   TList* subTasks = this->GetListOfTasks();
261   if (subTasks) {
262     subTasks->ForEach(AliDigitizer,Init)();
263   }  
264   return kTRUE;
265 }
266
267 ////////////////////////////////////////////////////////////////////////
268
269 void AliRunDigitizer::SetOutputFile(TString fn)
270 // the output will be to separate file, not to the signal file
271 {
272   fOutputFileName = fn;
273   (static_cast<AliStream*>(fInputStreams->At(0)))->ChangeMode("READ");
274   InitOutputGlobal();
275 }
276
277 ////////////////////////////////////////////////////////////////////////
278 Bool_t AliRunDigitizer::InitOutputGlobal()
279 {
280 // Creates the output file, called by InitEvent()
281
282   TString fn;
283   fn = fOutputDirName + '/' + fOutputFileName;
284   fOutput = new TFile(fn,"recreate");
285   if (GetDebug()>2) {
286     cerr<<"AliRunDigitizer::InitOutputGlobal(): file "<<fn.Data()<<" was opened"<<endl;
287   }
288   if (fOutput) return kTRUE;
289   Error("InitOutputGlobal","Could not create output file.");
290   return kFALSE;
291 }
292
293
294 ////////////////////////////////////////////////////////////////////////
295 void AliRunDigitizer::InitEvent()
296 {
297 // Creates TreeDxx in the output file, called from Digitize() once for 
298 //  each event. xx = fEvent
299
300   if (GetDebug()>2) 
301     cerr<<"AliRunDigitizer::InitEvent: fEvent = "<<fEvent<<endl;
302
303 // if fOutputFileName was not given, write output to signal file
304   if (fOutputFileName == "") {
305     fOutput = (static_cast<AliStream*>(fInputStreams->At(0)))->CurrentFile();
306   }
307   fOutput->cd();
308   char treeName[30];
309   sprintf(treeName,"TreeD%d",fEvent);
310   fTreeD = static_cast<TTree*>(fOutput->Get(treeName));
311   if (!fTreeD) {
312     fTreeD = new TTree(treeName,"Digits");
313     fTreeD->Write(0,TObject::kOverwrite);
314   }
315
316 // special tree for TPC
317   sprintf(treeName,"TreeD_75x40_100x60_%d",fEvent);
318   fTreeDTPC = static_cast<TTree*>(fOutput->Get(treeName));
319   if (!fTreeDTPC) {
320     fTreeDTPC = new TTree(treeName,"TPC_Digits");
321     fTreeDTPC->Write(0,TObject::kOverwrite);
322   }
323
324 // special tree for TRD
325   sprintf(treeName,"TreeD%d_TRD",fEvent);
326   fTreeDTRD = static_cast<TTree*>(fOutput->Get(treeName));
327   if (!fTreeDTRD) {
328     fTreeDTRD = new TTree(treeName,"TRD_Digits");
329     fTreeDTRD->Write(0,TObject::kOverwrite);
330   }
331
332 }
333
334 ////////////////////////////////////////////////////////////////////////
335 void AliRunDigitizer::FinishEvent()
336 {
337 // called at the end of loop over digitizers
338
339   fOutput->cd();
340   if (fCopyTreesFromInput > -1) {
341     char treeName[20];
342     Int_t i = fCopyTreesFromInput; 
343     sprintf(treeName,"TreeK%d",fCombination[i]);
344     fInputFiles[i]->Get(treeName)->Clone()->Write();
345     sprintf(treeName,"TreeH%d",fCombination[i]);
346     fInputFiles[i]->Get(treeName)->Clone()->Write();
347   }
348   fEvent++;
349   fNrOfEventsWritten++;
350   if (fTreeD) {
351     delete fTreeD;
352     fTreeD = 0;
353   }
354   if (fTreeDTPC) {
355     delete fTreeDTPC;
356     fTreeDTPC = 0;
357   }
358   if (fTreeDTRD) {
359     delete fTreeDTRD;
360     fTreeDTRD = 0;
361   }
362 }
363 ////////////////////////////////////////////////////////////////////////
364 void AliRunDigitizer::FinishGlobal()
365 {
366 // called at the end of Exec
367 // save unique objects to the output file
368
369   fOutput->cd();
370   this->Write();
371   if (fCopyTreesFromInput > -1) {
372     fInputFiles[fCopyTreesFromInput]->Get("TE")->Clone()->Write();
373     gAlice->Write();
374   }
375   fOutput->Close();
376 }
377
378
379 ////////////////////////////////////////////////////////////////////////
380 Int_t  AliRunDigitizer::GetNParticles(Int_t event)
381 {
382 // return number of particles in all input files for a given
383 // event (as numbered in the output file)
384 // return -1 if some file cannot be accessed
385
386   Int_t sum = 0;
387   Int_t sumI;
388   for (Int_t i = 0; i < fNinputs; i++) {
389     sumI = GetNParticles(GetInputEventNumber(event,i), i);
390     if (sumI < 0) return -1;
391     sum += sumI;
392   }
393   return sum;
394 }
395
396 ////////////////////////////////////////////////////////////////////////
397 Int_t  AliRunDigitizer::GetNParticles(Int_t event, Int_t input)
398 {
399 // return number of particles in input file input for a given
400 // event (as numbered in this input file)
401 // return -1 if some error
402
403 // Must be revised in the version with AliStream
404
405   return -1;
406
407 /*
408   TFile *file = ConnectInputFile(input);
409   if (!file) {
410     Error("GetNParticles","Cannot open input file");
411     return -1;
412   }
413
414 // find the header and get Nprimaries and Nsecondaries
415   TTree* tE = (TTree *)file->Get("TE") ;
416   if (!tE) {
417     Error("GetNParticles","input file does not contain TE");
418     return -1;
419   }
420   AliHeader* header;
421   header = 0;
422   tE->SetBranchAddress("Header", &header);
423   if (!tE->GetEntry(event)) {
424     Error("GetNParticles","event %d not found",event);
425     return -1;
426   }
427   if (GetDebug()>2) {
428     cerr<<"Nprimary: "<< header->GetNprimary()<<endl;
429     cerr<<"Nsecondary: "<<header->GetNsecondary()<<endl;
430   }
431   return header->GetNprimary() + header->GetNsecondary();
432 */
433 }
434
435 ////////////////////////////////////////////////////////////////////////
436 Int_t* AliRunDigitizer::GetInputEventNumbers(Int_t event)
437 {
438 // return pointer to an int array with input event numbers which were
439 // merged in the output event event
440
441 // simplified for now, implement later
442   Int_t * a = new Int_t[MAXSTREAMSTOMERGE];
443   for (Int_t i = 0; i < fNinputs; i++) {
444     a[i] = event;
445   }
446   return a;
447 }
448 ////////////////////////////////////////////////////////////////////////
449 Int_t AliRunDigitizer::GetInputEventNumber(Int_t event, Int_t input)
450 {
451 // return an event number of an eventInput from input file input
452 // which was merged to create output event event
453
454 // simplified for now, implement later
455   return event;
456 }
457 ////////////////////////////////////////////////////////////////////////
458 TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t event)
459 {
460 // return pointer to particle with index i (index with mask)
461
462 // decode the MASK
463   Int_t input = i/fkMASKSTEP;
464   return GetParticle(i,input,GetInputEventNumber(event,input));
465 }
466
467 ////////////////////////////////////////////////////////////////////////
468 TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t input, Int_t event)
469 {
470 // return pointer to particle with index i in the input file input
471 // (index without mask)
472 // event is the event number in the file input
473 // return 0 i fit does not exist
474
475 // Must be revised in the version with AliStream
476
477   return 0;
478 /*
479   TFile *file = ConnectInputFile(input);
480   if (!file) {
481     Error("GetParticle","Cannot open input file");
482     return 0;
483   }
484
485 // find the header and get Nprimaries and Nsecondaries
486   TTree* tE = (TTree *)file->Get("TE") ;
487   if (!tE) {
488     Error("GetParticle","input file does not contain TE");
489     return 0;
490   }
491   AliHeader* header;
492   header = 0;
493   tE->SetBranchAddress("Header", &header);
494   if (!tE->GetEntry(event)) {
495     Error("GetParticle","event %d not found",event);
496     return 0;
497   }
498   
499 // connect TreeK  
500   char treeName[30];
501   sprintf(treeName,"TreeK%d",event);  
502   TTree* tK = static_cast<TTree*>(file->Get(treeName));
503   if (!tK) {
504     Error("GetParticle","input file does not contain TreeK%d",event);
505     return 0;
506   }
507   TParticle *particleBuffer;
508   particleBuffer = 0;
509   tK->SetBranchAddress("Particles", &particleBuffer);
510
511
512 // algorithmic way of getting entry index
513 // (primary particles are filled after secondaries)
514   Int_t entry;
515   if (i<header->GetNprimary())
516     entry = i+header->GetNsecondary();
517   else 
518     entry = i-header->GetNprimary();
519   Int_t bytesRead = tK->GetEntry(entry);
520 //  new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
521   if (bytesRead)
522     return particleBuffer;
523   return  0;
524 */
525 }
526
527 ////////////////////////////////////////////////////////////////////////
528 void AliRunDigitizer::ExecuteTask(Option_t* option)
529 {
530 // overwrite ExecuteTask to do Digitize only
531
532   if (!IsActive()) return;
533   Digitize(option);
534   fHasExecuted = kTRUE;
535   return;
536 }