]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/AliReconstruction.cxx
Change from Int_t/SHort_t to smaller data types, the change from Float_t/Double_t...
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.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 the reconstruction                                      //
21 //                                                                           //
22 // Clusters and tracks are created for all detectors and all events by       //
23 // typing:                                                                   //
24 //                                                                           //
25 //   AliReconstruction rec;                                                  //
26 //   rec.Run();                                                              //
27 //                                                                           //
28 // The Run method returns kTRUE in case of successful execution.             //
29 //                                                                           //
30 // If the input to the reconstruction are not simulated digits but raw data, //
31 // this can be specified by an argument of the Run method or by the method   //
32 //                                                                           //
33 //   rec.SetInput("...");                                                    //
34 //                                                                           //
35 // The input formats and the corresponding argument are:                     //
36 // - DDL raw data files: directory name, ends with "/"                       //
37 // - raw data root file: root file name, extension ".root"                   //
38 // - raw data DATE file: DATE file name, any other non-empty string          //
39 // - MC root files     : empty string, default                               //
40 //                                                                           //
41 // By default all events are reconstructed. The reconstruction can be        //
42 // limited to a range of events by giving the index of the first and the     //
43 // last event as an argument to the Run method or by calling                 //
44 //                                                                           //
45 //   rec.SetEventRange(..., ...);                                            //
46 //                                                                           //
47 // The index -1 (default) can be used for the last event to indicate no      //
48 // upper limit of the event range.                                           //
49 //                                                                           //
50 // In case of raw-data reconstruction the user can modify the default        //
51 // number of events per digits/clusters/tracks file. In case the option      //
52 // is not used the number is set 1. In case the user provides 0, than        //
53 // the number of events is equal to the number of events inside the          //
54 // raw-data file (i.e. one digits/clusters/tracks file):                     //
55 //                                                                           //
56 //   rec.SetNumberOfEventsPerFile(...);                                      //
57 //                                                                           //
58 //                                                                           //
59 // The name of the galice file can be changed from the default               //
60 // "galice.root" by passing it as argument to the AliReconstruction          //
61 // constructor or by                                                         //
62 //                                                                           //
63 //   rec.SetGAliceFile("...");                                               //
64 //                                                                           //
65 // The local reconstruction can be switched on or off for individual         //
66 // detectors by                                                              //
67 //                                                                           //
68 //   rec.SetRunLocalReconstruction("...");                                   //
69 //                                                                           //
70 // The argument is a (case sensitive) string with the names of the           //
71 // detectors separated by a space. The special string "ALL" selects all      //
72 // available detectors. This is the default.                                 //
73 //                                                                           //
74 // The reconstruction of the primary vertex position can be switched off by  //
75 //                                                                           //
76 //   rec.SetRunVertexFinder(kFALSE);                                         //
77 //                                                                           //
78 // The tracking and the creation of ESD tracks can be switched on for        //
79 // selected detectors by                                                     //
80 //                                                                           //
81 //   rec.SetRunTracking("...");                                              //
82 //                                                                           //
83 // Uniform/nonuniform field tracking switches (default: uniform field)       //
84 //                                                                           //
85 //   rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
86 //                                                                           //
87 // The filling of additional ESD information can be steered by               //
88 //                                                                           //
89 //   rec.SetFillESD("...");                                                  //
90 //                                                                           //
91 // Again, for both methods the string specifies the list of detectors.       //
92 // The default is "ALL".                                                     //
93 //                                                                           //
94 // The call of the shortcut method                                           //
95 //                                                                           //
96 //   rec.SetRunReconstruction("...");                                        //
97 //                                                                           //
98 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
99 // SetFillESD with the same detector selecting string as argument.           //
100 //                                                                           //
101 // The reconstruction requires digits or raw data as input. For the creation //
102 // of digits and raw data have a look at the class AliSimulation.            //
103 //                                                                           //
104 // For debug purposes the method SetCheckPointLevel can be used. If the      //
105 // argument is greater than 0, files with ESD events will be written after   //
106 // selected steps of the reconstruction for each event:                      //
107 //   level 1: after tracking and after filling of ESD (final)                //
108 //   level 2: in addition after each tracking step                           //
109 //   level 3: in addition after the filling of ESD for each detector         //
110 // If a final check point file exists for an event, this event will be       //
111 // skipped in the reconstruction. The tracking and the filling of ESD for    //
112 // a detector will be skipped as well, if the corresponding check point      //
113 // file exists. The ESD event will then be loaded from the file instead.     //
114 //                                                                           //
115 ///////////////////////////////////////////////////////////////////////////////
116
117 #include <TArrayF.h>
118 #include <TFile.h>
119 #include <TSystem.h>
120 #include <TROOT.h>
121 #include <TPluginManager.h>
122 #include <TGeoManager.h>
123 #include <TLorentzVector.h>
124
125 #include "AliReconstruction.h"
126 #include "AliCodeTimer.h"
127 #include "AliReconstructor.h"
128 #include "AliLog.h"
129 #include "AliRunLoader.h"
130 #include "AliRun.h"
131 #include "AliRawReaderFile.h"
132 #include "AliRawReaderDate.h"
133 #include "AliRawReaderRoot.h"
134 #include "AliRawEventHeaderBase.h"
135 #include "AliESDEvent.h"
136 #include "AliESDMuonTrack.h"
137 #include "AliESDfriend.h"
138 #include "AliESDVertex.h"
139 #include "AliESDcascade.h"
140 #include "AliESDkink.h"
141 #include "AliESDtrack.h"
142 #include "AliESDCaloCluster.h"
143 #include "AliMultiplicity.h"
144 #include "AliTracker.h"
145 #include "AliVertexer.h"
146 #include "AliVertexerTracks.h"
147 #include "AliV0vertexer.h"
148 #include "AliCascadeVertexer.h"
149 #include "AliHeader.h"
150 #include "AliGenEventHeader.h"
151 #include "AliPID.h"
152 #include "AliESDpid.h"
153 #include "AliESDtrack.h"
154
155 #include "AliESDTagCreator.h"
156 #include "AliAODTagCreator.h"
157
158 #include "AliGeomManager.h"
159 #include "AliTrackPointArray.h"
160 #include "AliCDBManager.h"
161 #include "AliCDBEntry.h"
162 #include "AliAlignObj.h"
163
164 #include "AliCentralTrigger.h"
165 #include "AliCTPRawStream.h"
166
167 #include "AliAODEvent.h"
168 #include "AliAODHeader.h"
169 #include "AliAODTrack.h"
170 #include "AliAODVertex.h"
171 #include "AliAODCluster.h"
172
173 #include "AliQualAssDataMaker.h" 
174
175 //#include "TMemStatManager.h" // memory snapshots
176 #include "AliSysInfo.h" // memory snapshots
177
178 ClassImp(AliReconstruction)
179
180
181 //_____________________________________________________________________________
182 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
183
184 //_____________________________________________________________________________
185 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
186                                      const char* name, const char* title) :
187   TNamed(name, title),
188
189   fUniformField(kTRUE),
190   fRunVertexFinder(kTRUE),
191   fRunHLTTracking(kFALSE),
192   fRunMuonTracking(kFALSE),
193   fRunV0Finder(kTRUE),
194   fRunCascadeFinder(kTRUE),
195   fStopOnError(kFALSE),
196   fWriteAlignmentData(kFALSE),
197   fWriteESDfriend(kFALSE),
198   fWriteAOD(kFALSE),
199   fFillTriggerESD(kTRUE),
200
201   fCleanESD(kTRUE),
202   fDmax(50.),
203   fZmax(50.),
204
205   fRunLocalReconstruction("ALL"),
206   fRunTracking("ALL"),
207   fFillESD("ALL"),
208   fUseTrackingErrorsForAlignment(""),
209   fGAliceFileName(gAliceFilename),
210   fInput(""),
211   fEquipIdMap(""),
212   fFirstEvent(0),
213   fLastEvent(-1),
214   fNumberOfEventsPerFile(1),
215   fCheckPointLevel(0),
216   fOptions(),
217   fLoadAlignFromCDB(kTRUE),
218   fLoadAlignData("ALL"),
219   fESDPar(""),
220
221   fRunLoader(NULL),
222   fRawReader(NULL),
223
224   fVertexer(NULL),
225   fDiamondProfile(NULL),
226
227   fAlignObjArray(NULL),
228   fCDBUri(cdbUri),
229   fRemoteCDBUri(""),
230   fSpecCDBUri()
231 {
232 // create reconstruction object with default parameters
233   
234   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
235     fReconstructor[iDet] = NULL;
236     fLoader[iDet] = NULL;
237     fTracker[iDet] = NULL;
238     fQualAssDataMaker[iDet] = NULL;
239         fQACycles[iDet] = 999999;       
240   }
241   AliPID pid;
242 }
243
244 //_____________________________________________________________________________
245 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
246   TNamed(rec),
247
248   fUniformField(rec.fUniformField),
249   fRunVertexFinder(rec.fRunVertexFinder),
250   fRunHLTTracking(rec.fRunHLTTracking),
251   fRunMuonTracking(rec.fRunMuonTracking),
252   fRunV0Finder(rec.fRunV0Finder),
253   fRunCascadeFinder(rec.fRunCascadeFinder),
254   fStopOnError(rec.fStopOnError),
255   fWriteAlignmentData(rec.fWriteAlignmentData),
256   fWriteESDfriend(rec.fWriteESDfriend),
257   fWriteAOD(rec.fWriteAOD),
258   fFillTriggerESD(rec.fFillTriggerESD),
259
260   fCleanESD(rec.fCleanESD),
261   fDmax(rec.fDmax),
262   fZmax(rec.fZmax),
263
264   fRunLocalReconstruction(rec.fRunLocalReconstruction),
265   fRunTracking(rec.fRunTracking),
266   fFillESD(rec.fFillESD),
267   fUseTrackingErrorsForAlignment(rec.fUseTrackingErrorsForAlignment),
268   fGAliceFileName(rec.fGAliceFileName),
269   fInput(rec.fInput),
270   fEquipIdMap(rec.fEquipIdMap),
271   fFirstEvent(rec.fFirstEvent),
272   fLastEvent(rec.fLastEvent),
273   fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
274   fCheckPointLevel(0),
275   fOptions(),
276   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
277   fLoadAlignData(rec.fLoadAlignData),
278   fESDPar(rec.fESDPar),
279
280   fRunLoader(NULL),
281   fRawReader(NULL),
282
283   fVertexer(NULL),
284   fDiamondProfile(NULL),
285
286   fAlignObjArray(rec.fAlignObjArray),
287   fCDBUri(rec.fCDBUri),
288   fRemoteCDBUri(rec.fRemoteCDBUri),
289   fSpecCDBUri()
290 {
291 // copy constructor
292
293   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
294     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
295   }
296   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
297     fReconstructor[iDet] = NULL;
298     fLoader[iDet] = NULL;
299     fTracker[iDet] = NULL;
300     fQualAssDataMaker[iDet] = NULL;
301         fQACycles[iDet] = rec.fQACycles[iDet];  
302   }
303   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
304     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
305   }
306 }
307
308 //_____________________________________________________________________________
309 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
310 {
311 // assignment operator
312
313   this->~AliReconstruction();
314   new(this) AliReconstruction(rec);
315   return *this;
316 }
317
318 //_____________________________________________________________________________
319 AliReconstruction::~AliReconstruction()
320 {
321 // clean up
322
323   CleanUp();
324   fOptions.Delete();
325   fSpecCDBUri.Delete();
326
327   AliCodeTimer::Instance()->Print();
328 }
329
330 //_____________________________________________________________________________
331 void AliReconstruction::InitCDBStorage()
332 {
333 // activate a default CDB storage
334 // First check if we have any CDB storage set, because it is used 
335 // to retrieve the calibration and alignment constants
336
337   AliCDBManager* man = AliCDBManager::Instance();
338   if (man->IsDefaultStorageSet())
339   {
340     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
341     AliWarning("Default CDB storage has been already set !");
342     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
343     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
344     fCDBUri = "";
345   }
346   else {
347     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
348     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
349     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
350     man->SetDefaultStorage(fCDBUri);
351   }
352
353   // Remote storage (the Grid storage) is used if it is activated
354   // and if the object is not found in the default storage
355   // OBSOLETE: Removed
356   //  if (man->IsRemoteStorageSet())
357   //  {
358   //    AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
359   //    AliWarning("Remote CDB storage has been already set !");
360   //    AliWarning(Form("Ignoring the remote storage declared in AliReconstruction: %s",fRemoteCDBUri.Data()));
361   //    AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
362   //    fRemoteCDBUri = "";
363   //  }
364   //  else {
365   //    AliDebug(2,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
366   //    AliDebug(2, Form("Remote CDB storage is set to: %s",fRemoteCDBUri.Data()));
367   //    AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
368   //    man->SetRemoteStorage(fRemoteCDBUri);
369   //  }
370
371   // Now activate the detector specific CDB storage locations
372   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
373     TObject* obj = fSpecCDBUri[i];
374     if (!obj) continue;
375     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
376     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
377     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
378     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
379   }
380   man->Print();
381 }
382
383 //_____________________________________________________________________________
384 void AliReconstruction::SetDefaultStorage(const char* uri) {
385 // Store the desired default CDB storage location
386 // Activate it later within the Run() method
387
388   fCDBUri = uri;
389
390 }
391
392 //_____________________________________________________________________________
393 void AliReconstruction::SetRemoteStorage(const char* uri) {
394 // Store the desired remote CDB storage location
395 // Activate it later within the Run() method
396 // Remote storage (the Grid storage) is used if it is activated
397 // and if the object is not found in the default storage
398
399   fRemoteCDBUri = uri;
400
401 }
402
403 //_____________________________________________________________________________
404 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
405 // Store a detector-specific CDB storage location
406 // Activate it later within the Run() method
407
408   AliCDBPath aPath(calibType);
409   if(!aPath.IsValid()){
410         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
411         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
412                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
413                         aPath.SetPath(Form("%s/*", calibType));
414                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
415                         break;
416                 }
417         }
418         if(!aPath.IsValid()){
419                 AliError(Form("Not a valid path or detector: %s", calibType));
420                 return;
421         }
422   }
423
424 //  // check that calibType refers to a "valid" detector name
425 //  Bool_t isDetector = kFALSE;
426 //  for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
427 //    TString detName = fgkDetectorName[iDet];
428 //    if(aPath.GetLevel0() == detName) {
429 //      isDetector = kTRUE;
430 //      break;
431 //    }
432 //  }
433 //
434 //  if(!isDetector) {
435 //      AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
436 //      return;
437 //  }
438
439   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
440   if (obj) fSpecCDBUri.Remove(obj);
441   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
442
443 }
444
445
446
447
448 //_____________________________________________________________________________
449 Bool_t AliReconstruction::SetRunNumber()
450 {
451   // The method is called in Run() in order
452   // to set a correct run number.
453   // In case of raw data reconstruction the
454   // run number is taken from the raw data header
455
456   if(AliCDBManager::Instance()->GetRun() < 0) {
457     if (!fRunLoader) {
458       AliError("No run loader is found !"); 
459       return kFALSE;
460     }
461     // read run number from gAlice
462     if(fRunLoader->GetAliRun())
463       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
464     else {
465       if(fRawReader) {
466         if(fRawReader->NextEvent()) {
467           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
468           fRawReader->RewindEvents();
469         }
470         else {
471           AliError("No raw-data events found !");
472           return kFALSE;
473         }
474       }
475       else {
476         AliError("Neither gAlice nor RawReader objects are found !");
477         return kFALSE;
478       }
479     }
480     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
481   }
482   return kTRUE;
483 }
484
485 //_____________________________________________________________________________
486 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
487 {
488   // Read the alignment objects from CDB.
489   // Each detector is supposed to have the
490   // alignment objects in DET/Align/Data CDB path.
491   // All the detector objects are then collected,
492   // sorted by geometry level (starting from ALIC) and
493   // then applied to the TGeo geometry.
494   // Finally an overlaps check is performed.
495
496   // Load alignment data from CDB and fill fAlignObjArray 
497   if(fLoadAlignFromCDB){
498         
499     TString detStr = detectors;
500     TString loadAlObjsListOfDets = "";
501     
502     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
503       if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
504       loadAlObjsListOfDets += fgkDetectorName[iDet];
505       loadAlObjsListOfDets += " ";
506     } // end loop over detectors
507     loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
508     AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
509   }else{
510     // Check if the array with alignment objects was
511     // provided by the user. If yes, apply the objects
512     // to the present TGeo geometry
513     if (fAlignObjArray) {
514       if (gGeoManager && gGeoManager->IsClosed()) {
515         if (AliGeomManager::ApplyAlignObjsToGeom(*fAlignObjArray) == kFALSE) {
516           AliError("The misalignment of one or more volumes failed!"
517                    "Compare the list of simulated detectors and the list of detector alignment data!");
518           return kFALSE;
519         }
520       }
521       else {
522         AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
523         return kFALSE;
524       }
525     }
526   }
527   
528   delete fAlignObjArray; fAlignObjArray=0;
529
530   return kTRUE;
531 }
532
533 //_____________________________________________________________________________
534 void AliReconstruction::SetGAliceFile(const char* fileName)
535 {
536 // set the name of the galice file
537
538   fGAliceFileName = fileName;
539 }
540
541 //_____________________________________________________________________________
542 void AliReconstruction::SetOption(const char* detector, const char* option)
543 {
544 // set options for the reconstruction of a detector
545
546   TObject* obj = fOptions.FindObject(detector);
547   if (obj) fOptions.Remove(obj);
548   fOptions.Add(new TNamed(detector, option));
549 }
550
551
552 //_____________________________________________________________________________
553 Bool_t AliReconstruction::Run(const char* input)
554 {
555 // run the reconstruction
556
557   AliCodeTimerAuto("")
558   
559   // set the input
560   if (!input) input = fInput.Data();
561   TString fileName(input);
562   if (fileName.EndsWith("/")) {
563     fRawReader = new AliRawReaderFile(fileName);
564   } else if (fileName.EndsWith(".root")) {
565     fRawReader = new AliRawReaderRoot(fileName);
566   } else if (!fileName.IsNull()) {
567     fRawReader = new AliRawReaderDate(fileName);
568     fRawReader->SelectEvents(7);
569   }
570   if (!fEquipIdMap.IsNull() && fRawReader)
571     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
572
573    AliSysInfo::AddStamp("Start");
574   // get the run loader
575   if (!InitRunLoader()) return kFALSE;
576    AliSysInfo::AddStamp("LoadLoader");
577
578   // Initialize the CDB storage
579   InitCDBStorage();
580    AliSysInfo::AddStamp("LoadCDB");
581
582   // Set run number in CDBManager (if it is not already set by the user)
583   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
584
585   // Import ideal TGeo geometry and apply misalignment
586   if (!gGeoManager) {
587     TString geom(gSystem->DirName(fGAliceFileName));
588     geom += "/geometry.root";
589     AliGeomManager::LoadGeometry(geom.Data());
590     if (!gGeoManager) if (fStopOnError) return kFALSE;
591   }
592
593   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
594    AliSysInfo::AddStamp("LoadGeom");
595
596   // local reconstruction
597   if (!fRunLocalReconstruction.IsNull()) {
598     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
599       if (fStopOnError) {CleanUp(); return kFALSE;}
600     }
601   }
602 //  if (!fRunVertexFinder && fRunTracking.IsNull() && 
603 //      fFillESD.IsNull()) return kTRUE;
604
605   // get vertexer
606   if (fRunVertexFinder && !CreateVertexer()) {
607     if (fStopOnError) {
608       CleanUp(); 
609       return kFALSE;
610     }
611   }
612    AliSysInfo::AddStamp("Vertexer");
613
614   // get trackers
615   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
616     if (fStopOnError) {
617       CleanUp(); 
618       return kFALSE;
619     }      
620   }
621    AliSysInfo::AddStamp("LoadTrackers");
622
623   // get the possibly already existing ESD file and tree
624   AliESDEvent* esd = new AliESDEvent(); AliESDEvent* hltesd = new AliESDEvent();
625   TFile* fileOld = NULL;
626   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
627   if (!gSystem->AccessPathName("AliESDs.root")){
628     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
629     fileOld = TFile::Open("AliESDs.old.root");
630     if (fileOld && fileOld->IsOpen()) {
631       treeOld = (TTree*) fileOld->Get("esdTree");
632       if (treeOld)esd->ReadFromTree(treeOld);
633       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
634       if (hlttreeOld)   hltesd->ReadFromTree(hlttreeOld);
635     }
636   }
637
638   // create the ESD output file and tree
639   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
640   file->SetCompressionLevel(2);
641   if (!file->IsOpen()) {
642     AliError("opening AliESDs.root failed");
643     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
644   }
645
646   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
647   esd = new AliESDEvent();
648   esd->CreateStdContent();
649   esd->WriteToTree(tree);
650
651   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
652   hltesd = new AliESDEvent();
653   hltesd->CreateStdContent();
654   hltesd->WriteToTree(hlttree);
655
656   /* CKB Why?
657   delete esd; delete hltesd;
658   esd = NULL; hltesd = NULL;
659   */
660   // create the branch with ESD additions
661
662
663
664   AliESDfriend *esdf = 0; 
665   if (fWriteESDfriend) {
666     esdf = new AliESDfriend();
667     TBranch *br=tree->Branch("ESDfriend.","AliESDfriend", &esdf);
668     br->SetFile("AliESDfriends.root");
669     esd->AddObject(esdf);
670   }
671
672   
673   // Get the diamond profile from OCDB
674   AliCDBEntry* entry = AliCDBManager::Instance()
675         ->Get("GRP/Calib/MeanVertex");
676         
677   if(entry) {
678         fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());  
679   } else {
680         AliError("No diamond profile found in OCDB!");
681   }
682
683   AliVertexerTracks tVertexer(AliTracker::GetBz());
684   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
685
686   // loop over events
687  
688   if (fRawReader) fRawReader->RewindEvents();
689    TString detStr(fFillESD) ; 
690
691   ProcInfo_t ProcInfo;
692   gSystem->GetProcInfo(&ProcInfo);
693   AliInfo(Form("Current memory usage %d %d", ProcInfo.fMemResident, ProcInfo.fMemVirtual));
694   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
695     if (fRawReader) fRawReader->NextEvent();
696     if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
697       // copy old ESD to the new one
698       if (treeOld) {
699         esd->ReadFromTree(treeOld);
700         treeOld->GetEntry(iEvent);
701       }
702       tree->Fill();
703       if (hlttreeOld) {
704         esd->ReadFromTree(hlttreeOld);
705         hlttreeOld->GetEntry(iEvent);
706       }
707       hlttree->Fill();
708       continue;
709     }
710     
711     AliInfo(Form("processing event %d", iEvent));
712     fRunLoader->GetEvent(iEvent);
713
714     char aFileName[256];
715     sprintf(aFileName, "ESD_%d.%d_final.root", 
716             fRunLoader->GetHeader()->GetRun(), 
717             fRunLoader->GetHeader()->GetEventNrInRun());
718     if (!gSystem->AccessPathName(aFileName)) continue;
719
720     // local reconstruction
721     if (!fRunLocalReconstruction.IsNull()) {
722       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
723         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
724       }
725     }
726
727   
728     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
729     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
730     esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
731     hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
732     
733     // Set magnetic field from the tracker
734     esd->SetMagneticField(AliTracker::GetBz());
735     hltesd->SetMagneticField(AliTracker::GetBz());
736
737     
738     
739     // Fill raw-data error log into the ESD
740     if (fRawReader) FillRawDataErrorLog(iEvent,esd);
741
742     // vertex finder
743     if (fRunVertexFinder) {
744       if (!ReadESD(esd, "vertex")) {
745         if (!RunVertexFinder(esd)) {
746           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
747         }
748         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
749       }
750     }
751
752     // HLT tracking
753     if (!fRunTracking.IsNull()) {
754       if (fRunHLTTracking) {
755         hltesd->SetVertex(esd->GetVertex());
756         if (!RunHLTTracking(hltesd)) {
757           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
758         }
759       }
760     }
761
762     // Muon tracking
763     if (!fRunTracking.IsNull()) {
764       if (fRunMuonTracking) {
765         if (!RunMuonTracking(esd)) {
766           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
767         }
768       }
769     }
770
771     // barrel tracking
772     if (!fRunTracking.IsNull()) {
773       if (!ReadESD(esd, "tracking")) {
774         if (!RunTracking(esd)) {
775           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
776         }
777         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
778       }
779     }
780
781     // fill ESD
782     if (!fFillESD.IsNull()) {
783       if (!FillESD(esd, fFillESD)) {
784         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
785       }
786     }
787   
788     if (!fFillESD.IsNull()) 
789     RunQualAss(fFillESD.Data(), esd);
790
791     // fill Event header information from the RawEventHeader
792     if (fRawReader){FillRawEventHeaderESD(esd);}
793
794     // combined PID
795     AliESDpid::MakePID(esd);
796     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
797
798     if (fFillTriggerESD) {
799       if (!ReadESD(esd, "trigger")) {
800         if (!FillTriggerESD(esd)) {
801           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
802         }
803         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
804       }
805     }
806
807     file->cd();
808
809     //Try to improve the reconstructed primary vertex position using the tracks
810     AliESDVertex *pvtx=0;
811     Bool_t dovertex=kTRUE;
812     TObject* obj = fOptions.FindObject("ITS");
813     if (obj) {
814       TString optITS = obj->GetTitle();
815       if (optITS.Contains("cosmics") || optITS.Contains("COSMICS")) 
816         dovertex=kFALSE;
817     }
818     if(dovertex) pvtx=tVertexer.FindPrimaryVertex(esd);
819     if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
820     
821     if (pvtx)
822     if (pvtx->GetStatus()) {
823        // Store the improved primary vertex
824        esd->SetPrimaryVertex(pvtx);
825        // Propagate the tracks to the DCA to the improved primary vertex
826        Double_t somethingbig = 777.;
827        Double_t bz = esd->GetMagneticField();
828        Int_t nt=esd->GetNumberOfTracks();
829        while (nt--) {
830          AliESDtrack *t = esd->GetTrack(nt);
831          t->RelateToVertex(pvtx, bz, somethingbig);
832        } 
833     }
834
835     if (fRunV0Finder) {
836        // V0 finding
837        AliV0vertexer vtxer;
838        vtxer.Tracks2V0vertices(esd);
839
840        if (fRunCascadeFinder) {
841           // Cascade finding
842           AliCascadeVertexer cvtxer;
843           cvtxer.V0sTracks2CascadeVertices(esd);
844        }
845     }
846  
847     // write ESD
848     if (fCleanESD) CleanESD(esd);
849     if (fWriteESDfriend) {
850       esdf->~AliESDfriend();
851       new (esdf) AliESDfriend(); // Reset...
852       esd->GetESDfriend(esdf);
853     }
854     tree->Fill();
855
856     // write HLT ESD
857     hlttree->Fill();
858
859     if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
860     esd->Reset();
861     hltesd->Reset();
862     if (fWriteESDfriend) {
863       esdf->~AliESDfriend();
864       new (esdf) AliESDfriend(); // Reset...
865     }
866     // esdf->Reset();
867     // delete esdf; esdf = 0;
868   // ESD QA 
869  
870     gSystem->GetProcInfo(&ProcInfo);
871     AliInfo(Form("Event %d -> Current memory usage %d %d",iEvent, ProcInfo.fMemResident, ProcInfo.fMemVirtual));
872   }
873   
874   detStr = fFillESD ; 
875   // write quality assurance ESDs data (one entry for all events)
876   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
877         if (!IsSelected(fgkDetectorName[iDet], detStr)) 
878                 continue;
879     AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
880     if (!qadm) continue;
881     
882     qadm->EndOfCycle(AliQualAss::kRECPOINTS);
883     qadm->EndOfCycle(AliQualAss::kESDS);
884     qadm->Finish(AliQualAss::kRECPOINTS);
885     qadm->Finish(AliQualAss::kESDS) ; 
886     
887   }
888
889   tree->GetUserInfo()->Add(esd);
890   hlttree->GetUserInfo()->Add(hltesd);
891
892
893
894   if(fESDPar.Contains("ESD.par")){
895     AliInfo("Attaching ESD.par to Tree");
896     TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
897     tree->GetUserInfo()->Add(fn);
898   }
899
900
901   file->cd();
902   if (fWriteESDfriend)
903     tree->SetBranchStatus("ESDfriend*",0);
904   // we want to have only one tree version number
905   tree->Write(tree->GetName(),TObject::kOverwrite);
906   hlttree->Write();
907
908   if (fWriteAOD) {
909     TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
910     ESDFile2AODFile(file, aodFile);
911     aodFile->Close();
912   }
913
914   gROOT->cd();
915   CleanUp(file, fileOld);
916   
917   // Create tags for the events in the ESD tree (the ESD tree is always present)
918   // In case of empty events the tags will contain dummy values
919   AliESDTagCreator *esdtagCreator = new AliESDTagCreator();
920   esdtagCreator->CreateESDTags(fFirstEvent,fLastEvent);
921   if (fWriteAOD) {
922     AliAODTagCreator *aodtagCreator = new AliAODTagCreator();
923     aodtagCreator->CreateAODTags(fFirstEvent,fLastEvent);
924   }
925
926   return kTRUE;
927 }
928
929
930 //_____________________________________________________________________________
931 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
932 {
933 // run the local reconstruction
934   static Int_t eventNr=0;
935   AliCodeTimerAuto("")
936
937  //  AliCDBManager* man = AliCDBManager::Instance();
938 //   Bool_t origCache = man->GetCacheFlag();
939
940 //   TString detStr = detectors;
941 //   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
942 //     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
943 //     AliReconstructor* reconstructor = GetReconstructor(iDet);
944 //     if (!reconstructor) continue;
945 //     if (reconstructor->HasLocalReconstruction()) continue;
946
947 //     AliCodeTimerStart(Form("running reconstruction for %s", fgkDetectorName[iDet]));
948 //     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
949     
950 //     AliCodeTimerStart(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));                          
951 //     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
952
953 //     man->SetCacheFlag(kTRUE);
954 //     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
955 //     man->GetAll(calibPath); // entries are cached!
956
957 //     AliCodeTimerStop(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
958      
959 //     if (fRawReader) {
960 //       fRawReader->RewindEvents();
961 //       reconstructor->Reconstruct(fRunLoader, fRawReader);
962 //     } else {
963 //       reconstructor->Reconstruct(fRunLoader);
964 //     }
965      
966 //      AliCodeTimerStop(Form("running reconstruction for %s", fgkDetectorName[iDet]));
967     // AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr));
968
969 //     // unload calibration data
970 //     man->UnloadFromCache(calibPath);
971 //     //man->ClearCache();
972 //   }
973
974 //   man->SetCacheFlag(origCache);
975
976 //   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
977 //     AliError(Form("the following detectors were not found: %s",
978 //                   detStr.Data()));
979 //     if (fStopOnError) return kFALSE;
980 //   }
981
982           eventNr++;
983   return kTRUE;
984 }
985
986 //_____________________________________________________________________________
987 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
988 {
989 // run the local reconstruction
990   static Int_t eventNr=0;
991   AliCodeTimerAuto("")
992
993   TString detStr = detectors;
994   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
995     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
996     AliReconstructor* reconstructor = GetReconstructor(iDet);
997     if (!reconstructor) continue;
998     AliLoader* loader = fLoader[iDet];
999     if (!loader) {
1000       AliWarning(Form("No loader is defined for %s!",fgkDetectorName[iDet]));
1001       continue;
1002     }
1003
1004     // conversion of digits
1005     if (fRawReader && reconstructor->HasDigitConversion()) {
1006       AliInfo(Form("converting raw data digits into root objects for %s", 
1007                    fgkDetectorName[iDet]));
1008       AliCodeTimerAuto(Form("converting raw data digits into root objects for %s", 
1009                             fgkDetectorName[iDet]));
1010       loader->LoadDigits("update");
1011       loader->CleanDigits();
1012       loader->MakeDigitsContainer();
1013       TTree* digitsTree = loader->TreeD();
1014       reconstructor->ConvertDigits(fRawReader, digitsTree);
1015       loader->WriteDigits("OVERWRITE");
1016       loader->UnloadDigits();
1017     }
1018
1019     // local reconstruction
1020     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1021     AliCodeTimerAuto(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1022     loader->LoadRecPoints("update");
1023     loader->CleanRecPoints();
1024     loader->MakeRecPointsContainer();
1025     TTree* clustersTree = loader->TreeR();
1026     if (fRawReader && !reconstructor->HasDigitConversion()) {
1027       reconstructor->Reconstruct(fRawReader, clustersTree);
1028     } else {
1029       loader->LoadDigits("read");
1030       TTree* digitsTree = loader->TreeD();
1031       if (!digitsTree) {
1032         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
1033         if (fStopOnError) return kFALSE;
1034       } else {
1035         reconstructor->Reconstruct(digitsTree, clustersTree);
1036       }
1037       loader->UnloadDigits();
1038     }
1039
1040     AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
1041     if (qadm) {
1042       AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1043       AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1044           
1045       if (qadm->IsCycleDone() ) {
1046         
1047         qadm->EndOfCycle(AliQualAss::kRECPOINTS) ; 
1048         qadm->EndOfCycle(AliQualAss::kESDS) ; 
1049       qadm->StartOfCycle(AliQualAss::kRECPOINTS) ; 
1050       qadm->StartOfCycle(AliQualAss::kESDS, "same") ; 
1051    }
1052       qadm->Exec(AliQualAss::kRECPOINTS, clustersTree) ; 
1053       AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1054     }
1055
1056     loader->WriteRecPoints("OVERWRITE");
1057     loader->UnloadRecPoints();
1058      AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr));
1059   }
1060
1061   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1062     AliError(Form("the following detectors were not found: %s",
1063                   detStr.Data()));
1064     if (fStopOnError) return kFALSE;
1065   }
1066   eventNr++;
1067   return kTRUE;
1068 }
1069
1070 //_____________________________________________________________________________
1071 Bool_t AliReconstruction::RunVertexFinder(AliESDEvent*& esd)
1072 {
1073 // run the barrel tracking
1074
1075   AliCodeTimerAuto("")
1076
1077   AliESDVertex* vertex = NULL;
1078   Double_t vtxPos[3] = {0, 0, 0};
1079   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
1080   TArrayF mcVertex(3); 
1081   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
1082     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
1083     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
1084   }
1085
1086   if (fVertexer) {
1087     if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
1088     AliInfo("running the ITS vertex finder");
1089     if (fLoader[0]) fLoader[0]->LoadRecPoints();
1090     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
1091     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
1092     if(!vertex){
1093       AliWarning("Vertex not found");
1094       vertex = new AliESDVertex();
1095       vertex->SetName("default");
1096     }
1097     else {
1098       vertex->SetName("reconstructed");
1099     }
1100
1101   } else {
1102     AliInfo("getting the primary vertex from MC");
1103     vertex = new AliESDVertex(vtxPos, vtxErr);
1104   }
1105
1106   if (vertex) {
1107     vertex->GetXYZ(vtxPos);
1108     vertex->GetSigmaXYZ(vtxErr);
1109   } else {
1110     AliWarning("no vertex reconstructed");
1111     vertex = new AliESDVertex(vtxPos, vtxErr);
1112   }
1113   esd->SetVertex(vertex);
1114   // if SPD multiplicity has been determined, it is stored in the ESD
1115   AliMultiplicity *mult = fVertexer->GetMultiplicity();
1116   if(mult)esd->SetMultiplicity(mult);
1117
1118   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1119     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1120   }  
1121   delete vertex;
1122
1123   return kTRUE;
1124 }
1125
1126 //_____________________________________________________________________________
1127 Bool_t AliReconstruction::RunHLTTracking(AliESDEvent*& esd)
1128 {
1129 // run the HLT barrel tracking
1130
1131   AliCodeTimerAuto("")
1132
1133   if (!fRunLoader) {
1134     AliError("Missing runLoader!");
1135     return kFALSE;
1136   }
1137
1138   AliInfo("running HLT tracking");
1139
1140   // Get a pointer to the HLT reconstructor
1141   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1142   if (!reconstructor) return kFALSE;
1143
1144   // TPC + ITS
1145   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1146     TString detName = fgkDetectorName[iDet];
1147     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1148     reconstructor->SetOption(detName.Data());
1149     AliTracker *tracker = reconstructor->CreateTracker();
1150     if (!tracker) {
1151       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1152       if (fStopOnError) return kFALSE;
1153       continue;
1154     }
1155     Double_t vtxPos[3];
1156     Double_t vtxErr[3]={0.005,0.005,0.010};
1157     const AliESDVertex *vertex = esd->GetVertex();
1158     vertex->GetXYZ(vtxPos);
1159     tracker->SetVertex(vtxPos,vtxErr);
1160     if(iDet != 1) {
1161       fLoader[iDet]->LoadRecPoints("read");
1162       TTree* tree = fLoader[iDet]->TreeR();
1163       if (!tree) {
1164         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1165         return kFALSE;
1166       }
1167       tracker->LoadClusters(tree);
1168     }
1169     if (tracker->Clusters2Tracks(esd) != 0) {
1170       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1171       return kFALSE;
1172     }
1173     if(iDet != 1) {
1174       tracker->UnloadClusters();
1175     }
1176     delete tracker;
1177   }
1178
1179   return kTRUE;
1180 }
1181
1182 //_____________________________________________________________________________
1183 Bool_t AliReconstruction::RunMuonTracking(AliESDEvent*& esd)
1184 {
1185 // run the muon spectrometer tracking
1186
1187   AliCodeTimerAuto("")
1188
1189   if (!fRunLoader) {
1190     AliError("Missing runLoader!");
1191     return kFALSE;
1192   }
1193   Int_t iDet = 7; // for MUON
1194
1195   AliInfo("is running...");
1196
1197   // Get a pointer to the MUON reconstructor
1198   AliReconstructor *reconstructor = GetReconstructor(iDet);
1199   if (!reconstructor) return kFALSE;
1200
1201   
1202   TString detName = fgkDetectorName[iDet];
1203   AliDebug(1, Form("%s tracking", detName.Data()));
1204   AliTracker *tracker =  reconstructor->CreateTracker();
1205   if (!tracker) {
1206     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1207     return kFALSE;
1208   }
1209      
1210   // create Tracks
1211   fLoader[iDet]->LoadTracks("update");
1212   fLoader[iDet]->CleanTracks();
1213   fLoader[iDet]->MakeTracksContainer();
1214
1215   // read RecPoints
1216   fLoader[iDet]->LoadRecPoints("read");  
1217   tracker->LoadClusters(fLoader[iDet]->TreeR());
1218   
1219   Int_t rv = tracker->Clusters2Tracks(esd);
1220   
1221   fLoader[iDet]->UnloadRecPoints();
1222   
1223   if ( rv )
1224   {
1225     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1226     return kFALSE;
1227   }
1228   
1229   tracker->UnloadClusters();
1230   
1231   fLoader[iDet]->UnloadRecPoints();
1232
1233   fLoader[iDet]->WriteTracks("OVERWRITE");
1234   fLoader[iDet]->UnloadTracks();
1235
1236   delete tracker;
1237   
1238   return kTRUE;
1239 }
1240
1241
1242 //_____________________________________________________________________________
1243 Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd)
1244 {
1245 // run the barrel tracking
1246   static Int_t eventNr=0;
1247   AliCodeTimerAuto("")
1248
1249   AliInfo("running tracking");
1250
1251   //Fill the ESD with the T0 info (will be used by the TOF) 
1252   if (fReconstructor[11] && fLoader[11]) {
1253     fLoader[11]->LoadRecPoints("READ");
1254     TTree *treeR = fLoader[11]->TreeR();
1255     GetReconstructor(11)->FillESD((TTree *)NULL,treeR,esd);
1256   }
1257
1258   // pass 1: TPC + ITS inwards
1259   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1260     if (!fTracker[iDet]) continue;
1261     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1262
1263     // load clusters
1264     fLoader[iDet]->LoadRecPoints("read");
1265     TTree* tree = fLoader[iDet]->TreeR();
1266     if (!tree) {
1267       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1268       return kFALSE;
1269     }
1270     fTracker[iDet]->LoadClusters(tree);
1271
1272     // run tracking
1273     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1274       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1275       return kFALSE;
1276     }
1277     if (fCheckPointLevel > 1) {
1278       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1279     }
1280     // preliminary PID in TPC needed by the ITS tracker
1281     if (iDet == 1) {
1282       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1283       AliESDpid::MakePID(esd);
1284     } 
1285      AliSysInfo::AddStamp(Form("Tracking0%s_%d",fgkDetectorName[iDet],eventNr));
1286   }
1287
1288   // pass 2: ALL backwards
1289   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1290     if (!fTracker[iDet]) continue;
1291     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1292
1293     // load clusters
1294     if (iDet > 1) {     // all except ITS, TPC
1295       TTree* tree = NULL;
1296       fLoader[iDet]->LoadRecPoints("read");
1297       tree = fLoader[iDet]->TreeR();
1298       if (!tree) {
1299         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1300         return kFALSE;
1301       }
1302       fTracker[iDet]->LoadClusters(tree); 
1303        AliSysInfo::AddStamp(Form("LoadCluster0%s_%d",fgkDetectorName[iDet],eventNr));
1304     }
1305
1306     // run tracking
1307     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1308       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1309       //      return kFALSE;
1310     }
1311     if (fCheckPointLevel > 1) {
1312       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1313     }
1314
1315     // unload clusters
1316     if (iDet > 2) {     // all except ITS, TPC, TRD
1317       fTracker[iDet]->UnloadClusters();
1318       fLoader[iDet]->UnloadRecPoints();
1319     }
1320     // updated PID in TPC needed by the ITS tracker -MI
1321     if (iDet == 1) {
1322       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1323       AliESDpid::MakePID(esd);
1324     }
1325      AliSysInfo::AddStamp(Form("Tracking1%s_%d",fgkDetectorName[iDet],eventNr));
1326   }
1327
1328   // write space-points to the ESD in case alignment data output
1329   // is switched on
1330   if (fWriteAlignmentData)
1331     WriteAlignmentData(esd);
1332
1333   // pass 3: TRD + TPC + ITS refit inwards
1334   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1335     if (!fTracker[iDet]) continue;
1336     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1337
1338     // run tracking
1339     if (fTracker[iDet]->RefitInward(esd) != 0) {
1340       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1341       //      return kFALSE;
1342     }
1343     if (fCheckPointLevel > 1) {
1344       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1345     }
1346
1347     // unload clusters
1348     fTracker[iDet]->UnloadClusters();
1349     fLoader[iDet]->UnloadRecPoints();
1350      AliSysInfo::AddStamp(Form("Tracking2%s_%d",fgkDetectorName[iDet],eventNr));
1351   }
1352   //
1353   // Propagate track to the vertex - if not done by ITS
1354   //
1355   Int_t ntracks = esd->GetNumberOfTracks();
1356   for (Int_t itrack=0; itrack<ntracks; itrack++){
1357     const Double_t kRadius  = 3;   // beam pipe radius
1358     const Double_t kMaxStep = 5;   // max step
1359     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1360     Double_t       fieldZ   = AliTracker::GetBz();  //
1361     AliESDtrack * track = esd->GetTrack(itrack);
1362     if (!track) continue;
1363     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1364    AliTracker::PropagateTrackTo(track,kRadius,track->GetMass(),kMaxStep,kTRUE);
1365     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1366   }
1367   eventNr++;
1368   return kTRUE;
1369 }
1370
1371 //_____________________________________________________________________________
1372 Bool_t AliReconstruction::CleanESD(AliESDEvent *esd){
1373   //
1374   // Remove the data which are not needed for the physics analysis.
1375   //
1376
1377   AliInfo("Cleaning the ESD...");
1378   Int_t nTracks=esd->GetNumberOfTracks();
1379   AliInfo(Form("Number of ESD tracks before cleaning %d",nTracks));
1380
1381   Float_t cleanPars[]={fDmax,fZmax};
1382   Bool_t rc=esd->Clean(cleanPars);
1383
1384   nTracks=esd->GetNumberOfTracks();
1385   AliInfo(Form("Number of ESD tracks after cleaning %d",nTracks));
1386
1387   return rc;
1388 }
1389
1390 //_____________________________________________________________________________
1391 Bool_t AliReconstruction::FillESD(AliESDEvent*& esd, const TString& detectors)
1392 {
1393 // fill the event summary data
1394
1395   AliCodeTimerAuto("")
1396     static Int_t eventNr=0; 
1397   TString detStr = detectors;
1398   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1399     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1400     AliReconstructor* reconstructor = GetReconstructor(iDet);
1401     if (!reconstructor) continue;
1402
1403     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1404       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1405       TTree* clustersTree = NULL;
1406       if (fLoader[iDet]) {
1407         fLoader[iDet]->LoadRecPoints("read");
1408         clustersTree = fLoader[iDet]->TreeR();
1409         if (!clustersTree) {
1410           AliError(Form("Can't get the %s clusters tree", 
1411                         fgkDetectorName[iDet]));
1412           if (fStopOnError) return kFALSE;
1413         }
1414       }
1415       if (fRawReader && !reconstructor->HasDigitConversion()) {
1416         reconstructor->FillESD(fRawReader, clustersTree, esd);
1417       } else {
1418         TTree* digitsTree = NULL;
1419         if (fLoader[iDet]) {
1420           fLoader[iDet]->LoadDigits("read");
1421           digitsTree = fLoader[iDet]->TreeD();
1422           if (!digitsTree) {
1423             AliError(Form("Can't get the %s digits tree", 
1424                           fgkDetectorName[iDet]));
1425             if (fStopOnError) return kFALSE;
1426           }
1427         }
1428         reconstructor->FillESD(digitsTree, clustersTree, esd);
1429         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1430       }
1431       if (fLoader[iDet]) {
1432         fLoader[iDet]->UnloadRecPoints();
1433       }
1434
1435       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1436     }
1437   }
1438
1439   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1440     AliError(Form("the following detectors were not found: %s", 
1441                   detStr.Data()));
1442     if (fStopOnError) return kFALSE;
1443   }
1444    AliSysInfo::AddStamp(Form("FillESD%d",eventNr));
1445   eventNr++;
1446   return kTRUE;
1447 }
1448
1449 //_____________________________________________________________________________
1450 Bool_t AliReconstruction::FillTriggerESD(AliESDEvent*& esd)
1451 {
1452   // Reads the trigger decision which is
1453   // stored in Trigger.root file and fills
1454   // the corresponding esd entries
1455
1456   AliCodeTimerAuto("")
1457   
1458   AliInfo("Filling trigger information into the ESD");
1459
1460   if (fRawReader) {
1461     AliCTPRawStream input(fRawReader);
1462     if (!input.Next()) {
1463       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1464       return kFALSE;
1465     }
1466     esd->SetTriggerMask(input.GetClassMask());
1467     esd->SetTriggerCluster(input.GetClusterMask());
1468   }
1469   else {
1470     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1471     if (runloader) {
1472       if (!runloader->LoadTrigger()) {
1473         AliCentralTrigger *aCTP = runloader->GetTrigger();
1474         esd->SetTriggerMask(aCTP->GetClassMask());
1475         esd->SetTriggerCluster(aCTP->GetClusterMask());
1476       }
1477       else {
1478         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1479         return kFALSE;
1480       }
1481     }
1482     else {
1483       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1484       return kFALSE;
1485     }
1486   }
1487
1488   return kTRUE;
1489 }
1490
1491
1492
1493
1494
1495 //_____________________________________________________________________________
1496 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESDEvent*& esd)
1497 {
1498   // 
1499   // Filling information from RawReader Header
1500   // 
1501
1502   AliInfo("Filling information from RawReader Header");
1503   esd->SetBunchCrossNumber(0);
1504   esd->SetOrbitNumber(0);
1505   esd->SetPeriodNumber(0);
1506   esd->SetTimeStamp(0);
1507   esd->SetEventType(0);
1508   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1509   if (eventHeader){
1510
1511     const UInt_t *id = eventHeader->GetP("Id");
1512     esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1513     esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1514     esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1515
1516     esd->SetTimeStamp((eventHeader->Get("Timestamp")));  
1517     esd->SetEventType((eventHeader->Get("Type")));
1518   }
1519
1520   return kTRUE;
1521 }
1522
1523
1524 //_____________________________________________________________________________
1525 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1526 {
1527 // check whether detName is contained in detectors
1528 // if yes, it is removed from detectors
1529
1530   // check if all detectors are selected
1531   if ((detectors.CompareTo("ALL") == 0) ||
1532       detectors.BeginsWith("ALL ") ||
1533       detectors.EndsWith(" ALL") ||
1534       detectors.Contains(" ALL ")) {
1535     detectors = "ALL";
1536     return kTRUE;
1537   }
1538
1539   // search for the given detector
1540   Bool_t result = kFALSE;
1541   if ((detectors.CompareTo(detName) == 0) ||
1542       detectors.BeginsWith(detName+" ") ||
1543       detectors.EndsWith(" "+detName) ||
1544       detectors.Contains(" "+detName+" ")) {
1545     detectors.ReplaceAll(detName, "");
1546     result = kTRUE;
1547   }
1548
1549   // clean up the detectors string
1550   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1551   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1552   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1553
1554   return result;
1555 }
1556
1557 //_____________________________________________________________________________
1558 Bool_t AliReconstruction::InitRunLoader()
1559 {
1560 // get or create the run loader
1561
1562   if (gAlice) delete gAlice;
1563   gAlice = NULL;
1564
1565   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1566     // load all base libraries to get the loader classes
1567     TString libs = gSystem->GetLibraries();
1568     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1569       TString detName = fgkDetectorName[iDet];
1570       if (detName == "HLT") continue;
1571       if (libs.Contains("lib" + detName + "base.so")) continue;
1572       gSystem->Load("lib" + detName + "base.so");
1573     }
1574     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1575     if (!fRunLoader) {
1576       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1577       CleanUp();
1578       return kFALSE;
1579     }
1580     fRunLoader->CdGAFile();
1581     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1582       if (fRunLoader->LoadgAlice() == 0) {
1583         gAlice = fRunLoader->GetAliRun();
1584         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1585       }
1586     }
1587     if (!gAlice && !fRawReader) {
1588       AliError(Form("no gAlice object found in file %s",
1589                     fGAliceFileName.Data()));
1590       CleanUp();
1591       return kFALSE;
1592     }
1593
1594     //PH This is a temporary fix to give access to the kinematics
1595     //PH that is needed for the labels of ITS clusters
1596     fRunLoader->LoadKinematics();
1597
1598   } else {               // galice.root does not exist
1599     if (!fRawReader) {
1600       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1601       CleanUp();
1602       return kFALSE;
1603     }
1604     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1605                                     AliConfig::GetDefaultEventFolderName(),
1606                                     "recreate");
1607     if (!fRunLoader) {
1608       AliError(Form("could not create run loader in file %s", 
1609                     fGAliceFileName.Data()));
1610       CleanUp();
1611       return kFALSE;
1612     }
1613     fRunLoader->MakeTree("E");
1614     Int_t iEvent = 0;
1615     while (fRawReader->NextEvent()) {
1616       fRunLoader->SetEventNumber(iEvent);
1617       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
1618                                      iEvent, iEvent);
1619       fRunLoader->MakeTree("H");
1620       fRunLoader->TreeE()->Fill();
1621       iEvent++;
1622     }
1623     fRawReader->RewindEvents();
1624     if (fNumberOfEventsPerFile > 0)
1625       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1626     else
1627       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1628     fRunLoader->WriteHeader("OVERWRITE");
1629     fRunLoader->CdGAFile();
1630     fRunLoader->Write(0, TObject::kOverwrite);
1631 //    AliTracker::SetFieldMap(???);
1632   }
1633
1634   return kTRUE;
1635 }
1636
1637 //_____________________________________________________________________________
1638 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1639 {
1640 // get the reconstructor object and the loader for a detector
1641
1642   if (fReconstructor[iDet]) return fReconstructor[iDet];
1643
1644   // load the reconstructor object
1645   TPluginManager* pluginManager = gROOT->GetPluginManager();
1646   TString detName = fgkDetectorName[iDet];
1647   TString recName = "Ali" + detName + "Reconstructor";
1648   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1649
1650   AliReconstructor* reconstructor = NULL;
1651   // first check if a plugin is defined for the reconstructor
1652   TPluginHandler* pluginHandler = 
1653     pluginManager->FindHandler("AliReconstructor", detName);
1654   // if not, add a plugin for it
1655   if (!pluginHandler) {
1656     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1657     TString libs = gSystem->GetLibraries();
1658     if (libs.Contains("lib" + detName + "base.so") ||
1659         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1660       pluginManager->AddHandler("AliReconstructor", detName, 
1661                                 recName, detName + "rec", recName + "()");
1662     } else {
1663       pluginManager->AddHandler("AliReconstructor", detName, 
1664                                 recName, detName, recName + "()");
1665     }
1666     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1667   }
1668   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1669     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1670   }
1671   if (reconstructor) {
1672     TObject* obj = fOptions.FindObject(detName.Data());
1673     if (obj) reconstructor->SetOption(obj->GetTitle());
1674     reconstructor->Init();
1675     fReconstructor[iDet] = reconstructor;
1676   }
1677
1678   // get or create the loader
1679   if (detName != "HLT") {
1680     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1681     if (!fLoader[iDet]) {
1682       AliConfig::Instance()
1683         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
1684                                 detName, detName);
1685       // first check if a plugin is defined for the loader
1686       pluginHandler = 
1687         pluginManager->FindHandler("AliLoader", detName);
1688       // if not, add a plugin for it
1689       if (!pluginHandler) {
1690         TString loaderName = "Ali" + detName + "Loader";
1691         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1692         pluginManager->AddHandler("AliLoader", detName, 
1693                                   loaderName, detName + "base", 
1694                                   loaderName + "(const char*, TFolder*)");
1695         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1696       }
1697       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1698         fLoader[iDet] = 
1699           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1700                                                  fRunLoader->GetEventFolder());
1701       }
1702       if (!fLoader[iDet]) {   // use default loader
1703         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1704       }
1705       if (!fLoader[iDet]) {
1706         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1707         if (fStopOnError) return NULL;
1708       } else {
1709         fRunLoader->AddLoader(fLoader[iDet]);
1710         fRunLoader->CdGAFile();
1711         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1712         fRunLoader->Write(0, TObject::kOverwrite);
1713       }
1714     }
1715   }
1716       
1717   return reconstructor;
1718 }
1719
1720 //_____________________________________________________________________________
1721 Bool_t AliReconstruction::CreateVertexer()
1722 {
1723 // create the vertexer
1724
1725   fVertexer = NULL;
1726   AliReconstructor* itsReconstructor = GetReconstructor(0);
1727   if (itsReconstructor) {
1728     fVertexer = itsReconstructor->CreateVertexer();
1729   }
1730   if (!fVertexer) {
1731     AliWarning("couldn't create a vertexer for ITS");
1732     if (fStopOnError) return kFALSE;
1733   }
1734
1735   return kTRUE;
1736 }
1737
1738 //_____________________________________________________________________________
1739 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1740 {
1741 // create the trackers
1742
1743   TString detStr = detectors;
1744   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1745     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1746     AliReconstructor* reconstructor = GetReconstructor(iDet);
1747     if (!reconstructor) continue;
1748     TString detName = fgkDetectorName[iDet];
1749     if (detName == "HLT") {
1750       fRunHLTTracking = kTRUE;
1751       continue;
1752     }
1753     if (detName == "MUON") {
1754       fRunMuonTracking = kTRUE;
1755       continue;
1756     }
1757
1758
1759     fTracker[iDet] = reconstructor->CreateTracker();
1760     if (!fTracker[iDet] && (iDet < 7)) {
1761       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1762       if (fStopOnError) return kFALSE;
1763     }
1764   }
1765
1766   return kTRUE;
1767 }
1768
1769 //_____________________________________________________________________________
1770 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1771 {
1772 // delete trackers and the run loader and close and delete the file
1773
1774   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1775     delete fReconstructor[iDet];
1776     fReconstructor[iDet] = NULL;
1777     fLoader[iDet] = NULL;
1778     delete fTracker[iDet];
1779     fTracker[iDet] = NULL;
1780     delete fQualAssDataMaker[iDet];
1781     fQualAssDataMaker[iDet] = NULL;
1782   }
1783   delete fVertexer;
1784   fVertexer = NULL;
1785   delete fDiamondProfile;
1786   fDiamondProfile = NULL;
1787
1788   delete fRunLoader;
1789   fRunLoader = NULL;
1790   delete fRawReader;
1791   fRawReader = NULL;
1792
1793   if (file) {
1794     file->Close();
1795     delete file;
1796   }
1797
1798   if (fileOld) {
1799     fileOld->Close();
1800     delete fileOld;
1801     gSystem->Unlink("AliESDs.old.root");
1802   }
1803 }
1804
1805 //_____________________________________________________________________________
1806
1807 Bool_t AliReconstruction::ReadESD(AliESDEvent*& esd, const char* recStep) const
1808 {
1809 // read the ESD event from a file
1810
1811   if (!esd) return kFALSE;
1812   char fileName[256];
1813   sprintf(fileName, "ESD_%d.%d_%s.root", 
1814           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1815   if (gSystem->AccessPathName(fileName)) return kFALSE;
1816
1817   AliInfo(Form("reading ESD from file %s", fileName));
1818   AliDebug(1, Form("reading ESD from file %s", fileName));
1819   TFile* file = TFile::Open(fileName);
1820   if (!file || !file->IsOpen()) {
1821     AliError(Form("opening %s failed", fileName));
1822     delete file;
1823     return kFALSE;
1824   }
1825
1826   gROOT->cd();
1827   delete esd;
1828   esd = (AliESDEvent*) file->Get("ESD");
1829   file->Close();
1830   delete file;
1831   return kTRUE;
1832
1833 }
1834
1835
1836
1837 //_____________________________________________________________________________
1838 void AliReconstruction::WriteESD(AliESDEvent* esd, const char* recStep) const
1839 {
1840 // write the ESD event to a file
1841
1842   if (!esd) return;
1843   char fileName[256];
1844   sprintf(fileName, "ESD_%d.%d_%s.root", 
1845           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1846
1847   AliDebug(1, Form("writing ESD to file %s", fileName));
1848   TFile* file = TFile::Open(fileName, "recreate");
1849   if (!file || !file->IsOpen()) {
1850     AliError(Form("opening %s failed", fileName));
1851   } else {
1852     esd->Write("ESD");
1853     file->Close();
1854   }
1855   delete file;
1856 }
1857
1858
1859
1860
1861
1862 //_____________________________________________________________________________
1863 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1864 {
1865   // write all files from the given esd file to an aod file
1866
1867   // create an AliAOD object 
1868   AliAODEvent *aod = new AliAODEvent();
1869   aod->CreateStdContent();
1870   
1871   // go to the file
1872   aodFile->cd();
1873   
1874   // create the tree
1875   TTree *aodTree = new TTree("aodTree", "AliAOD tree");
1876   aodTree->Branch(aod->GetList());
1877
1878   // connect to ESD
1879   TTree *t = (TTree*) esdFile->Get("esdTree");
1880   AliESDEvent *esd = new AliESDEvent();
1881   esd->ReadFromTree(t);
1882
1883   Int_t nEvents = t->GetEntries();
1884
1885   // set arrays and pointers
1886   Float_t posF[3];
1887   Double_t pos[3];
1888   Double_t p[3];
1889   Double_t covVtx[6];
1890   Double_t covTr[21];
1891   Double_t pid[10];
1892
1893   // loop over events and fill them
1894   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1895     t->GetEntry(iEvent);
1896
1897     // Multiplicity information needed by the header (to be revised!)
1898     Int_t nTracks   = esd->GetNumberOfTracks();
1899     Int_t nPosTracks = 0;
1900     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) 
1901       if (esd->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
1902
1903     // Access the header
1904     AliAODHeader *header = aod->GetHeader();
1905
1906     // fill the header
1907     header->SetRunNumber       (esd->GetRunNumber()       );
1908     header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
1909     header->SetOrbitNumber     (esd->GetOrbitNumber()     );
1910     header->SetPeriodNumber    (esd->GetPeriodNumber()    );
1911     header->SetTriggerMask     (esd->GetTriggerMask()     ); 
1912     header->SetTriggerCluster  (esd->GetTriggerCluster()  );
1913     header->SetEventType       (esd->GetEventType()       );
1914     header->SetMagneticField   (esd->GetMagneticField()   );
1915     header->SetZDCN1Energy     (esd->GetZDCN1Energy()     );
1916     header->SetZDCP1Energy     (esd->GetZDCP1Energy()     );
1917     header->SetZDCN2Energy     (esd->GetZDCN2Energy()     );
1918     header->SetZDCP2Energy     (esd->GetZDCP2Energy()     );
1919     header->SetZDCEMEnergy     (esd->GetZDCEMEnergy()     );
1920     header->SetRefMultiplicity   (nTracks);
1921     header->SetRefMultiplicityPos(nPosTracks);
1922     header->SetRefMultiplicityNeg(nTracks - nPosTracks);
1923     header->SetMuonMagFieldScale(-999.); // FIXME
1924     header->SetCentrality(-999.);        // FIXME
1925
1926     Int_t nV0s      = esd->GetNumberOfV0s();
1927     Int_t nCascades = esd->GetNumberOfCascades();
1928     Int_t nKinks    = esd->GetNumberOfKinks();
1929     Int_t nVertices = nV0s + nCascades + nKinks;
1930     
1931     aod->ResetStd(nTracks, nVertices);
1932     AliAODTrack *aodTrack;
1933     
1934     // Array to take into account the tracks already added to the AOD
1935     Bool_t * usedTrack = NULL;
1936     if (nTracks>0) {
1937       usedTrack = new Bool_t[nTracks];
1938       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
1939     }
1940     // Array to take into account the V0s already added to the AOD
1941     Bool_t * usedV0 = NULL;
1942     if (nV0s>0) {
1943       usedV0 = new Bool_t[nV0s];
1944       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
1945     }
1946     // Array to take into account the kinks already added to the AOD
1947     Bool_t * usedKink = NULL;
1948     if (nKinks>0) {
1949       usedKink = new Bool_t[nKinks];
1950       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
1951     }
1952
1953     // Access to the AOD container of vertices
1954     TClonesArray &vertices = *(aod->GetVertices());
1955     Int_t jVertices=0;
1956
1957     // Access to the AOD container of tracks
1958     TClonesArray &tracks = *(aod->GetTracks());
1959     Int_t jTracks=0; 
1960   
1961     // Add primary vertex. The primary tracks will be defined
1962     // after the loops on the composite objects (V0, cascades, kinks)
1963     const AliESDVertex *vtx = esd->GetPrimaryVertex();
1964       
1965     vtx->GetXYZ(pos); // position
1966     vtx->GetCovMatrix(covVtx); //covariance matrix
1967
1968     AliAODVertex * primary = new(vertices[jVertices++])
1969       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1970          
1971     // Create vertices starting from the most complex objects
1972       
1973     // Cascades
1974     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
1975       AliESDcascade *cascade = esd->GetCascade(nCascade);
1976       
1977       cascade->GetXYZcascade(pos[0], pos[1], pos[2]);  // Bo: bug correction
1978       cascade->GetPosCovXi(covVtx);
1979      
1980       // Add the cascade vertex
1981       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
1982                                                                         covVtx,
1983                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
1984                                                                         primary,
1985                                                                         nCascade,
1986                                                                         AliAODVertex::kCascade);
1987
1988       primary->AddDaughter(vcascade);
1989
1990       // Add the V0 from the cascade. The ESD class have to be optimized...
1991       // Now we have to search for the corresponding V0 in the list of V0s
1992       // using the indeces of the positive and negative tracks
1993
1994       Int_t posFromV0 = cascade->GetPindex();
1995       Int_t negFromV0 = cascade->GetNindex();
1996
1997       AliESDv0 * v0 = 0x0;
1998       Int_t indV0 = -1;
1999
2000       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2001
2002         v0 = esd->GetV0(iV0);
2003         Int_t posV0 = v0->GetPindex();
2004         Int_t negV0 = v0->GetNindex();
2005
2006         if (posV0==posFromV0 && negV0==negFromV0) {
2007           indV0 = iV0;
2008           break;
2009         }
2010       }
2011
2012       AliAODVertex * vV0FromCascade = 0x0;
2013
2014       if (indV0>-1 && !usedV0[indV0] ) {
2015         
2016         // the V0 exists in the array of V0s and is not used
2017
2018         usedV0[indV0] = kTRUE;
2019         
2020         v0->GetXYZ(pos[0], pos[1], pos[2]);
2021         v0->GetPosCov(covVtx);
2022         
2023         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2024                                                                  covVtx,
2025                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2026                                                                  vcascade,
2027                                                                  indV0,
2028                                                                  AliAODVertex::kV0);
2029       } else {
2030
2031         // the V0 doesn't exist in the array of V0s or was used
2032         cerr << "Error: event " << iEvent << " cascade " << nCascade
2033              << " The V0 " << indV0 
2034              << " doesn't exist in the array of V0s or was used!" << endl;
2035
2036         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2037         cascade->GetPosCov(covVtx);
2038       
2039         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2040                                                                  covVtx,
2041                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2042                                                                  vcascade,
2043                                                                  indV0,
2044                                                                  AliAODVertex::kV0);
2045         vcascade->AddDaughter(vV0FromCascade);
2046       }
2047
2048       // Add the positive tracks from the V0
2049
2050       if (! usedTrack[posFromV0]) {
2051
2052         usedTrack[posFromV0] = kTRUE;
2053
2054         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2055         esdTrack->GetPxPyPz(p);
2056         esdTrack->GetXYZ(pos);
2057         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2058         esdTrack->GetESDpid(pid);
2059         
2060         vV0FromCascade->AddDaughter(aodTrack =
2061                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2062                                            esdTrack->GetLabel(), 
2063                                            p, 
2064                                            kTRUE,
2065                                            pos,
2066                                            kFALSE,
2067                                            covTr, 
2068                                            (Short_t)esdTrack->Charge(),
2069                                            esdTrack->GetITSClusterMap(), 
2070                                            pid,
2071                                            vV0FromCascade,
2072                                            kTRUE,  // check if this is right
2073                                            kFALSE, // check if this is right
2074                                            AliAODTrack::kSecondary)
2075                 );
2076         aodTrack->ConvertAliPIDtoAODPID();
2077       }
2078       else {
2079         cerr << "Error: event " << iEvent << " cascade " << nCascade
2080              << " track " << posFromV0 << " has already been used!" << endl;
2081       }
2082
2083       // Add the negative tracks from the V0
2084
2085       if (!usedTrack[negFromV0]) {
2086         
2087         usedTrack[negFromV0] = kTRUE;
2088         
2089         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2090         esdTrack->GetPxPyPz(p);
2091         esdTrack->GetXYZ(pos);
2092         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2093         esdTrack->GetESDpid(pid);
2094         
2095         vV0FromCascade->AddDaughter(aodTrack =
2096                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2097                                            esdTrack->GetLabel(),
2098                                            p,
2099                                            kTRUE,
2100                                            pos,
2101                                            kFALSE,
2102                                            covTr, 
2103                                            (Short_t)esdTrack->Charge(),
2104                                            esdTrack->GetITSClusterMap(), 
2105                                            pid,
2106                                            vV0FromCascade,
2107                                            kTRUE,  // check if this is right
2108                                            kFALSE, // check if this is right
2109                                            AliAODTrack::kSecondary)
2110                 );
2111         aodTrack->ConvertAliPIDtoAODPID();
2112       }
2113       else {
2114         cerr << "Error: event " << iEvent << " cascade " << nCascade
2115              << " track " << negFromV0 << " has already been used!" << endl;
2116       }
2117
2118       // Add the bachelor track from the cascade
2119
2120       Int_t bachelor = cascade->GetBindex();
2121       
2122       if(!usedTrack[bachelor]) {
2123       
2124         usedTrack[bachelor] = kTRUE;
2125         
2126         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2127         esdTrack->GetPxPyPz(p);
2128         esdTrack->GetXYZ(pos);
2129         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2130         esdTrack->GetESDpid(pid);
2131
2132         vcascade->AddDaughter(aodTrack =
2133                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2134                                            esdTrack->GetLabel(),
2135                                            p,
2136                                            kTRUE,
2137                                            pos,
2138                                            kFALSE,
2139                                            covTr, 
2140                                            (Short_t)esdTrack->Charge(),
2141                                            esdTrack->GetITSClusterMap(), 
2142                                            pid,
2143                                            vcascade,
2144                                            kTRUE,  // check if this is right
2145                                            kFALSE, // check if this is right
2146                                            AliAODTrack::kSecondary)
2147                 );
2148         aodTrack->ConvertAliPIDtoAODPID();
2149      }
2150       else {
2151         cerr << "Error: event " << iEvent << " cascade " << nCascade
2152              << " track " << bachelor << " has already been used!" << endl;
2153       }
2154
2155       // Add the primary track of the cascade (if any)
2156
2157     } // end of the loop on cascades
2158     
2159     // V0s
2160         
2161     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2162
2163       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2164
2165       AliESDv0 *v0 = esd->GetV0(nV0);
2166       
2167       v0->GetXYZ(pos[0], pos[1], pos[2]);
2168       v0->GetPosCov(covVtx);
2169
2170       AliAODVertex * vV0 = 
2171         new(vertices[jVertices++]) AliAODVertex(pos,
2172                                                 covVtx,
2173                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2174                                                 primary,
2175                                                 nV0,
2176                                                 AliAODVertex::kV0);
2177       primary->AddDaughter(vV0);
2178
2179       Int_t posFromV0 = v0->GetPindex();
2180       Int_t negFromV0 = v0->GetNindex();
2181       
2182       // Add the positive tracks from the V0
2183
2184       if (!usedTrack[posFromV0]) {
2185         
2186         usedTrack[posFromV0] = kTRUE;
2187
2188         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2189         esdTrack->GetPxPyPz(p);
2190         esdTrack->GetXYZ(pos);
2191         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2192         esdTrack->GetESDpid(pid);
2193         
2194         vV0->AddDaughter(aodTrack =
2195                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2196                                            esdTrack->GetLabel(), 
2197                                            p, 
2198                                            kTRUE,
2199                                            pos,
2200                                            kFALSE,
2201                                            covTr, 
2202                                            (Short_t)esdTrack->Charge(),
2203                                            esdTrack->GetITSClusterMap(), 
2204                                            pid,
2205                                            vV0,
2206                                            kTRUE,  // check if this is right
2207                                            kFALSE, // check if this is right
2208                                            AliAODTrack::kSecondary)
2209                 );
2210         aodTrack->ConvertAliPIDtoAODPID();
2211       }
2212       else {
2213         cerr << "Error: event " << iEvent << " V0 " << nV0
2214              << " track " << posFromV0 << " has already been used!" << endl;
2215       }
2216
2217       // Add the negative tracks from the V0
2218
2219       if (!usedTrack[negFromV0]) {
2220
2221         usedTrack[negFromV0] = kTRUE;
2222
2223         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2224         esdTrack->GetPxPyPz(p);
2225         esdTrack->GetXYZ(pos);
2226         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2227         esdTrack->GetESDpid(pid);
2228
2229         vV0->AddDaughter(aodTrack =
2230                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2231                                            esdTrack->GetLabel(),
2232                                            p,
2233                                            kTRUE,
2234                                            pos,
2235                                            kFALSE,
2236                                            covTr, 
2237                                            (Short_t)esdTrack->Charge(),
2238                                            esdTrack->GetITSClusterMap(), 
2239                                            pid,
2240                                            vV0,
2241                                            kTRUE,  // check if this is right
2242                                            kFALSE, // check if this is right
2243                                            AliAODTrack::kSecondary)
2244                 );
2245         aodTrack->ConvertAliPIDtoAODPID();
2246       }
2247       else {
2248         cerr << "Error: event " << iEvent << " V0 " << nV0
2249              << " track " << negFromV0 << " has already been used!" << endl;
2250       }
2251
2252     } // end of the loop on V0s
2253     
2254     // Kinks: it is a big mess the access to the information in the kinks
2255     // The loop is on the tracks in order to find the mother and daugther of each kink
2256
2257
2258     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2259
2260
2261       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2262
2263       Int_t ikink = esdTrack->GetKinkIndex(0);
2264
2265       if (ikink) {
2266         // Negative kink index: mother, positive: daughter
2267
2268         // Search for the second track of the kink
2269
2270         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2271
2272           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2273
2274           Int_t jkink = esdTrack1->GetKinkIndex(0);
2275
2276           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2277
2278             // The two tracks are from the same kink
2279           
2280             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2281
2282             Int_t imother = -1;
2283             Int_t idaughter = -1;
2284
2285             if (ikink<0 && jkink>0) {
2286
2287               imother = iTrack;
2288               idaughter = jTrack;
2289             }
2290             else if (ikink>0 && jkink<0) {
2291
2292               imother = jTrack;
2293               idaughter = iTrack;
2294             }
2295             else {
2296               cerr << "Error: Wrong combination of kink indexes: "
2297               << ikink << " " << jkink << endl;
2298               continue;
2299             }
2300
2301             // Add the mother track
2302
2303             AliAODTrack * mother = NULL;
2304
2305             if (!usedTrack[imother]) {
2306         
2307               usedTrack[imother] = kTRUE;
2308         
2309               AliESDtrack *esdTrack = esd->GetTrack(imother);
2310               esdTrack->GetPxPyPz(p);
2311               esdTrack->GetXYZ(pos);
2312               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2313               esdTrack->GetESDpid(pid);
2314
2315               mother = 
2316                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2317                                            esdTrack->GetLabel(),
2318                                            p,
2319                                            kTRUE,
2320                                            pos,
2321                                            kFALSE,
2322                                            covTr, 
2323                                            (Short_t)esdTrack->Charge(),
2324                                            esdTrack->GetITSClusterMap(), 
2325                                            pid,
2326                                            primary,
2327                                            kTRUE, // check if this is right
2328                                            kTRUE, // check if this is right
2329                                            AliAODTrack::kPrimary);
2330               primary->AddDaughter(mother);
2331               mother->ConvertAliPIDtoAODPID();
2332             }
2333             else {
2334               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2335               << " track " << imother << " has already been used!" << endl;
2336             }
2337
2338             // Add the kink vertex
2339             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2340
2341             AliAODVertex * vkink = 
2342             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2343                                                     NULL,
2344                                                     0.,
2345                                                     mother,
2346                                                     esdTrack->GetID(), // This is the track ID of the mother's track!
2347                                                     AliAODVertex::kKink);
2348             // Add the daughter track
2349
2350             AliAODTrack * daughter = NULL;
2351
2352             if (!usedTrack[idaughter]) {
2353         
2354               usedTrack[idaughter] = kTRUE;
2355         
2356               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2357               esdTrack->GetPxPyPz(p);
2358               esdTrack->GetXYZ(pos);
2359               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2360               esdTrack->GetESDpid(pid);
2361
2362               daughter = 
2363                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2364                                            esdTrack->GetLabel(),
2365                                            p,
2366                                            kTRUE,
2367                                            pos,
2368                                            kFALSE,
2369                                            covTr, 
2370                                            (Short_t)esdTrack->Charge(),
2371                                            esdTrack->GetITSClusterMap(), 
2372                                            pid,
2373                                            vkink,
2374                                            kTRUE, // check if this is right
2375                                            kTRUE, // check if this is right
2376                                            AliAODTrack::kPrimary);
2377               vkink->AddDaughter(daughter);
2378               daughter->ConvertAliPIDtoAODPID();
2379             }
2380             else {
2381               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2382               << " track " << idaughter << " has already been used!" << endl;
2383             }
2384
2385
2386           }
2387         }
2388
2389       }      
2390
2391     }
2392
2393     
2394     // Tracks (primary and orphan)
2395       
2396     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2397         
2398
2399       if (usedTrack[nTrack]) continue;
2400
2401       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2402       esdTrack->GetPxPyPz(p);
2403       esdTrack->GetXYZ(pos);
2404       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2405       esdTrack->GetESDpid(pid);
2406
2407       Float_t impactXY, impactZ;
2408
2409       esdTrack->GetImpactParameters(impactXY,impactZ);
2410
2411       if (impactXY<3) {
2412         // track inside the beam pipe
2413       
2414         primary->AddDaughter(aodTrack =
2415             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2416                                          esdTrack->GetLabel(),
2417                                          p,
2418                                          kTRUE,
2419                                          pos,
2420                                          kFALSE,
2421                                          covTr, 
2422                                          (Short_t)esdTrack->Charge(),
2423                                          esdTrack->GetITSClusterMap(), 
2424                                          pid,
2425                                          primary,
2426                                          kTRUE, // check if this is right
2427                                          kTRUE, // check if this is right
2428                                          AliAODTrack::kPrimary)
2429             );
2430         aodTrack->ConvertAliPIDtoAODPID();
2431       }
2432       else {
2433         // outside the beam pipe: orphan track
2434             aodTrack =
2435             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2436                                          esdTrack->GetLabel(),
2437                                          p,
2438                                          kTRUE,
2439                                          pos,
2440                                          kFALSE,
2441                                          covTr, 
2442                                          (Short_t)esdTrack->Charge(),
2443                                          esdTrack->GetITSClusterMap(), 
2444                                          pid,
2445                                          NULL,
2446                                          kFALSE, // check if this is right
2447                                          kFALSE, // check if this is right
2448                                          AliAODTrack::kOrphan);
2449             aodTrack->ConvertAliPIDtoAODPID();
2450       } 
2451     } // end of loop on tracks
2452
2453     // muon tracks
2454     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2455     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2456       
2457       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);     
2458       p[0] = esdMuTrack->Px(); 
2459       p[1] = esdMuTrack->Py(); 
2460       p[2] = esdMuTrack->Pz();
2461       pos[0] = primary->GetX(); 
2462       pos[1] = primary->GetY(); 
2463       pos[2] = primary->GetZ();
2464       
2465       // has to be changed once the muon pid is provided by the ESD
2466       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2467       
2468       primary->AddDaughter(aodTrack =
2469           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2470                                              0, // no label provided
2471                                              p,
2472                                              kTRUE,
2473                                              pos,
2474                                              kFALSE,
2475                                              NULL, // no covariance matrix provided
2476                                              esdMuTrack->Charge(),
2477                                              0, // ITSClusterMap is set below
2478                                              pid,
2479                                              primary,
2480                                              kFALSE,  // muon tracks are not used to fit the primary vtx
2481                                              kFALSE,  // not used for vertex fit
2482                                              AliAODTrack::kPrimary)
2483           );
2484
2485       aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
2486       Int_t track2Trigger = esdMuTrack->GetMatchTrigger();
2487       aodTrack->SetMatchTrigger(track2Trigger);
2488       if (track2Trigger) 
2489         aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
2490       else 
2491         aodTrack->SetChi2MatchTrigger(0.);
2492     }
2493     
2494     // Access to the AOD container of clusters
2495     TClonesArray &clusters = *(aod->GetClusters());
2496     Int_t jClusters=0;
2497
2498     // Calo Clusters
2499     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2500
2501     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2502
2503       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2504
2505       Int_t id = cluster->GetID();
2506       Int_t label = -1;
2507       Float_t energy = cluster->E();
2508       cluster->GetPosition(posF);
2509       AliAODVertex *prodVertex = primary;
2510       AliAODTrack *primTrack = NULL;
2511       Char_t ttype=AliAODCluster::kUndef;
2512
2513       if (cluster->GetClusterType() == AliESDCaloCluster::kPHOSCluster){
2514         ttype = AliAODCluster::kPHOSNeutral;
2515       }
2516       else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1){
2517         ttype = AliAODCluster::kEMCALClusterv1;
2518       }      
2519       else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALPseudoCluster){
2520         ttype = AliAODCluster::kEMCALPseudoCluster;
2521       }    
2522
2523       new(clusters[jClusters++]) AliAODCluster(id,
2524                                                label,
2525                                                energy,
2526                                                pos,
2527                                                NULL, // no covariance matrix provided
2528                                                NULL, // no pid for clusters provided
2529                                                prodVertex,
2530                                                primTrack,
2531                                                ttype);
2532
2533     } // end of loop on calo clusters
2534
2535     // tracklets
2536     const AliMultiplicity *mult = esd->GetMultiplicity();
2537     if (mult) {
2538       if (mult->GetNumberOfTracklets()>0) {
2539         aod->GetTracklets()->CreateContainer(mult->GetNumberOfTracklets());
2540
2541         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
2542           aod->GetTracklets()->SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n));
2543         }
2544       }
2545     } else {
2546       Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
2547     }
2548
2549     delete [] usedTrack;
2550     delete [] usedV0;
2551     delete [] usedKink;
2552
2553     // fill the tree for this event
2554     aodTree->Fill();
2555   } // end of event loop
2556
2557   aodTree->GetUserInfo()->Add(aod);
2558
2559   // close ESD file
2560   esdFile->Close();
2561
2562   // write the tree to the specified file
2563   aodFile = aodTree->GetCurrentFile();
2564   aodFile->cd();
2565   aodTree->Write();
2566
2567   return;
2568 }
2569
2570 void AliReconstruction::WriteAlignmentData(AliESDEvent* esd)
2571 {
2572   // Write space-points which are then used in the alignment procedures
2573   // For the moment only ITS, TRD and TPC
2574
2575   // Load TOF clusters
2576   if (fTracker[3]){
2577     fLoader[3]->LoadRecPoints("read");
2578     TTree* tree = fLoader[3]->TreeR();
2579     if (!tree) {
2580       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2581       return;
2582     }
2583     fTracker[3]->LoadClusters(tree);
2584   }
2585   Int_t ntracks = esd->GetNumberOfTracks();
2586   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2587     {
2588       AliESDtrack *track = esd->GetTrack(itrack);
2589       Int_t nsp = 0;
2590       Int_t idx[200];
2591       for (Int_t iDet = 3; iDet >= 0; iDet--)
2592         nsp += track->GetNcls(iDet);
2593       if (nsp) {
2594         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2595         track->SetTrackPointArray(sp);
2596         Int_t isptrack = 0;
2597         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2598           AliTracker *tracker = fTracker[iDet];
2599           if (!tracker) continue;
2600           Int_t nspdet = track->GetNcls(iDet);
2601           if (nspdet <= 0) continue;
2602           track->GetClusters(iDet,idx);
2603           AliTrackPoint p;
2604           Int_t isp = 0;
2605           Int_t isp2 = 0;
2606           while (isp < nspdet) {
2607             Bool_t isvalid;
2608             if(IsSelected(fgkDetectorName[iDet],fUseTrackingErrorsForAlignment)) {
2609               isvalid = tracker->GetTrackPointTrackingError(idx[isp2],p,track);
2610             } else {
2611               isvalid = tracker->GetTrackPoint(idx[isp2],p); 
2612             } 
2613             isp2++;
2614             const Int_t kNTPCmax = 159;
2615             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2616             if (!isvalid) continue;
2617             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2618           }
2619         }       
2620       }
2621     }
2622   if (fTracker[3]){
2623     fTracker[3]->UnloadClusters();
2624     fLoader[3]->UnloadRecPoints();
2625   }
2626 }
2627
2628 //_____________________________________________________________________________
2629 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESDEvent* esd)
2630 {
2631   // The method reads the raw-data error log
2632   // accumulated within the rawReader.
2633   // It extracts the raw-data errors related to
2634   // the current event and stores them into
2635   // a TClonesArray inside the esd object.
2636
2637   if (!fRawReader) return;
2638
2639   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2640
2641     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2642     if (!log) continue;
2643     if (iEvent != log->GetEventNumber()) continue;
2644
2645     esd->AddRawDataErrorLog(log);
2646   }
2647
2648 }
2649
2650 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2651   // Dump a file content into a char in TNamed
2652   ifstream in;
2653   in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2654   Int_t kBytes = (Int_t)in.tellg();
2655   printf("Size: %d \n",kBytes);
2656   TNamed *fn = 0;
2657   if(in.good()){
2658     char* memblock = new char [kBytes];
2659     in.seekg (0, ios::beg);
2660     in.read (memblock, kBytes);
2661     in.close();
2662     TString fData(memblock,kBytes);
2663     fn = new TNamed(fName,fData);
2664     printf("fData Size: %d \n",fData.Sizeof());
2665     printf("fName Size: %d \n",fName.Sizeof());
2666     printf("fn    Size: %d \n",fn->Sizeof());
2667     delete[] memblock;
2668   }
2669   else{
2670     AliInfo(Form("Could not Open %s\n",fPath.Data()));
2671   }
2672
2673   return fn;
2674 }
2675
2676 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2677   // This is not really needed in AliReconstruction at the moment
2678   // but can serve as a template
2679
2680   TList *fList = fTree->GetUserInfo();
2681   TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2682   printf("fn Size: %d \n",fn->Sizeof());
2683
2684   TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2685   const char* cdata = fn->GetTitle();
2686   printf("fTmp Size %d\n",fTmp.Sizeof());
2687
2688   int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2689   printf("calculated size %d\n",size);
2690   ofstream out(fName.Data(),ios::out | ios::binary);
2691   out.write(cdata,size);
2692   out.close();
2693
2694 }
2695
2696 //_____________________________________________________________________________
2697 AliQualAssDataMaker * AliReconstruction::GetQualAssDataMaker(Int_t iDet)
2698 {
2699 // get the quality assurance data maker object and the loader for a detector
2700
2701   if (fQualAssDataMaker[iDet]) 
2702     return fQualAssDataMaker[iDet];
2703
2704   // load the QA data maker object
2705   TPluginManager* pluginManager = gROOT->GetPluginManager();
2706   TString detName = fgkDetectorName[iDet];
2707   TString qadmName = "Ali" + detName + "QualAssDataMaker";
2708   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) 
2709     return NULL;
2710
2711   AliQualAssDataMaker * qadm = NULL;
2712   // first check if a plugin is defined for the quality assurance data maker
2713   TPluginHandler* pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2714   // if not, add a plugin for it
2715   if (!pluginHandler) {
2716     AliDebug(1, Form("defining plugin for %s", qadmName.Data()));
2717     TString libs = gSystem->GetLibraries();
2718     if (libs.Contains("lib" + detName + "base.so") ||
2719         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
2720       pluginManager->AddHandler("AliQualAssDataMaker", detName, 
2721                                 qadmName, detName + "qadm", qadmName + "()");
2722     } else {
2723       pluginManager->AddHandler("AliQualAssDataMaker", detName, 
2724                                 qadmName, detName, qadmName + "()");
2725     }
2726     pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2727   }
2728   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
2729     qadm = (AliQualAssDataMaker *) pluginHandler->ExecPlugin(0);
2730   }
2731   if (qadm) {
2732     AliInfo(Form("Initializing quality assurance data maker for %s", fgkDetectorName[iDet]));
2733     qadm->Init(AliQualAss::kRECPOINTS, fRunLoader->GetHeader()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2734     qadm->Init(AliQualAss::kESDS, fRunLoader->GetHeader()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2735     qadm->StartOfCycle(AliQualAss::kRECPOINTS);
2736         qadm->StartOfCycle(AliQualAss::kESDS, "same") ;         
2737     fQualAssDataMaker[iDet] = qadm;
2738   }
2739
2740   return qadm;
2741 }
2742
2743 //_____________________________________________________________________________
2744 Bool_t AliReconstruction::RunQualAss(const char* detectors, AliESDEvent *& esd)
2745 {
2746   // run the Quality Assurance data producer
2747
2748   AliCodeTimerAuto("")
2749   TString detStr = detectors;
2750   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
2751    if (!IsSelected(fgkDetectorName[iDet], detStr)) 
2752      continue;
2753    AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
2754    if (!qadm) 
2755      continue;
2756    AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2757    AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2758     
2759       qadm->Exec(AliQualAss::kESDS, esd) ; 
2760       qadm->Increment() ; 
2761
2762    AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2763  }
2764  if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
2765    AliError(Form("the following detectors were not found: %s",
2766                  detStr.Data()));
2767    if (fStopOnError) 
2768      return kFALSE;
2769  }
2770  
2771  return kTRUE;
2772   
2773
2774 }
2775
2776 //_____________________________________________________________________________
2777 Int_t AliReconstruction::GetDetIndex(const char* detector)
2778 {
2779   // return the detector index corresponding to detector
2780   Int_t index = -1 ; 
2781   for (index = 0; index < fgkNDetectors ; index++) {
2782     if ( strcmp(detector, fgkDetectorName[index]) == 0 )
2783         break ; 
2784   }     
2785   return index ; 
2786 }