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