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