New code to record the timing information of various methods in simulation and recons...
[u/mrichter/AliRoot.git] / STEER / AliSimulation.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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // class for running generation, simulation and digitization                 //
21 //                                                                           //
22 // Hits, sdigits and digits are created for all detectors by typing:         //
23 //                                                                           //
24 //   AliSimulation sim;                                                      //
25 //   sim.Run();                                                              //
26 //                                                                           //
27 // The Run method returns kTRUE in case of successful execution.             //
28 // The number of events can be given as argument to the Run method or it     //
29 // can be set by                                                             //
30 //                                                                           //
31 //   sim.SetNumberOfEvents(n);                                               //
32 //                                                                           //
33 // The name of the configuration file can be passed as argument to the       //
34 // AliSimulation constructor or can be specified by                          //
35 //                                                                           //
36 //   sim.SetConfigFile("...");                                               //
37 //                                                                           //
38 // The generation of particles and the simulation of detector hits can be    //
39 // switched on or off by                                                     //
40 //                                                                           //
41 //   sim.SetRunGeneration(kTRUE);   // generation of primary particles       //
42 //   sim.SetRunSimulation(kFALSE);  // but no tracking                       //
43 //                                                                           //
44 // For which detectors sdigits and digits will be created, can be steered    //
45 // by                                                                        //
46 //                                                                           //
47 //   sim.SetMakeSDigits("ALL");     // make sdigits for all detectors        //
48 //   sim.SetMakeDigits("ITS TPC");  // make digits only for ITS and TPC      //
49 //                                                                           //
50 // The argument is a (case sensitive) string with the names of the           //
51 // detectors separated by a space. An empty string ("") can be used to       //
52 // disable the creation of sdigits or digits. The special string "ALL"       //
53 // selects all available detectors. This is the default.                     //
54 //                                                                           //
55 // The creation of digits from hits instead of from sdigits can be selected  //
56 // by                                                                        //
57 //                                                                           //
58 //   sim.SetMakeDigitsFromHits("TRD");                                       //
59 //                                                                           //
60 // The argument is again a string with the selected detectors. Be aware that //
61 // this feature is not available for all detectors and that merging is not   //
62 // possible, when digits are created directly from hits.                     //
63 //                                                                           //
64 // Background events can be merged by calling                                //
65 //                                                                           //
66 //   sim.MergeWith("background/galice.root", 2);                             //
67 //                                                                           //
68 // The first argument is the file name of the background galice file. The    //
69 // second argument is the number of signal events per background event.      //
70 // By default this number is calculated from the number of available         //
71 // background events. MergeWith can be called several times to merge more    //
72 // than two event streams. It is assumed that the sdigits were already       //
73 // produced for the background events.                                       //
74 //                                                                           //
75 // The output of raw data can be switched on by calling                      //
76 //                                                                           //
77 //   sim.SetWriteRawData("MUON");   // write raw data for MUON               //
78 //                                                                           //
79 // The default output format of the raw data are DDL files. They are         //
80 // converted to a DATE file, if a file name is given as second argument.     //
81 // For this conversion the program "dateStream" is required. If the file     //
82 // name has the extension ".root", the DATE file is converted to a root      //
83 // file. The program "alimdc" is used for this purpose. For the conversion   //
84 // to DATE and root format the two conversion programs have to be installed. //
85 // Only the raw data in the final format is kept if the third argument is    //
86 // kTRUE.                                                                    //
87 //                                                                           //
88 // The methods RunSimulation, RunSDigitization, RunDigitization,             //
89 // RunHitsDigitization and WriteRawData can be used to run only parts of     //
90 // the full simulation chain. The creation of raw data DDL files and their   //
91 // conversion to the DATE or root format can be run directly by calling      //
92 // the methods WriteRawFiles, ConvertRawFilesToDate and ConvertDateToRoot.   //
93 //                                                                           //
94 // The default number of events per file, which is usually set in the        //
95 // config file, can be changed for individual detectors and data types       //
96 // by calling                                                                //
97 //                                                                           //
98 //   sim.SetEventsPerFile("PHOS", "Reconstructed Points", 3);                //
99 //                                                                           //
100 // The first argument is the detector, the second one the data type and the  //
101 // last one the number of events per file. Valid data types are "Hits",      //
102 // "Summable Digits", "Digits", "Reconstructed Points" and "Tracks".         //
103 // The number of events per file has to be set before the simulation of      //
104 // hits. Otherwise it has no effect.                                         //
105 //                                                                           //
106 ///////////////////////////////////////////////////////////////////////////////
107
108 #include <TVirtualMCApplication.h>
109 #include <TGeoManager.h>
110 #include <TObjString.h>
111 #include <TSystem.h>
112 #include <TFile.h>
113
114 #include "AliCodeTimer.h"
115 #include "AliCDBStorage.h"
116 #include "AliCDBEntry.h"
117 #include "AliCDBManager.h"
118 #include "AliGeomManager.h"
119 #include "AliAlignObj.h"
120 #include "AliCentralTrigger.h"
121 #include "AliDAQ.h"
122 #include "AliDigitizer.h"
123 #include "AliGenerator.h"
124 #include "AliLog.h"
125 #include "AliModule.h"
126 #include "AliRun.h"
127 #include "AliRunDigitizer.h"
128 #include "AliRunLoader.h"
129 #include "AliSimulation.h"
130 #include "AliVertexGenFile.h"
131 #include "AliCentralTrigger.h"
132 #include "AliCTPRawData.h"
133 #include "AliRawReaderFile.h"
134 #include "AliESD.h"
135 #include "AliHeader.h"
136 #include "AliGenEventHeader.h"
137 #include "AliMC.h"
138
139 ClassImp(AliSimulation)
140
141 AliSimulation *AliSimulation::fgInstance = 0;
142
143 //_____________________________________________________________________________
144 AliSimulation::AliSimulation(const char* configFileName, const char* cdbUri,
145                              const char* name, const char* title) :
146   TNamed(name, title),
147
148   fRunGeneration(kTRUE),
149   fRunSimulation(kTRUE),
150   fLoadAlignFromCDB(kTRUE),
151   fLoadAlObjsListOfDets("ALL"),
152   fMakeSDigits("ALL"),
153   fMakeDigits("ALL"),
154   fMakeTrigger(""),
155   fMakeDigitsFromHits(""),
156   fWriteRawData(""),
157   fRawDataFileName(""),
158   fDeleteIntermediateFiles(kFALSE),
159   fStopOnError(kFALSE),
160
161   fNEvents(1),
162   fConfigFileName(configFileName),
163   fGAliceFileName("galice.root"),
164   fEventsPerFile(),
165   fBkgrdFileNames(NULL),
166   fAlignObjArray(NULL),
167   fUseBkgrdVertex(kTRUE),
168   fRegionOfInterest(kFALSE),
169   fCDBUri(cdbUri),
170   fSpecCDBUri(),
171   fEmbeddingFlag(kFALSE)
172 {
173 // create simulation object with default parameters
174   fgInstance = this;
175   SetGAliceFile("galice.root");
176 }
177
178 //_____________________________________________________________________________
179 AliSimulation::AliSimulation(const AliSimulation& sim) :
180   TNamed(sim),
181
182   fRunGeneration(sim.fRunGeneration),
183   fRunSimulation(sim.fRunSimulation),
184   fLoadAlignFromCDB(sim.fLoadAlignFromCDB),
185   fLoadAlObjsListOfDets(sim.fLoadAlObjsListOfDets),
186   fMakeSDigits(sim.fMakeSDigits),
187   fMakeDigits(sim.fMakeDigits),
188   fMakeTrigger(sim.fMakeTrigger),
189   fMakeDigitsFromHits(sim.fMakeDigitsFromHits),
190   fWriteRawData(sim.fWriteRawData),
191   fRawDataFileName(""),
192   fDeleteIntermediateFiles(kFALSE),
193   fStopOnError(sim.fStopOnError),
194
195   fNEvents(sim.fNEvents),
196   fConfigFileName(sim.fConfigFileName),
197   fGAliceFileName(sim.fGAliceFileName),
198   fEventsPerFile(),
199   fBkgrdFileNames(NULL),
200   fAlignObjArray(NULL),
201   fUseBkgrdVertex(sim.fUseBkgrdVertex),
202   fRegionOfInterest(sim.fRegionOfInterest),
203   fCDBUri(sim.fCDBUri),
204   fSpecCDBUri(),
205   fEmbeddingFlag(sim.fEmbeddingFlag)
206 {
207 // copy constructor
208
209   for (Int_t i = 0; i < sim.fEventsPerFile.GetEntriesFast(); i++) {
210     if (!sim.fEventsPerFile[i]) continue;
211     fEventsPerFile.Add(sim.fEventsPerFile[i]->Clone());
212   }
213
214   fBkgrdFileNames = new TObjArray;
215   for (Int_t i = 0; i < sim.fBkgrdFileNames->GetEntriesFast(); i++) {
216     if (!sim.fBkgrdFileNames->At(i)) continue;
217     fBkgrdFileNames->Add(sim.fBkgrdFileNames->At(i)->Clone());
218   }
219
220   for (Int_t i = 0; i < sim.fSpecCDBUri.GetEntriesFast(); i++) {
221     if (sim.fSpecCDBUri[i]) fSpecCDBUri.Add(sim.fSpecCDBUri[i]->Clone());
222   }
223   fgInstance = this;
224 }
225
226 //_____________________________________________________________________________
227 AliSimulation& AliSimulation::operator = (const AliSimulation& sim)
228 {
229 // assignment operator
230
231   this->~AliSimulation();
232   new(this) AliSimulation(sim);
233   return *this;
234 }
235
236 //_____________________________________________________________________________
237 AliSimulation::~AliSimulation()
238 {
239 // clean up
240
241   fEventsPerFile.Delete();
242 //  if(fAlignObjArray) fAlignObjArray->Delete(); // fAlignObjArray->RemoveAll() ???
243 //  delete fAlignObjArray; fAlignObjArray=0;
244
245   if (fBkgrdFileNames) {
246     fBkgrdFileNames->Delete();
247     delete fBkgrdFileNames;
248   }
249
250   fSpecCDBUri.Delete();
251   if (fgInstance==this) fgInstance = 0;
252
253   AliCodeTimer::Instance()->Print();
254 }
255
256
257 //_____________________________________________________________________________
258 void AliSimulation::SetNumberOfEvents(Int_t nEvents)
259 {
260 // set the number of events for one run
261
262   fNEvents = nEvents;
263 }
264
265 //_____________________________________________________________________________
266 void AliSimulation::InitCDBStorage()
267 {
268 // activate a default CDB storage
269 // First check if we have any CDB storage set, because it is used 
270 // to retrieve the calibration and alignment constants
271
272   AliCDBManager* man = AliCDBManager::Instance();
273   if (man->IsDefaultStorageSet())
274   {
275     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
276     AliWarning("Default CDB storage has been already set !");
277     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
278     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
279     fCDBUri = "";
280   }
281   else {
282     AliDebug(2,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
283     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
284     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
285     man->SetDefaultStorage(fCDBUri);
286   }
287
288   // Now activate the detector specific CDB storage locations
289   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
290     TObject* obj = fSpecCDBUri[i];
291     if (!obj) continue;
292     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
293     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
294     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
295     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
296   }
297   man->Print();
298 }
299
300 //_____________________________________________________________________________
301 void AliSimulation::SetDefaultStorage(const char* uri) {
302 // Store the desired default CDB storage location
303 // Activate it later within the Run() method
304
305   fCDBUri = uri;
306
307 }
308
309 //_____________________________________________________________________________
310 void AliSimulation::SetSpecificStorage(const char* calibType, const char* uri) {
311 // Store a detector-specific CDB storage location
312 // Activate it later within the Run() method
313
314   AliCDBPath aPath(calibType);
315   if(!aPath.IsValid()){
316         AliError(Form("Not a valid path: %s", calibType));
317         return;
318   }
319
320   TObject* obj = fSpecCDBUri.FindObject(calibType);
321   if (obj) fSpecCDBUri.Remove(obj);
322   fSpecCDBUri.Add(new TNamed(calibType, uri));
323
324 }
325
326 //_____________________________________________________________________________
327 void AliSimulation::SetConfigFile(const char* fileName)
328 {
329 // set the name of the config file
330
331   fConfigFileName = fileName;
332 }
333
334 //_____________________________________________________________________________
335 void AliSimulation::SetGAliceFile(const char* fileName)
336 {
337 // set the name of the galice file
338 // the path is converted to an absolute one if it is relative
339
340   fGAliceFileName = fileName;
341   if (!gSystem->IsAbsoluteFileName(fGAliceFileName)) {
342     char* absFileName = gSystem->ConcatFileName(gSystem->WorkingDirectory(),
343                                                 fGAliceFileName);
344     fGAliceFileName = absFileName;
345     delete[] absFileName;
346   }
347
348   AliDebug(2, Form("galice file name set to %s", fileName));
349 }
350
351 //_____________________________________________________________________________
352 void AliSimulation::SetEventsPerFile(const char* detector, const char* type, 
353                                      Int_t nEvents)
354 {
355 // set the number of events per file for the given detector and data type
356 // ("Hits", "Summable Digits", "Digits", "Reconstructed Points" or "Tracks")
357
358   TNamed* obj = new TNamed(detector, type);
359   obj->SetUniqueID(nEvents);
360   fEventsPerFile.Add(obj);
361 }
362
363 //_____________________________________________________________________________
364 Bool_t AliSimulation::MisalignGeometry(AliRunLoader *runLoader)
365 {
366   // Read the alignment objects from CDB.
367   // Each detector is supposed to have the
368   // alignment objects in DET/Align/Data CDB path.
369   // All the detector objects are then collected,
370   // sorted by geometry level (starting from ALIC) and
371   // then applied to the TGeo geometry.
372   // Finally an overlaps check is performed.
373
374   if (!AliGeomManager::GetGeometry() || !AliGeomManager::GetGeometry()->IsClosed()) {
375     AliError("Can't apply the misalignment! Geometry is not loaded or it is still opened!");
376     return kFALSE;
377   }  
378   Bool_t delRunLoader = kFALSE;
379   if (!runLoader) {
380     runLoader = LoadRun("READ");
381     if (!runLoader) return kFALSE;
382     delRunLoader = kTRUE;
383   }
384
385   // Export ideal geometry 
386   AliGeomManager::GetGeometry()->Export("geometry.root");
387
388   // Load alignment data from CDB and apply to geometry through AliGeomManager
389   if(fLoadAlignFromCDB){
390     
391     TString detStr = fLoadAlObjsListOfDets;
392     TString loadAlObjsListOfDets = "";
393     
394     TObjArray* detArray = runLoader->GetAliRun()->Detectors();
395     for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
396       AliModule* det = (AliModule*) detArray->At(iDet);
397       if (!det || !det->IsActive()) continue;
398       if (IsSelected(det->GetName(), detStr)) {
399         //add det to list of dets to be aligned from CDB
400         loadAlObjsListOfDets += det->GetName();
401         loadAlObjsListOfDets += " ";
402       }
403     } // end loop over detectors
404     AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
405   }else{
406     // Check if the array with alignment objects was
407     // provided by the user. If yes, apply the objects
408     // to the present TGeo geometry
409     if (fAlignObjArray) {
410       if (AliGeomManager::ApplyAlignObjsToGeom(*fAlignObjArray) == kFALSE) {
411         AliError("The misalignment of one or more volumes failed!"
412                  "Compare the list of simulated detectors and the list of detector alignment data!");
413         if (delRunLoader) delete runLoader;
414         return kFALSE;
415       }
416     }
417   }
418
419   // Update the internal geometry of modules (ITS needs it)
420   TString detStr = fLoadAlObjsListOfDets;
421   TObjArray* detArray = runLoader->GetAliRun()->Detectors();
422   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
423
424     AliModule* det = (AliModule*) detArray->At(iDet);
425     if (!det || !det->IsActive()) continue;
426     if (IsSelected(det->GetName(), detStr)) {
427       det->UpdateInternalGeometry();
428     }
429   } // end loop over detectors
430
431
432   if (delRunLoader) delete runLoader;
433
434   return kTRUE;
435 }
436
437
438 //_____________________________________________________________________________
439 Bool_t AliSimulation::SetRunNumber()
440 {
441   // Set the CDB manager run number
442   // The run number is retrieved from gAlice
443
444   if(AliCDBManager::Instance()->GetRun() < 0) {
445     AliRunLoader* runLoader = LoadRun("READ");
446     if (!runLoader) return kFALSE;
447     else {
448       AliCDBManager::Instance()->SetRun(runLoader->GetAliRun()->GetRunNumber());
449       AliInfo(Form("Run number: %d",AliCDBManager::Instance()->GetRun()));
450       delete runLoader;
451     }
452   }
453   return kTRUE;
454 }
455
456 //_____________________________________________________________________________
457 void AliSimulation::MergeWith(const char* fileName, Int_t nSignalPerBkgrd)
458 {
459 // add a file with background events for merging
460
461   TObjString* fileNameStr = new TObjString(fileName);
462   fileNameStr->SetUniqueID(nSignalPerBkgrd);
463   if (!fBkgrdFileNames) fBkgrdFileNames = new TObjArray;
464   fBkgrdFileNames->Add(fileNameStr);
465 }
466
467 void AliSimulation::EmbedInto(const char* fileName, Int_t nSignalPerBkgrd)
468 {
469 // add a file with background events for embeddin
470   MergeWith(fileName, nSignalPerBkgrd);
471   fEmbeddingFlag = kTRUE;
472 }
473
474 //_____________________________________________________________________________
475 Bool_t AliSimulation::Run(Int_t nEvents)
476 {
477 // run the generation, simulation and digitization
478
479   AliCodeTimerAuto("")
480   
481   InitCDBStorage();
482
483   if (nEvents > 0) fNEvents = nEvents;
484
485   // generation and simulation -> hits
486   if (fRunGeneration) {
487     if (!RunSimulation()) if (fStopOnError) return kFALSE;
488   }
489
490   // Set run number in CDBManager (if it is not already set in RunSimulation)
491   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
492
493   // If RunSimulation was not called, load the geometry and misalign it
494   if (!AliGeomManager::GetGeometry()) {
495     // Initialize the geometry manager
496     AliGeomManager::LoadGeometry("geometry.root");
497     if (!AliGeomManager::GetGeometry()) if (fStopOnError) return kFALSE;
498     // Misalign geometry
499     if(!MisalignGeometry()) if (fStopOnError) return kFALSE;
500   }
501
502   // hits -> summable digits
503   if (!fMakeSDigits.IsNull()) {
504     if (!RunSDigitization(fMakeSDigits)) if (fStopOnError) return kFALSE;
505   }
506
507   // summable digits -> digits
508   if (!fMakeDigits.IsNull()) {
509     if (!RunDigitization(fMakeDigits, fMakeDigitsFromHits)) {
510       if (fStopOnError) return kFALSE;
511     }
512   }
513
514   // hits -> digits
515   if (!fMakeDigitsFromHits.IsNull()) {
516     if (fBkgrdFileNames && (fBkgrdFileNames->GetEntriesFast() > 0)) {
517       AliWarning(Form("Merging and direct creation of digits from hits " 
518                  "was selected for some detectors. "
519                  "No merging will be done for the following detectors: %s",
520                  fMakeDigitsFromHits.Data()));
521     }
522     if (!RunHitsDigitization(fMakeDigitsFromHits)) {
523       if (fStopOnError) return kFALSE;
524     }
525   }
526
527   // digits -> trigger
528   if (!RunTrigger(fMakeTrigger)) {
529     if (fStopOnError) return kFALSE;
530   }
531
532   // digits -> raw data
533   if (!fWriteRawData.IsNull()) {
534     if (!WriteRawData(fWriteRawData, fRawDataFileName, 
535                       fDeleteIntermediateFiles)) {
536       if (fStopOnError) return kFALSE;
537     }
538   }
539
540   return kTRUE;
541 }
542
543 //_____________________________________________________________________________
544 Bool_t AliSimulation::RunTrigger(const char* descriptors)
545 {
546   // run the trigger
547
548   AliCodeTimerAuto("")
549
550    AliRunLoader* runLoader = LoadRun("READ");
551    if (!runLoader) return kFALSE;
552    TString des = descriptors;
553
554    if (des.IsNull()) {
555      if (gAlice->GetTriggerDescriptor() != "") {
556        des = gAlice->GetTriggerDescriptor();
557      }
558      else {
559        AliWarning("No trigger descriptor is specified. Skipping the trigger simulation...");
560        return kTRUE;
561      }
562    }
563
564    runLoader->MakeTree( "GG" );
565    AliCentralTrigger* aCTP = runLoader->GetTrigger();
566   // Load Descriptors
567    aCTP->LoadDescriptor( des );
568
569   // digits -> trigger
570    if( !aCTP->RunTrigger( runLoader ) ) {
571       if (fStopOnError) {
572     //  delete aCTP;
573          return kFALSE;
574       }
575    }
576
577    delete runLoader;
578
579    return kTRUE;
580 }
581
582 //_____________________________________________________________________________
583 Bool_t AliSimulation::WriteTriggerRawData()
584 {
585   // Writes the CTP (trigger) DDL raw data
586   // Details of the format are given in the
587   // trigger TDR - pages 134 and 135.
588   AliCTPRawData writer;
589   writer.RawData();
590
591   return kTRUE;
592 }
593
594 //_____________________________________________________________________________
595 Bool_t AliSimulation::RunSimulation(Int_t nEvents)
596 {
597 // run the generation and simulation
598
599   AliCodeTimerAuto("")
600
601   if (!gAlice) {
602     AliError("no gAlice object. Restart aliroot and try again.");
603     return kFALSE;
604   }
605   if (gAlice->Modules()->GetEntries() > 0) {
606     AliError("gAlice was already run. Restart aliroot and try again.");
607     return kFALSE;
608   }
609
610   AliInfo(Form("initializing gAlice with config file %s",
611           fConfigFileName.Data()));
612   StdoutToAliInfo(StderrToAliError(
613     gAlice->Init(fConfigFileName.Data());
614   ););
615
616   // Get the trigger descriptor string
617   // Either from AliSimulation or from
618   // gAlice
619   if (fMakeTrigger.IsNull()) {
620     if (gAlice->GetTriggerDescriptor() != "")
621       fMakeTrigger = gAlice->GetTriggerDescriptor();
622   }
623   else
624     gAlice->SetTriggerDescriptor(fMakeTrigger.Data());
625
626   // Set run number in CDBManager
627   AliInfo(Form("Run number: %d",AliCDBManager::Instance()->GetRun()));
628
629   AliRunLoader* runLoader = gAlice->GetRunLoader();
630   if (!runLoader) {
631              AliError(Form("gAlice has no run loader object. "
632                              "Check your config file: %s", fConfigFileName.Data()));
633              return kFALSE;
634   }
635   SetGAliceFile(runLoader->GetFileName());
636  
637   // Misalign geometry
638 #if ROOT_VERSION_CODE < 331527
639   AliGeomManager::SetGeometry(gGeoManager);
640   MisalignGeometry(runLoader);
641 #endif
642
643 //   AliRunLoader* runLoader = gAlice->GetRunLoader();
644 //   if (!runLoader) {
645 //     AliError(Form("gAlice has no run loader object. "
646 //                   "Check your config file: %s", fConfigFileName.Data()));
647 //     return kFALSE;
648 //   }
649 //   SetGAliceFile(runLoader->GetFileName());
650
651   if (!gAlice->Generator()) {
652     AliError(Form("gAlice has no generator object. "
653                   "Check your config file: %s", fConfigFileName.Data()));
654     return kFALSE;
655   }
656   if (nEvents <= 0) nEvents = fNEvents;
657
658   // get vertex from background file in case of merging
659   if (fUseBkgrdVertex &&
660       fBkgrdFileNames && (fBkgrdFileNames->GetEntriesFast() > 0)) {
661     Int_t signalPerBkgrd = GetNSignalPerBkgrd(nEvents);
662     const char* fileName = ((TObjString*)
663                             (fBkgrdFileNames->At(0)))->GetName();
664     AliInfo(Form("The vertex will be taken from the background "
665                  "file %s with nSignalPerBackground = %d", 
666                  fileName, signalPerBkgrd));
667     AliVertexGenFile* vtxGen = new AliVertexGenFile(fileName, signalPerBkgrd);
668     gAlice->Generator()->SetVertexGenerator(vtxGen);
669   }
670
671   if (!fRunSimulation) {
672     gAlice->Generator()->SetTrackingFlag(0);
673   }
674
675   // set the number of events per file for given detectors and data types
676   for (Int_t i = 0; i < fEventsPerFile.GetEntriesFast(); i++) {
677     if (!fEventsPerFile[i]) continue;
678     const char* detName = fEventsPerFile[i]->GetName();
679     const char* typeName = fEventsPerFile[i]->GetTitle();
680     TString loaderName(detName);
681     loaderName += "Loader";
682     AliLoader* loader = runLoader->GetLoader(loaderName);
683     if (!loader) {
684       AliError(Form("RunSimulation", "no loader for %s found\n"
685                     "Number of events per file not set for %s %s", 
686                     detName, typeName, detName));
687       continue;
688     }
689     AliDataLoader* dataLoader = 
690       loader->GetDataLoader(typeName);
691     if (!dataLoader) {
692       AliError(Form("no data loader for %s found\n"
693                     "Number of events per file not set for %s %s", 
694                     typeName, detName, typeName));
695       continue;
696     }
697     dataLoader->SetNumberOfEventsPerFile(fEventsPerFile[i]->GetUniqueID());
698     AliDebug(1, Form("number of events per file set to %d for %s %s",
699                      fEventsPerFile[i]->GetUniqueID(), detName, typeName));
700   }
701
702   AliInfo("running gAlice");
703   StdoutToAliInfo(StderrToAliError(
704     gAlice->Run(nEvents);
705   ););
706
707   delete runLoader;
708
709
710   return kTRUE;
711 }
712
713 //_____________________________________________________________________________
714 Bool_t AliSimulation::RunSDigitization(const char* detectors)
715 {
716 // run the digitization and produce summable digits
717
718   AliCodeTimerAuto("")
719
720   AliRunLoader* runLoader = LoadRun();
721   if (!runLoader) return kFALSE;
722
723   TString detStr = detectors;
724   TObjArray* detArray = runLoader->GetAliRun()->Detectors();
725   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
726     AliModule* det = (AliModule*) detArray->At(iDet);
727     if (!det || !det->IsActive()) continue;
728     if (IsSelected(det->GetName(), detStr)) {
729       AliInfo(Form("creating summable digits for %s", det->GetName()));
730       AliCodeTimerAuto(Form("creating summable digits for %s", det->GetName()));
731       det->Hits2SDigits();
732     }
733   }
734
735   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
736     AliError(Form("the following detectors were not found: %s",
737                   detStr.Data()));
738     if (fStopOnError) return kFALSE;
739   }
740
741   delete runLoader;
742
743   return kTRUE;
744 }
745
746
747 //_____________________________________________________________________________
748 Bool_t AliSimulation::RunDigitization(const char* detectors, 
749                                       const char* excludeDetectors)
750 {
751 // run the digitization and produce digits from sdigits
752
753   AliCodeTimerAuto("")
754
755   while (AliRunLoader::GetRunLoader()) delete AliRunLoader::GetRunLoader();
756   if (gAlice) delete gAlice;
757   gAlice = NULL;
758
759   Int_t nStreams = 1;
760   if (fBkgrdFileNames) nStreams = fBkgrdFileNames->GetEntriesFast() + 1;
761   Int_t signalPerBkgrd = GetNSignalPerBkgrd();
762   AliRunDigitizer* manager = new AliRunDigitizer(nStreams, signalPerBkgrd);
763   // manager->SetEmbeddingFlag(fEmbeddingFlag);
764   manager->SetInputStream(0, fGAliceFileName.Data());
765   for (Int_t iStream = 1; iStream < nStreams; iStream++) {
766     const char* fileName = ((TObjString*)
767                             (fBkgrdFileNames->At(iStream-1)))->GetName();
768     manager->SetInputStream(iStream, fileName);
769   }
770
771   TString detStr = detectors;
772   TString detExcl = excludeDetectors;
773   manager->GetInputStream(0)->ImportgAlice();
774   AliRunLoader* runLoader = 
775     AliRunLoader::GetRunLoader(manager->GetInputStream(0)->GetFolderName());
776   TObjArray* detArray = runLoader->GetAliRun()->Detectors();
777   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
778     AliModule* det = (AliModule*) detArray->At(iDet);
779     if (!det || !det->IsActive()) continue;
780     if (IsSelected(det->GetName(), detStr) && 
781         !IsSelected(det->GetName(), detExcl)) {
782       AliDigitizer* digitizer = det->CreateDigitizer(manager);
783       
784       if (!digitizer) {
785         AliError(Form("no digitizer for %s", det->GetName()));
786         if (fStopOnError) return kFALSE;
787       } else {
788         digitizer->SetRegionOfInterest(fRegionOfInterest);
789       }
790     }
791   }
792
793   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
794     AliError(Form("the following detectors were not found: %s", 
795                   detStr.Data()));
796     if (fStopOnError) return kFALSE;
797   }
798
799   if (!manager->GetListOfTasks()->IsEmpty()) {
800     AliInfo("executing digitization");
801     manager->Exec("");
802   }
803
804   delete manager;
805
806   return kTRUE;
807 }
808
809 //_____________________________________________________________________________
810 Bool_t AliSimulation::RunHitsDigitization(const char* detectors)
811 {
812 // run the digitization and produce digits from hits
813
814   AliCodeTimerAuto("")
815
816   AliRunLoader* runLoader = LoadRun("READ");
817   if (!runLoader) return kFALSE;
818
819   TString detStr = detectors;
820   TObjArray* detArray = runLoader->GetAliRun()->Detectors();
821   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
822     AliModule* det = (AliModule*) detArray->At(iDet);
823     if (!det || !det->IsActive()) continue;
824     if (IsSelected(det->GetName(), detStr)) {
825       AliInfo(Form("creating digits from hits for %s", det->GetName()));
826       det->Hits2Digits();
827     }
828   }
829
830   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
831     AliError(Form("the following detectors were not found: %s", 
832                   detStr.Data()));
833     if (fStopOnError) return kFALSE;
834   }
835
836   delete runLoader;
837   //PH Temporary fix to avoid interference with the PHOS loder/getter
838   //PH The problem has to be solved in more general way 09/06/05
839
840   return kTRUE;
841 }
842
843 //_____________________________________________________________________________
844 Bool_t AliSimulation::WriteRawData(const char* detectors, 
845                                    const char* fileName,
846                                    Bool_t deleteIntermediateFiles)
847 {
848 // convert the digits to raw data
849 // First DDL raw data files for the given detectors are created.
850 // If a file name is given, the DDL files are then converted to a DATE file.
851 // If deleteIntermediateFiles is true, the DDL raw files are deleted 
852 // afterwards.
853 // If the file name has the extension ".root", the DATE file is converted
854 // to a root file.
855 // If deleteIntermediateFiles is true, the DATE file is deleted afterwards.
856
857   AliCodeTimerAuto("")
858
859   if (!WriteRawFiles(detectors)) {
860     if (fStopOnError) return kFALSE;
861   }
862
863   TString dateFileName(fileName);
864   if (!dateFileName.IsNull()) {
865     Bool_t rootOutput = dateFileName.EndsWith(".root");
866     if (rootOutput) dateFileName += ".date";
867     if (!ConvertRawFilesToDate(dateFileName)) {
868       if (fStopOnError) return kFALSE;
869     }
870     if (deleteIntermediateFiles) {
871       AliRunLoader* runLoader = LoadRun("READ");
872       if (runLoader) for (Int_t iEvent = 0; 
873                           iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
874         char command[256];
875         sprintf(command, "rm -r raw%d", iEvent);
876         gSystem->Exec(command);
877       }
878     }
879
880     if (rootOutput) {
881       if (!ConvertDateToRoot(dateFileName, fileName)) {
882         if (fStopOnError) return kFALSE;
883       }
884       if (deleteIntermediateFiles) {
885         gSystem->Unlink(dateFileName);
886       }
887     }
888   }
889
890   return kTRUE;
891 }
892
893 //_____________________________________________________________________________
894 Bool_t AliSimulation::WriteRawFiles(const char* detectors)
895 {
896 // convert the digits to raw data DDL files
897
898   AliCodeTimerAuto("")
899   
900   AliRunLoader* runLoader = LoadRun("READ");
901   if (!runLoader) return kFALSE;
902
903   // write raw data to DDL files
904   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
905     AliInfo(Form("processing event %d", iEvent));
906     runLoader->GetEvent(iEvent);
907     TString baseDir = gSystem->WorkingDirectory();
908     char dirName[256];
909     sprintf(dirName, "raw%d", iEvent);
910     gSystem->MakeDirectory(dirName);
911     if (!gSystem->ChangeDirectory(dirName)) {
912       AliError(Form("couldn't change to directory %s", dirName));
913       if (fStopOnError) return kFALSE; else continue;
914     }
915
916     TString detStr = detectors;
917     TObjArray* detArray = runLoader->GetAliRun()->Detectors();
918     for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
919       AliModule* det = (AliModule*) detArray->At(iDet);
920       if (!det || !det->IsActive()) continue;
921       if (IsSelected(det->GetName(), detStr)) {
922         AliInfo(Form("creating raw data from digits for %s", det->GetName()));
923         det->Digits2Raw();
924       }
925     }
926
927     if (!WriteTriggerRawData())
928       if (fStopOnError) return kFALSE;
929
930     gSystem->ChangeDirectory(baseDir);
931     if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
932       AliError(Form("the following detectors were not found: %s", 
933                     detStr.Data()));
934       if (fStopOnError) return kFALSE;
935     }
936   }
937
938   delete runLoader;
939   
940   return kTRUE;
941 }
942
943 //_____________________________________________________________________________
944 Bool_t AliSimulation::ConvertRawFilesToDate(const char* dateFileName)
945 {
946 // convert raw data DDL files to a DATE file with the program "dateStream"
947
948   AliCodeTimerAuto("")
949   
950   char* path = gSystem->Which(gSystem->Getenv("PATH"), "dateStream");
951   if (!path) {
952     AliError("the program dateStream was not found");
953     if (fStopOnError) return kFALSE;
954   } else {
955     delete[] path;
956   }
957
958   AliRunLoader* runLoader = LoadRun("READ");
959   if (!runLoader) return kFALSE;
960
961   AliInfo(Form("converting raw data DDL files to DATE file %s", dateFileName));
962   char command[256];
963   // Note the option -s. It is used in order to avoid
964   // the generation of SOR/EOR events.
965   sprintf(command, "dateStream -s -D -o %s -# %d -C", 
966           dateFileName, runLoader->GetNumberOfEvents());
967   FILE* pipe = gSystem->OpenPipe(command, "w");
968
969   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
970     fprintf(pipe, "GDC\n");
971     Float_t ldc = 0;
972     Int_t prevLDC = -1;
973
974     // loop over detectors and DDLs
975     for (Int_t iDet = 0; iDet < AliDAQ::kNDetectors; iDet++) {
976       for (Int_t iDDL = 0; iDDL < AliDAQ::NumberOfDdls(iDet); iDDL++) {
977
978         Int_t ddlID = AliDAQ::DdlID(iDet,iDDL);
979         Int_t ldcID = Int_t(ldc + 0.0001);
980         ldc += AliDAQ::NumberOfLdcs(iDet) / AliDAQ::NumberOfDdls(iDet);
981
982         char rawFileName[256];
983         sprintf(rawFileName, "raw%d/%s", 
984                 iEvent, AliDAQ::DdlFileName(iDet,iDDL));
985
986         // check existence and size of raw data file
987         FILE* file = fopen(rawFileName, "rb");
988         if (!file) continue;
989         fseek(file, 0, SEEK_END);
990         unsigned long size = ftell(file);
991         fclose(file);
992         if (!size) continue;
993
994         if (ldcID != prevLDC) {
995           fprintf(pipe, " LDC Id %d\n", ldcID);
996           prevLDC = ldcID;
997         }
998         fprintf(pipe, "  Equipment Id %d Payload %s\n", ddlID, rawFileName);
999       }
1000     }
1001   }
1002
1003   Int_t result = gSystem->ClosePipe(pipe);
1004
1005   delete runLoader;
1006   return (result == 0);
1007 }
1008
1009 //_____________________________________________________________________________
1010 Bool_t AliSimulation::ConvertDateToRoot(const char* dateFileName,
1011                                         const char* rootFileName)
1012 {
1013 // convert a DATE file to a root file with the program "alimdc"
1014
1015   // ALIMDC setup
1016   const Int_t kDBSize = 2000000000;
1017   const Int_t kTagDBSize = 1000000000;
1018   const Bool_t kFilter = kFALSE;
1019   const Int_t kCompression = 0;
1020
1021   char* path = gSystem->Which(gSystem->Getenv("PATH"), "alimdc");
1022   if (!path) {
1023     AliError("the program alimdc was not found");
1024     if (fStopOnError) return kFALSE;
1025   } else {
1026     delete[] path;
1027   }
1028
1029   AliInfo(Form("converting DATE file %s to root file %s", 
1030                dateFileName, rootFileName));
1031
1032   const char* rawDBFS[2] = { "/tmp/mdc1", "/tmp/mdc2" };
1033   const char* tagDBFS    = "/tmp/mdc1/tags";
1034
1035   // User defined file system locations
1036   if (gSystem->Getenv("ALIMDC_RAWDB1")) 
1037     rawDBFS[0] = gSystem->Getenv("ALIMDC_RAWDB1");
1038   if (gSystem->Getenv("ALIMDC_RAWDB2")) 
1039     rawDBFS[1] = gSystem->Getenv("ALIMDC_RAWDB2");
1040   if (gSystem->Getenv("ALIMDC_TAGDB")) 
1041     tagDBFS = gSystem->Getenv("ALIMDC_TAGDB");
1042
1043   gSystem->Exec(Form("rm -rf %s",rawDBFS[0]));
1044   gSystem->Exec(Form("rm -rf %s",rawDBFS[1]));
1045   gSystem->Exec(Form("rm -rf %s",tagDBFS));
1046
1047   gSystem->Exec(Form("mkdir %s",rawDBFS[0]));
1048   gSystem->Exec(Form("mkdir %s",rawDBFS[1]));
1049   gSystem->Exec(Form("mkdir %s",tagDBFS));
1050
1051   Int_t result = gSystem->Exec(Form("alimdc %d %d %d %d %s", 
1052                                     kDBSize, kTagDBSize, kFilter, kCompression, dateFileName));
1053   gSystem->Exec(Form("mv %s/*.root %s", rawDBFS[0], rootFileName));
1054
1055   gSystem->Exec(Form("rm -rf %s",rawDBFS[0]));
1056   gSystem->Exec(Form("rm -rf %s",rawDBFS[1]));
1057   gSystem->Exec(Form("rm -rf %s",tagDBFS));
1058
1059   return (result == 0);
1060 }
1061
1062
1063 //_____________________________________________________________________________
1064 AliRunLoader* AliSimulation::LoadRun(const char* mode) const
1065 {
1066 // delete existing run loaders, open a new one and load gAlice
1067
1068   while (AliRunLoader::GetRunLoader()) delete AliRunLoader::GetRunLoader();
1069   AliRunLoader* runLoader = 
1070     AliRunLoader::Open(fGAliceFileName.Data(), 
1071                        AliConfig::GetDefaultEventFolderName(), mode);
1072   if (!runLoader) {
1073     AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1074     return NULL;
1075   }
1076   runLoader->LoadgAlice();
1077   gAlice = runLoader->GetAliRun();
1078   if (!gAlice) {
1079     AliError(Form("no gAlice object found in file %s", 
1080                   fGAliceFileName.Data()));
1081     return NULL;
1082   }
1083   return runLoader;
1084 }
1085
1086 //_____________________________________________________________________________
1087 Int_t AliSimulation::GetNSignalPerBkgrd(Int_t nEvents) const
1088 {
1089 // get or calculate the number of signal events per background event
1090
1091   if (!fBkgrdFileNames) return 1;
1092   Int_t nBkgrdFiles = fBkgrdFileNames->GetEntriesFast();
1093   if (nBkgrdFiles == 0) return 1;
1094
1095   // get the number of signal events
1096   if (nEvents <= 0) {
1097     AliRunLoader* runLoader = 
1098         AliRunLoader::Open(fGAliceFileName.Data(), "SIGNAL");
1099     if (!runLoader) return 1;
1100     
1101     nEvents = runLoader->GetNumberOfEvents();
1102     delete runLoader;
1103   }
1104
1105   Int_t result = 0;
1106   for (Int_t iBkgrdFile = 0; iBkgrdFile < nBkgrdFiles; iBkgrdFile++) {
1107     // get the number of background events
1108     const char* fileName = ((TObjString*)
1109                             (fBkgrdFileNames->At(iBkgrdFile)))->GetName();
1110     AliRunLoader* runLoader =
1111       AliRunLoader::Open(fileName, "BKGRD");
1112     if (!runLoader) continue;
1113     Int_t nBkgrdEvents = runLoader->GetNumberOfEvents();
1114     delete runLoader;
1115   
1116     // get or calculate the number of signal per background events
1117     Int_t nSignalPerBkgrd = fBkgrdFileNames->At(iBkgrdFile)->GetUniqueID();
1118     if (nSignalPerBkgrd <= 0) {
1119       nSignalPerBkgrd = (nEvents-1) / nBkgrdEvents + 1;
1120     } else if (result && (result != nSignalPerBkgrd)) {
1121       AliInfo(Form("the number of signal events per background event "
1122                    "will be changed from %d to %d for stream %d", 
1123                    nSignalPerBkgrd, result, iBkgrdFile+1));
1124       nSignalPerBkgrd = result;
1125     }
1126
1127     if (!result) result = nSignalPerBkgrd;
1128     if (nSignalPerBkgrd * nBkgrdEvents < nEvents) {
1129       AliWarning(Form("not enough background events (%d) for %d signal events "
1130                       "using %d signal per background events for stream %d",
1131                       nBkgrdEvents, nEvents, nSignalPerBkgrd, iBkgrdFile+1));
1132     }
1133   }
1134
1135   return result;
1136 }
1137
1138 //_____________________________________________________________________________
1139 Bool_t AliSimulation::IsSelected(TString detName, TString& detectors) const
1140 {
1141 // check whether detName is contained in detectors
1142 // if yes, it is removed from detectors
1143
1144   // check if all detectors are selected
1145   if ((detectors.CompareTo("ALL") == 0) ||
1146       detectors.BeginsWith("ALL ") ||
1147       detectors.EndsWith(" ALL") ||
1148       detectors.Contains(" ALL ")) {
1149     detectors = "ALL";
1150     return kTRUE;
1151   }
1152
1153   // search for the given detector
1154   Bool_t result = kFALSE;
1155   if ((detectors.CompareTo(detName) == 0) ||
1156       detectors.BeginsWith(detName+" ") ||
1157       detectors.EndsWith(" "+detName) ||
1158       detectors.Contains(" "+detName+" ")) {
1159     detectors.ReplaceAll(detName, "");
1160     result = kTRUE;
1161   }
1162
1163   // clean up the detectors string
1164   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1165   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1166   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1167
1168   return result;
1169 }
1170
1171 Bool_t AliSimulation::ConvertRaw2SDigits(const char* rawDirectory, const char* esdFileName) 
1172 {
1173 //
1174 // Steering routine  to convert raw data in directory rawDirectory/ to fake SDigits. 
1175 // These can be used for embedding of MC tracks into RAW data using the standard 
1176 // merging procedure.
1177 //
1178 // If an ESD file is given the reconstructed vertex is taken from it and stored in the event header.
1179 //
1180     if (!gAlice) {
1181         AliError("no gAlice object. Restart aliroot and try again.");
1182         return kFALSE;
1183     }
1184     if (gAlice->Modules()->GetEntries() > 0) {
1185         AliError("gAlice was already run. Restart aliroot and try again.");
1186         return kFALSE;
1187     }
1188     
1189     AliInfo(Form("initializing gAlice with config file %s",fConfigFileName.Data()));
1190     StdoutToAliInfo(StderrToAliError(gAlice->Init(fConfigFileName.Data());););
1191 //
1192 //  Initialize CDB     
1193     InitCDBStorage();
1194     AliCDBManager* man = AliCDBManager::Instance();
1195     man->SetRun(0); // Should this come from rawdata header ?
1196     
1197     Int_t iDet;
1198     //
1199     // Get the runloader
1200     AliRunLoader* runLoader = gAlice->GetRunLoader();
1201     //
1202     // Open esd file if available
1203     TFile* esdFile = TFile::Open(esdFileName);
1204     Bool_t esdOK = (esdFile != 0);
1205     AliESD* esd = new AliESD;
1206     TTree* treeESD = 0;
1207     if (esdOK) {
1208         treeESD = (TTree*) esdFile->Get("esdTree");
1209         if (!treeESD) {
1210             AliWarning("No ESD tree found");
1211             esdOK = kFALSE;
1212         } else {
1213             treeESD->SetBranchAddress("ESD", &esd);
1214         }
1215     }
1216     //
1217     // Create the RawReader
1218     AliRawReaderFile* rawReader = new AliRawReaderFile(rawDirectory);
1219     //
1220     // Get list of detectors
1221     TObjArray* detArray = runLoader->GetAliRun()->Detectors();
1222     //
1223     // Get Header
1224     AliHeader* header = runLoader->GetHeader();
1225     //
1226     // Event loop
1227     Int_t nev = 0;
1228     while(kTRUE) {
1229         if (!(rawReader->NextEvent())) break;
1230         //
1231         // Detector loop
1232         for (iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
1233             AliModule* det = (AliModule*) detArray->At(iDet);
1234             AliInfo(Form("Calling Raw2SDigits for %s\n", det->GetName()));
1235             det->Raw2SDigits(rawReader);
1236             rawReader->Reset();
1237         } // detectors
1238
1239         //
1240         //  If ESD information available obtain reconstructed vertex and store in header.
1241         if (esdOK) {
1242             treeESD->GetEvent(nev);
1243             const AliESDVertex* esdVertex = esd->GetPrimaryVertex();
1244             Double_t position[3];
1245             esdVertex->GetXYZ(position);
1246             AliGenEventHeader* mcHeader = new  AliGenEventHeader("ESD");
1247             TArrayF mcV;
1248             mcV.Set(3);
1249             for (Int_t i = 0; i < 3; i++) mcV[i] = position[i];
1250             mcHeader->SetPrimaryVertex(mcV);
1251             header->Reset(0,nev);
1252             header->SetGenEventHeader(mcHeader);
1253             printf("***** Saved vertex %f %f %f \n", position[0], position[1], position[2]);
1254         }
1255         nev++;
1256 //
1257 //      Finish the event
1258         runLoader->TreeE()->Fill();
1259         runLoader->SetNextEvent();
1260     } // events
1261  
1262     delete rawReader;
1263 //
1264 //  Finish the run 
1265     runLoader->CdGAFile();
1266     runLoader->WriteHeader("OVERWRITE");
1267     runLoader->WriteRunLoader();
1268
1269     return kTRUE;
1270 }