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