correct bug in the initialization
[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.2  2001/07/28 10:44:32  hristov
19 Loop variable declared once; typos corrected
20
21 Revision 1.1  2001/07/27 12:59:00  jchudoba
22 Manager class for merging/digitization
23
24 */
25
26 ////////////////////////////////////////////////////////////////////////
27 //
28 // AliRunDigitizer.cxx
29 //
30 // Manager object for merging/digitization
31 //
32 // Instance of this class manages the digitization and/or merging of
33 // Sdigits into Digits. 
34 //
35 // Only one instance of this class is created in the macro:
36 //   AliRunDigitizer * manager = new AliRunDigitizer();
37 // Then instances of specific detector digitizers are created:
38 //   AliMUONDigitizer *dMUON  = new AliMUONDigitizer(manager)
39 // and the manager is configured (you have to specify input files 
40 // and an output file). The manager generates a combination of input 
41 // events according to an option set by SetCombinationType. Then it
42 // connects appropriate trees from the input files, creates TreeD 
43 // in the output and runs once per event Digitize method of all existing 
44 // AliDetDigitizers (without any option). AliDetDigitizers ask manager
45 // for a TTree with input (manager->GetNextTreeH(TTree *) 
46 // or manager->GetNextTreeS(TTree *),
47 // merge all inputs, digitize it, and save it in the TreeD 
48 // obtained by manager->GetTreeD(). Output events are stored with 
49 // numbers from 0, this default can be changed by 
50 // manager->SetFirstOutputEventNr(Int_t) method. The particle numbers
51 // in the output are shifted by MASK, which is taken from manager.
52 //
53 // Single input file is permitted. Maximum MAXFILESTOMERGE can be merged.
54 // Input from the memory (on-the-fly merging) is not yet 
55 // supported, as well as access to the input data by invoking methods
56 // on the output data.
57 //
58 // Access to the geometrical data is via gAlice for now (supposing the 
59 // same geometry in all input files), gAlice is taken from the defined 
60 // input file.
61 //
62 // Example with MUON digitizer:
63 //
64 //   AliRunDigitizer * manager = new AliRunDigitizer();
65 //   manager->SetInput("1track_10events_phi45_60.root");
66 //   manager->SetInput("1track_10events_phi120_135.root");
67 //   manager->SetOutputDir("/home/z2/jchudoba/batch/jobtmp");
68 //   manager->SetOutputFile("digits.root");
69 //   AliMUONDigitizer *dMUON  = new AliMUONDigitizer(manager);
70 //   manager->SetNrOfEventsToWrite(3);
71 //   manager->SetCopyTreesFromInput(1);
72 //   manager->Digitize();
73 //
74 //////////////////////////////////////////////////////////////////////// 
75
76 // system includes
77
78 #include <iostream.h>
79
80 // ROOT includes
81
82 #include "TFile.h"
83 #include "TTree.h"
84
85 // AliROOT includes
86
87 #include "AliRunDigitizer.h"
88 #include "AliDigitizer.h"
89 #include "AliRun.h"
90 #include "AliHeader.h"
91 #include "TParticle.h"
92
93 ClassImp(AliRunDigitizer)
94
95 ////////////////////////////////////////////////////////////////////////
96
97   AliRunDigitizer::AliRunDigitizer() : TNamed("AliRunDigitizer","")
98 {
99 // default ctor
100
101   Int_t i;
102
103   for (i=0;i<MAXDETECTORS;i++) fDigitizers[i]=0;
104   fNDigitizers = 0;
105   fNinputs = 0;
106   fOutputFileName = "digits.root";
107   fOutputDirName = "/tmp/";
108   fCombination.Set(MAXFILESTOMERGE);
109   for (i=0;i<MAXFILESTOMERGE;i++) {
110     fArrayTreeS[i]=fArrayTreeH[i]=fArrayTreeTPCS[i]=NULL;
111     fCombination[i]=-1;
112   }
113   fkMASKSTEP = 10000000;
114   fkMASK[0] = 0;
115   for (i=1;i<MAXFILESTOMERGE;i++) {
116     fkMASK[i] = fkMASK[i-1] + fkMASKSTEP;
117   }
118   fInputFileNames = new TClonesArray("TObjString",1);
119   fInputFiles = new TClonesArray("TFile",1);
120   fMinNEvents = 99999999;
121   fCombinationType=1;
122   fCombinationFileName = "combination.txt";
123   fOutput = 0;
124   fEvent = 0;
125   fNrOfEventsToWrite = 0;
126   fNrOfEventsWritten = 0;
127   fCopyTreesFromInput = -1;
128   fDebug = 3;
129   if (GetDebug()>2) 
130     cerr<<"AliRunDigitizer::AliRunDigitizer() called"<<endl;
131 }
132
133 ////////////////////////////////////////////////////////////////////////
134
135 AliRunDigitizer::~AliRunDigitizer() {
136 // dtor
137
138   if (fInputFiles) delete fInputFiles;
139   if (fInputFileNames) delete fInputFileNames;
140 }
141
142 ////////////////////////////////////////////////////////////////////////
143 void AliRunDigitizer::AddDigitizer(AliDigitizer *digitizer)
144 {
145 // add digitizer to the list of active digitizers
146
147   if (fNDigitizers >= MAXDETECTORS) {
148     cerr<<"Too many detectors to digitize. Increase value of MAXDETECTORS"
149         <<" constant in AliRunDigitizer.h and recompile or decrease the"
150         <<" the number of detectors"<<endl;
151   } else {
152     fDigitizers[fNDigitizers++] = digitizer;
153   }
154 }
155
156 ////////////////////////////////////////////////////////////////////////
157
158 Bool_t AliRunDigitizer::SetInput(char *inputFile)
159 {
160 // receives the name of the input file. Opens it and stores pointer 
161 // to it, returns kFALSE if fails
162 // if inputFile = MEMORY, uses pointers from gAlice - not yet implemented
163
164 // input cannot be output - open READ only
165   TFile *file = new((*fInputFiles)[fNinputs]) TFile(inputFile,"READ"); 
166
167   if (GetDebug()>2) cerr<<"AliRunDigitizer::SetInput: file = "<<file<<endl;
168   if (!file->IsOpen()) {
169     cerr<<"ERROR: AliRunDigitizer::SetInput: cannot open file "
170         <<inputFile<<endl;
171     return kFALSE;
172   }
173
174 // find Header and get nr of events there
175   TTree * te = (TTree *) file->Get("TE") ;
176   if (!te) {
177     cerr<<"ERROR: AliRunDigitizer::SetInput:  input file does "
178         <<"not contain TE"<<endl;
179     return kFALSE;
180   }
181   Int_t nEntries = (Int_t) te->GetEntries();
182
183   if (GetDebug()>2) cerr<<"AliRunDigitizer::SetInput: nEntries = "
184                         <<nEntries<<endl;
185   if (nEntries < 1) {
186     cerr<<"ERROR: AliRunDigitizer::SetInput:  input file does "
187         <<"not contain any event"<<endl;
188     return kFALSE;
189   }
190
191   if (nEntries < fMinNEvents) fNrOfEventsToWrite = fMinNEvents = nEntries;
192
193 // find gAlice object if it is a first input file
194 //  this is only temporary solution, we need gAlice to access detector
195 //  geometry parameters. Unfortunately I have to include AliRun header file.
196   if (fNinputs == 0) {
197     gAlice = (AliRun*)file->Get("gAlice");
198     if (GetDebug() > 2) cerr<<"gAlice taken from the first input: "
199                             <<gAlice<<endl;
200     if (!gAlice) {
201       cerr<<"ERROR: AliRunDigitizer::SetInput:  first input file "
202           <<"does not contain gAlice object"<<endl;
203       return kFALSE;
204     }
205   }
206     
207 // store this file name if it is OK
208   new((*fInputFileNames)[fNinputs])  TObjString(inputFile);
209   fNinputs++;
210   if (GetDebug() > 2) cerr<<"fNinputs = "<<fNinputs<<endl;
211   return kTRUE;
212 }
213
214 ////////////////////////////////////////////////////////////////////////
215 Bool_t AliRunDigitizer::MakeCombination()
216 {
217 // make a new combination of events from different files
218
219   Int_t type = fCombinationType;
220   if (fNrOfEventsWritten >= fNrOfEventsToWrite) return kFALSE;
221
222   switch (type) {
223         
224   case 1:
225 // type = 1:  1-1-1 - take the same event number from each file
226     if (fCombination[0]<fMinNEvents-1) {        
227       for (Int_t i=0;i<fNinputs;i++) fCombination[i]++;
228       return kTRUE;
229     }
230     if (GetDebug()>2)
231       cerr<<"AliRunDigitizer::MakeCombination: Warning: "
232           <<"maximum number of Events in an input file "
233           <<"was reached."<<endl;
234     break;
235
236   case 2:
237 // type = 2: read from the file combinations.ascii
238 //  not yet implemented 100% correctly - requires 4 entries in the row
239     static FILE *fp ;
240     static Int_t linesRead;
241     if (!fp) {
242       fp = fopen(fCombinationFileName.Data(),"r");
243       linesRead = 0;
244     }
245     if (!fp) {
246       cerr<<"AliRunDigitizer::MakeCombination ERROR: "
247           <<"Cannot open input file with combinations."<<endl;
248       return kFALSE;
249     }
250     Int_t var[4], nInputs;
251     char line[80];
252 // since I do not close or rewind the file, the position should be correct
253     if (fgets(line,80,fp)) {
254       nInputs = sscanf(&line[0],"%d%d%d%d",&var[0],&var[1],&var[2],&var[3]);
255       if (nInputs != fNinputs) {
256         cerr<<"AliRunDigitizer::MakeCombination ERROR: "
257             <<"Nr. of input files is different from nr "
258             <<"integers defining the combination"<<endl;
259         return kFALSE;
260       }
261       while(nInputs--) {
262         fCombination[nInputs] = var[nInputs];
263       }
264       return kTRUE;
265     } else {
266       cerr<<"AliRunDigitizer::MakeCombination ERROR: "
267           <<"no more input in the file with combinations"<<endl;
268       return kFALSE;
269     }
270
271   default:
272     cerr<<"AliRunDigitizer::MakeCombination: ERROR: "
273         <<"wrong type of required combination type: "<<type<<endl;
274   }
275   return kFALSE;
276 }
277
278 ////////////////////////////////////////////////////////////////////////
279 void AliRunDigitizer::PrintCombination()
280 {
281 // debug method to print current combination
282
283   cerr<<"AliRunDigitizer::PrintCombination: Current events combination:"<<endl;
284   for (Int_t i=0;i<fNinputs;i++) {
285     cerr<<"File: "<<((TObjString *)fInputFileNames->At(i))->GetString()<<"\tEvent: "<<fCombination[i]<<endl;
286   }
287 }
288
289 ////////////////////////////////////////////////////////////////////////
290 void AliRunDigitizer::Digitize()
291 {
292 // get a new combination of inputs, connect input trees and loop 
293 // over all digitizers
294
295   if (!InitGlobal()) {
296     cerr<<"False from InitGlobal"<<endl;
297     return;
298   }
299 // take gAlice from the first input file. It is needed to access
300 //  geometry data
301   while (MakeCombination()) {
302     if (GetDebug()>2) PrintCombination();
303     ConnectInputTrees();
304     InitPerEvent();
305 // loop over all registered digitizers and let them do the work
306     for (Int_t i=0;i<fNDigitizers; i++) {
307       fDigitizers[i]->Digitize();
308     }
309     FinishPerEvent();
310   }
311   FinishGlobal();
312 }
313
314 ////////////////////////////////////////////////////////////////////////
315 void AliRunDigitizer::ConnectInputTrees()
316 {
317 // fill arrays fArrayTreeS, fArrayTreeH and fArrayTreeTPCS with 
318 // pointers to the correct events according fCombination values
319 // null pointers can be in the output, AliDigitizer has to check it
320
321   TFile *file;
322   TTree *tree;
323   char treeName[20];
324   for (Int_t i=0;i<fNinputs;i++) {
325     file = (TFile*)(*fInputFiles)[i];
326     sprintf(treeName,"TreeS%d",fCombination[i]);
327     tree = (TTree *) file->Get(treeName);
328     fArrayTreeS[i] = tree;
329     sprintf(treeName,"TreeH%d",fCombination[i]);
330     tree = (TTree *) file->Get(treeName);
331     fArrayTreeH[i] = tree;
332     sprintf(treeName,"TreeD_75x40_100x60_%d",fCombination[i]);
333     tree = (TTree *) file->Get(treeName);
334     fArrayTreeTPCS[i] = tree;
335    }
336 }
337
338 ////////////////////////////////////////////////////////////////////////
339 Bool_t AliRunDigitizer::InitGlobal()
340 {
341 // called once before Digitize() is called, initialize digitizers and output
342
343   if (!InitOutputGlobal()) return kFALSE;
344   for (Int_t i=0;i<fNDigitizers; i++) {
345     if (!fDigitizers[i]->Init()) return kFALSE;
346   }
347   return kTRUE;
348 }
349
350 ////////////////////////////////////////////////////////////////////////
351 Bool_t AliRunDigitizer::InitOutputGlobal()
352 {
353 // Creates the output file, called once by InitGlobal()
354
355   TString fn;
356   fn = fOutputDirName + '/' + fOutputFileName;
357   fOutput = new TFile(fn,"recreate");
358   if (GetDebug()>2) {
359     cerr<<"file "<<fn.Data()<<" was opened"<<endl;
360     cerr<<"fOutput = "<<fOutput<<endl;
361   }
362   if (fOutput) return kTRUE;
363   Error("InitOutputGlobal","Could not create output file.");
364   return kFALSE;
365 }
366
367
368 ////////////////////////////////////////////////////////////////////////
369 void AliRunDigitizer::InitPerEvent()
370 {
371 // Creates TreeDxx in the output file, called from Digitize() once for 
372 //  each event. xx = fEvent
373
374   if (GetDebug()>2) 
375     cerr<<"AliRunDigitizer::InitPerEvent: fEvent = "<<fEvent<<endl;
376   fOutput->cd();
377   char hname[30];
378   sprintf(hname,"TreeD%d",fEvent);
379   fTreeD = new TTree(hname,"Digits");
380 //  fTreeD->Write();   // Do I have to write it here???
381 }
382
383
384 ////////////////////////////////////////////////////////////////////////
385 TTree* AliRunDigitizer::GetNextTreeH(TTree *current) const
386 {
387 // returns next one after the current one
388 // returns the first if the current is NULL
389 // returns NULL if the current is the last one
390
391   if (fNinputs <= 0) return 0;
392   if (current == 0) return fArrayTreeH[0];
393   for (Int_t i=0; i<fNinputs-1; i++) {
394     if (current == fArrayTreeH[i]) return fArrayTreeH[i+1];
395   }
396   return 0;
397 }
398 ////////////////////////////////////////////////////////////////////////
399 TTree* AliRunDigitizer::GetNextTreeS(TTree *current) const
400 {
401 // returns next one after the current one
402 // returns the first if the current is NULL
403 // returns NULL if the current is the last one
404
405   if (fNinputs <= 0) return 0;
406   if (current == 0) return fArrayTreeS[0];
407   for (Int_t i=0; i<fNinputs-1; i++) {
408     if (current == fArrayTreeS[i]) return fArrayTreeS[i+1];
409   }
410   return 0;
411 }
412 ////////////////////////////////////////////////////////////////////////
413 TTree* AliRunDigitizer::GetNextTreeTPCS(TTree *current) const
414 {
415 // returns next one after the current one
416 // returns the first if the current is NULL
417 // returns NULL if the current is the last one
418
419   if (fNinputs <= 0) return 0;
420   if (current == 0) return fArrayTreeTPCS[0];
421   for (Int_t i=0; i<fNinputs-1; i++) {
422     if (current == fArrayTreeTPCS[i]) return fArrayTreeTPCS[i+1];
423   }
424   return 0;
425 }
426
427 ////////////////////////////////////////////////////////////////////////
428 Int_t AliRunDigitizer::GetNextMask(Int_t current) const
429 {
430 // returns next one after the current one
431 // returns the first if the current is negative
432 // returns negative if the current is the last one
433
434   if (fNinputs <= 0) return -1;
435   if (current < 0) return fkMASK[0]; 
436   for (Int_t i=0; i<fNinputs-1; i++) {
437     if (current == fkMASK[i]) return fkMASK[i+1];
438   }
439   return -1;
440 }
441 ////////////////////////////////////////////////////////////////////////
442 void AliRunDigitizer::FinishPerEvent()
443 {
444 // called at the end of loop over digitizers
445
446   fOutput->cd();
447   if (fCopyTreesFromInput > -1) {
448     char treeName[20];
449     Int_t i = fCopyTreesFromInput; 
450     sprintf(treeName,"TreeK%d",fCombination[i]);
451     ((TFile*)fInputFiles->At(i))->Get(treeName)->Clone()->Write();
452     sprintf(treeName,"TreeH%d",fCombination[i]);
453     ((TFile*)fInputFiles->At(i))->Get(treeName)->Clone()->Write();
454   }
455   fEvent++;
456   fNrOfEventsWritten++;
457   if (fTreeD) {
458     delete fTreeD;
459     fTreeD = 0;
460   }
461 }
462 ////////////////////////////////////////////////////////////////////////
463 void AliRunDigitizer::FinishGlobal()
464 {
465 // called at the end of Exec
466 // save unique objects to the output file
467
468   fOutput->cd();
469   this->Write();
470   if (fCopyTreesFromInput > -1) {
471     ((TFile*)fInputFiles->At(fCopyTreesFromInput))->Get("TE")
472       ->Clone()->Write();
473     gAlice->Write();
474   }
475   fOutput->Close();
476 }
477
478
479 ////////////////////////////////////////////////////////////////////////
480 Int_t  AliRunDigitizer::GetNParticles(Int_t event)
481 {
482 // return number of particles in all input files for a given
483 // event (as numbered in the output file)
484 // return -1 if some file cannot be accessed
485
486   Int_t sum = 0;
487   Int_t sumI;
488   for (Int_t i = 0; i < fNinputs; i++) {
489     sumI = GetNParticles(GetInputEventNumber(event,i), i);
490     if (sumI < 0) return -1;
491     sum += sumI;
492   }
493   return sum;
494 }
495
496 ////////////////////////////////////////////////////////////////////////
497 Int_t  AliRunDigitizer::GetNParticles(Int_t event, Int_t input)
498 {
499 // return number of particles in input file input for a given
500 // event (as numbered in this input file)
501 // return -1 if some error
502
503   TFile *file = ConnectInputFile(input);
504   if (!file) {
505     Error("GetNParticles","Cannot open input file");
506     return -1;
507   }
508
509 // find the header and get Nprimaries and Nsecondaries
510   TTree* tE = (TTree *)file->Get("TE") ;
511   if (!tE) {
512     Error("GetNParticles","input file does not contain TE");
513     return -1;
514   }
515   AliHeader* header;
516   header = 0;
517   tE->SetBranchAddress("Header", &header);
518   if (!tE->GetEntry(event)) {
519     Error("GetNParticles","event %d not found",event);
520     return -1;
521   }
522   if (GetDebug()>2) {
523     cerr<<"Nprimary: "<< header->GetNprimary()<<endl;
524     cerr<<"Nsecondary: "<<header->GetNsecondary()<<endl;
525   }
526   return header->GetNprimary() + header->GetNsecondary();
527 }
528
529 ////////////////////////////////////////////////////////////////////////
530 Int_t* AliRunDigitizer::GetInputEventNumbers(Int_t event)
531 {
532 // return pointer to an int array with input event numbers which were
533 // merged in the output event event
534
535 // simplified for now, implement later
536   Int_t a[MAXFILESTOMERGE];
537   for (Int_t i = 0; i < fNinputs; i++) {
538     a[i] = event;
539   }
540   return a;
541 }
542 ////////////////////////////////////////////////////////////////////////
543 Int_t AliRunDigitizer::GetInputEventNumber(Int_t event, Int_t input)
544 {
545 // return an event number of an eventInput from input file input
546 // which was merged to create output event event
547
548 // simplified for now, implement later
549   return event;
550 }
551 ////////////////////////////////////////////////////////////////////////
552 TFile* AliRunDigitizer::ConnectInputFile(Int_t input)
553 {
554 // open input file with index input
555 // 1st file has index 0
556 // return 0x0 if fails
557
558 // check if file with index input is already open
559   TFile *file = 0;
560   if (fInputFiles->GetEntriesFast() > input)
561     file = static_cast<TFile *>(fInputFiles->At(input));
562
563   if (!file) {
564     // find the file name and open it
565     TObjString *fn = static_cast<TObjString *>(fInputFileNames->At(input));
566     file = new((*fInputFiles)[input]) TFile((fn->GetString()).Data(),"READ");
567     if (!file) {
568       Error("ConnectInputFile","Cannot open input file");
569       return 0;
570     }
571     if (!file->IsOpen()) {
572       Error("ConnectInputFile","Cannot open input file");
573       return 0;
574     }
575   }
576   return file;
577 }
578
579 ////////////////////////////////////////////////////////////////////////
580 TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t event)
581 {
582 // return pointer to particle with index i (index with mask)
583
584 // decode the MASK
585   Int_t input = i/fkMASKSTEP;
586   return GetParticle(i,input,GetInputEventNumber(event,input));
587 }
588
589 ////////////////////////////////////////////////////////////////////////
590 TParticle* AliRunDigitizer::GetParticle(Int_t i, Int_t input, Int_t event)
591 {
592 // return pointer to particle with index i in the input file input
593 // (index without mask)
594 // event is the event number in the file input
595 // return 0 i fit does not exist
596
597   TFile *file = ConnectInputFile(input);
598   if (!file) {
599     Error("GetParticle","Cannot open input file");
600     return 0;
601   }
602
603 // find the header and get Nprimaries and Nsecondaries
604   TTree* tE = (TTree *)file->Get("TE") ;
605   if (!tE) {
606     Error("GetParticle","input file does not contain TE");
607     return 0;
608   }
609   AliHeader* header;
610   header = 0;
611   tE->SetBranchAddress("Header", &header);
612   if (!tE->GetEntry(event)) {
613     Error("GetParticle","event %d not found",event);
614     return 0;
615   }
616   
617 // connect TreeK  
618   char treeName[30];
619   sprintf(treeName,"TreeK%d",event);  
620   TTree* tK = static_cast<TTree*>(file->Get(treeName));
621   if (!tK) {
622     Error("GetParticle","input file does not contain TreeK%d",event);
623     return 0;
624   }
625   TParticle *particleBuffer;
626   particleBuffer = 0;
627   tK->SetBranchAddress("Particles", &particleBuffer);
628
629
630 // algorithmic way of getting entry index
631 // (primary particles are filled after secondaries)
632   Int_t entry;
633   if (i<header->GetNprimary())
634     entry = i+header->GetNsecondary();
635   else 
636     entry = i-header->GetNprimary();
637 //  tK->GetEntry(0);          // do I need this???
638   Int_t bytesRead = tK->GetEntry(entry);
639 //  new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
640   if (bytesRead)
641     return particleBuffer;
642   return  0;
643 }
644 ////////////////////////////////////////////////////////////////////////
645