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