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