166879342ed6260155bb3825b2ce4035d3752eb4
[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     qadm->EndOfCycle(AliQualAss::kRECPOINTS);
882     qadm->EndOfCycle(AliQualAss::kESDS);
883     qadm->Finish(AliQualAss::kRECPOINTS);
884     qadm->Finish(AliQualAss::kESDS) ; 
885   }
886
887   tree->GetUserInfo()->Add(esd);
888   hlttree->GetUserInfo()->Add(hltesd);
889
890
891
892   if(fESDPar.Contains("ESD.par")){
893     AliInfo("Attaching ESD.par to Tree");
894     TNamed *fn = CopyFileToTNamed(fESDPar.Data(),"ESD.par");
895     tree->GetUserInfo()->Add(fn);
896   }
897
898
899   file->cd();
900   if (fWriteESDfriend)
901     tree->SetBranchStatus("ESDfriend*",0);
902   tree->Write();
903   hlttree->Write();
904
905   if (fWriteAOD) {
906     TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
907     ESDFile2AODFile(file, aodFile);
908     aodFile->Close();
909   }
910
911   gROOT->cd();
912   CleanUp(file, fileOld);
913   
914   // Create tags for the events in the ESD tree (the ESD tree is always present)
915   // In case of empty events the tags will contain dummy values
916   AliESDTagCreator *esdtagCreator = new AliESDTagCreator();
917   esdtagCreator->CreateESDTags(fFirstEvent,fLastEvent);
918   if (fWriteAOD) {
919     AliAODTagCreator *aodtagCreator = new AliAODTagCreator();
920     aodtagCreator->CreateAODTags(fFirstEvent,fLastEvent);
921   }
922
923   return kTRUE;
924 }
925
926
927 //_____________________________________________________________________________
928 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
929 {
930 // run the local reconstruction
931   static Int_t eventNr=0;
932   AliCodeTimerAuto("")
933
934  //  AliCDBManager* man = AliCDBManager::Instance();
935 //   Bool_t origCache = man->GetCacheFlag();
936
937 //   TString detStr = detectors;
938 //   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
939 //     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
940 //     AliReconstructor* reconstructor = GetReconstructor(iDet);
941 //     if (!reconstructor) continue;
942 //     if (reconstructor->HasLocalReconstruction()) continue;
943
944 //     AliCodeTimerStart(Form("running reconstruction for %s", fgkDetectorName[iDet]));
945 //     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
946     
947 //     AliCodeTimerStart(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));                          
948 //     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
949
950 //     man->SetCacheFlag(kTRUE);
951 //     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
952 //     man->GetAll(calibPath); // entries are cached!
953
954 //     AliCodeTimerStop(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
955      
956 //     if (fRawReader) {
957 //       fRawReader->RewindEvents();
958 //       reconstructor->Reconstruct(fRunLoader, fRawReader);
959 //     } else {
960 //       reconstructor->Reconstruct(fRunLoader);
961 //     }
962      
963 //      AliCodeTimerStop(Form("running reconstruction for %s", fgkDetectorName[iDet]));
964     // AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr));
965
966 //     // unload calibration data
967 //     man->UnloadFromCache(calibPath);
968 //     //man->ClearCache();
969 //   }
970
971 //   man->SetCacheFlag(origCache);
972
973 //   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
974 //     AliError(Form("the following detectors were not found: %s",
975 //                   detStr.Data()));
976 //     if (fStopOnError) return kFALSE;
977 //   }
978
979           eventNr++;
980   return kTRUE;
981 }
982
983 //_____________________________________________________________________________
984 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
985 {
986 // run the local reconstruction
987   static Int_t eventNr=0;
988   AliCodeTimerAuto("")
989
990   TString detStr = detectors;
991   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
992     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
993     AliReconstructor* reconstructor = GetReconstructor(iDet);
994     if (!reconstructor) continue;
995     AliLoader* loader = fLoader[iDet];
996     if (!loader) {
997       AliWarning(Form("No loader is defined for %s!",fgkDetectorName[iDet]));
998       continue;
999     }
1000
1001     // conversion of digits
1002     if (fRawReader && reconstructor->HasDigitConversion()) {
1003       AliInfo(Form("converting raw data digits into root objects for %s", 
1004                    fgkDetectorName[iDet]));
1005       AliCodeTimerAuto(Form("converting raw data digits into root objects for %s", 
1006                             fgkDetectorName[iDet]));
1007       loader->LoadDigits("update");
1008       loader->CleanDigits();
1009       loader->MakeDigitsContainer();
1010       TTree* digitsTree = loader->TreeD();
1011       reconstructor->ConvertDigits(fRawReader, digitsTree);
1012       loader->WriteDigits("OVERWRITE");
1013       loader->UnloadDigits();
1014     }
1015
1016     // local reconstruction
1017     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1018     AliCodeTimerAuto(Form("running reconstruction for %s", fgkDetectorName[iDet]));
1019     loader->LoadRecPoints("update");
1020     loader->CleanRecPoints();
1021     loader->MakeRecPointsContainer();
1022     TTree* clustersTree = loader->TreeR();
1023     if (fRawReader && !reconstructor->HasDigitConversion()) {
1024       reconstructor->Reconstruct(fRawReader, clustersTree);
1025     } else {
1026       loader->LoadDigits("read");
1027       TTree* digitsTree = loader->TreeD();
1028       if (!digitsTree) {
1029         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
1030         if (fStopOnError) return kFALSE;
1031       } else {
1032         reconstructor->Reconstruct(digitsTree, clustersTree);
1033       }
1034       loader->UnloadDigits();
1035     }
1036
1037     AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
1038     if (qadm) {
1039       AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1040       AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1041           
1042     if (qadm->IsCycleDone() ) {
1043       qadm->EndOfCycle(AliQualAss::kRECPOINTS) ; 
1044           qadm->EndOfCycle(AliQualAss::kESDS) ; 
1045       qadm->StartOfCycle(AliQualAss::kRECPOINTS) ; 
1046           qadm->StartOfCycle(AliQualAss::kESDS, "same") ; 
1047    }
1048       qadm->Exec(AliQualAss::kRECPOINTS, clustersTree) ; 
1049       AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
1050     }
1051
1052     loader->WriteRecPoints("OVERWRITE");
1053     loader->UnloadRecPoints();
1054      AliSysInfo::AddStamp(Form("LRec%s_%d",fgkDetectorName[iDet],eventNr));
1055   }
1056
1057   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1058     AliError(Form("the following detectors were not found: %s",
1059                   detStr.Data()));
1060     if (fStopOnError) return kFALSE;
1061   }
1062   eventNr++;
1063   return kTRUE;
1064 }
1065
1066 //_____________________________________________________________________________
1067 Bool_t AliReconstruction::RunVertexFinder(AliESDEvent*& esd)
1068 {
1069 // run the barrel tracking
1070
1071   AliCodeTimerAuto("")
1072
1073   AliESDVertex* vertex = NULL;
1074   Double_t vtxPos[3] = {0, 0, 0};
1075   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
1076   TArrayF mcVertex(3); 
1077   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
1078     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
1079     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
1080   }
1081
1082   if (fVertexer) {
1083     if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
1084     AliInfo("running the ITS vertex finder");
1085     if (fLoader[0]) fLoader[0]->LoadRecPoints();
1086     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
1087     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
1088     if(!vertex){
1089       AliWarning("Vertex not found");
1090       vertex = new AliESDVertex();
1091       vertex->SetName("default");
1092     }
1093     else {
1094       vertex->SetName("reconstructed");
1095     }
1096
1097   } else {
1098     AliInfo("getting the primary vertex from MC");
1099     vertex = new AliESDVertex(vtxPos, vtxErr);
1100   }
1101
1102   if (vertex) {
1103     vertex->GetXYZ(vtxPos);
1104     vertex->GetSigmaXYZ(vtxErr);
1105   } else {
1106     AliWarning("no vertex reconstructed");
1107     vertex = new AliESDVertex(vtxPos, vtxErr);
1108   }
1109   esd->SetVertex(vertex);
1110   // if SPD multiplicity has been determined, it is stored in the ESD
1111   AliMultiplicity *mult = fVertexer->GetMultiplicity();
1112   if(mult)esd->SetMultiplicity(mult);
1113
1114   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1115     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1116   }  
1117   delete vertex;
1118
1119   return kTRUE;
1120 }
1121
1122 //_____________________________________________________________________________
1123 Bool_t AliReconstruction::RunHLTTracking(AliESDEvent*& esd)
1124 {
1125 // run the HLT barrel tracking
1126
1127   AliCodeTimerAuto("")
1128
1129   if (!fRunLoader) {
1130     AliError("Missing runLoader!");
1131     return kFALSE;
1132   }
1133
1134   AliInfo("running HLT tracking");
1135
1136   // Get a pointer to the HLT reconstructor
1137   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1138   if (!reconstructor) return kFALSE;
1139
1140   // TPC + ITS
1141   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1142     TString detName = fgkDetectorName[iDet];
1143     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1144     reconstructor->SetOption(detName.Data());
1145     AliTracker *tracker = reconstructor->CreateTracker();
1146     if (!tracker) {
1147       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1148       if (fStopOnError) return kFALSE;
1149       continue;
1150     }
1151     Double_t vtxPos[3];
1152     Double_t vtxErr[3]={0.005,0.005,0.010};
1153     const AliESDVertex *vertex = esd->GetVertex();
1154     vertex->GetXYZ(vtxPos);
1155     tracker->SetVertex(vtxPos,vtxErr);
1156     if(iDet != 1) {
1157       fLoader[iDet]->LoadRecPoints("read");
1158       TTree* tree = fLoader[iDet]->TreeR();
1159       if (!tree) {
1160         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1161         return kFALSE;
1162       }
1163       tracker->LoadClusters(tree);
1164     }
1165     if (tracker->Clusters2Tracks(esd) != 0) {
1166       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1167       return kFALSE;
1168     }
1169     if(iDet != 1) {
1170       tracker->UnloadClusters();
1171     }
1172     delete tracker;
1173   }
1174
1175   return kTRUE;
1176 }
1177
1178 //_____________________________________________________________________________
1179 Bool_t AliReconstruction::RunMuonTracking(AliESDEvent*& esd)
1180 {
1181 // run the muon spectrometer tracking
1182
1183   AliCodeTimerAuto("")
1184
1185   if (!fRunLoader) {
1186     AliError("Missing runLoader!");
1187     return kFALSE;
1188   }
1189   Int_t iDet = 7; // for MUON
1190
1191   AliInfo("is running...");
1192
1193   // Get a pointer to the MUON reconstructor
1194   AliReconstructor *reconstructor = GetReconstructor(iDet);
1195   if (!reconstructor) return kFALSE;
1196
1197   
1198   TString detName = fgkDetectorName[iDet];
1199   AliDebug(1, Form("%s tracking", detName.Data()));
1200   AliTracker *tracker =  reconstructor->CreateTracker();
1201   if (!tracker) {
1202     AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1203     return kFALSE;
1204   }
1205      
1206   // create Tracks
1207   fLoader[iDet]->LoadTracks("update");
1208   fLoader[iDet]->CleanTracks();
1209   fLoader[iDet]->MakeTracksContainer();
1210
1211   // read RecPoints
1212   fLoader[iDet]->LoadRecPoints("read");  
1213   tracker->LoadClusters(fLoader[iDet]->TreeR());
1214   
1215   Int_t rv = tracker->Clusters2Tracks(esd);
1216   
1217   fLoader[iDet]->UnloadRecPoints();
1218   
1219   if ( rv )
1220   {
1221     AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1222     return kFALSE;
1223   }
1224   
1225   tracker->UnloadClusters();
1226   
1227   fLoader[iDet]->UnloadRecPoints();
1228
1229   fLoader[iDet]->WriteTracks("OVERWRITE");
1230   fLoader[iDet]->UnloadTracks();
1231
1232   delete tracker;
1233   
1234   return kTRUE;
1235 }
1236
1237
1238 //_____________________________________________________________________________
1239 Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd)
1240 {
1241 // run the barrel tracking
1242   static Int_t eventNr=0;
1243   AliCodeTimerAuto("")
1244
1245   AliInfo("running tracking");
1246
1247   //Fill the ESD with the T0 info (will be used by the TOF) 
1248   if (fReconstructor[11] && fLoader[11]) {
1249     fLoader[11]->LoadRecPoints("READ");
1250     TTree *treeR = fLoader[11]->TreeR();
1251     GetReconstructor(11)->FillESD((TTree *)NULL,treeR,esd);
1252   }
1253
1254   // pass 1: TPC + ITS inwards
1255   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1256     if (!fTracker[iDet]) continue;
1257     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1258
1259     // load clusters
1260     fLoader[iDet]->LoadRecPoints("read");
1261     TTree* tree = fLoader[iDet]->TreeR();
1262     if (!tree) {
1263       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1264       return kFALSE;
1265     }
1266     fTracker[iDet]->LoadClusters(tree);
1267
1268     // run tracking
1269     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1270       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1271       return kFALSE;
1272     }
1273     if (fCheckPointLevel > 1) {
1274       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1275     }
1276     // preliminary PID in TPC needed by the ITS tracker
1277     if (iDet == 1) {
1278       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1279       AliESDpid::MakePID(esd);
1280     } 
1281      AliSysInfo::AddStamp(Form("Tracking0%s_%d",fgkDetectorName[iDet],eventNr));
1282   }
1283
1284   // pass 2: ALL backwards
1285   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1286     if (!fTracker[iDet]) continue;
1287     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1288
1289     // load clusters
1290     if (iDet > 1) {     // all except ITS, TPC
1291       TTree* tree = NULL;
1292       fLoader[iDet]->LoadRecPoints("read");
1293       tree = fLoader[iDet]->TreeR();
1294       if (!tree) {
1295         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1296         return kFALSE;
1297       }
1298       fTracker[iDet]->LoadClusters(tree); 
1299        AliSysInfo::AddStamp(Form("LoadCluster0%s_%d",fgkDetectorName[iDet],eventNr));
1300     }
1301
1302     // run tracking
1303     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1304       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1305       //      return kFALSE;
1306     }
1307     if (fCheckPointLevel > 1) {
1308       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1309     }
1310
1311     // unload clusters
1312     if (iDet > 2) {     // all except ITS, TPC, TRD
1313       fTracker[iDet]->UnloadClusters();
1314       fLoader[iDet]->UnloadRecPoints();
1315     }
1316     // updated PID in TPC needed by the ITS tracker -MI
1317     if (iDet == 1) {
1318       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1319       AliESDpid::MakePID(esd);
1320     }
1321      AliSysInfo::AddStamp(Form("Tracking1%s_%d",fgkDetectorName[iDet],eventNr));
1322   }
1323
1324   // write space-points to the ESD in case alignment data output
1325   // is switched on
1326   if (fWriteAlignmentData)
1327     WriteAlignmentData(esd);
1328
1329   // pass 3: TRD + TPC + ITS refit inwards
1330   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1331     if (!fTracker[iDet]) continue;
1332     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1333
1334     // run tracking
1335     if (fTracker[iDet]->RefitInward(esd) != 0) {
1336       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1337       //      return kFALSE;
1338     }
1339     if (fCheckPointLevel > 1) {
1340       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1341     }
1342
1343     // unload clusters
1344     fTracker[iDet]->UnloadClusters();
1345     fLoader[iDet]->UnloadRecPoints();
1346      AliSysInfo::AddStamp(Form("Tracking2%s_%d",fgkDetectorName[iDet],eventNr));
1347   }
1348   //
1349   // Propagate track to the vertex - if not done by ITS
1350   //
1351   Int_t ntracks = esd->GetNumberOfTracks();
1352   for (Int_t itrack=0; itrack<ntracks; itrack++){
1353     const Double_t kRadius  = 3;   // beam pipe radius
1354     const Double_t kMaxStep = 5;   // max step
1355     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1356     Double_t       fieldZ   = AliTracker::GetBz();  //
1357     AliESDtrack * track = esd->GetTrack(itrack);
1358     if (!track) continue;
1359     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1360    AliTracker::PropagateTrackTo(track,kRadius,track->GetMass(),kMaxStep,kTRUE);
1361     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1362   }
1363   eventNr++;
1364   return kTRUE;
1365 }
1366
1367 //_____________________________________________________________________________
1368 Bool_t AliReconstruction::CleanESD(AliESDEvent *esd){
1369   //
1370   // Remove the data which are not needed for the physics analysis.
1371   //
1372
1373   AliInfo("Cleaning the ESD...");
1374   Int_t nTracks=esd->GetNumberOfTracks();
1375   AliInfo(Form("Number of ESD tracks before cleaning %d",nTracks));
1376
1377   Float_t cleanPars[]={fDmax,fZmax};
1378   Bool_t rc=esd->Clean(cleanPars);
1379
1380   nTracks=esd->GetNumberOfTracks();
1381   AliInfo(Form("Number of ESD tracks after cleaning %d",nTracks));
1382
1383   return rc;
1384 }
1385
1386 //_____________________________________________________________________________
1387 Bool_t AliReconstruction::FillESD(AliESDEvent*& esd, const TString& detectors)
1388 {
1389 // fill the event summary data
1390
1391   AliCodeTimerAuto("")
1392     static Int_t eventNr=0; 
1393   TString detStr = detectors;
1394   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1395     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1396     AliReconstructor* reconstructor = GetReconstructor(iDet);
1397     if (!reconstructor) continue;
1398
1399     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1400       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1401       TTree* clustersTree = NULL;
1402       if (fLoader[iDet]) {
1403         fLoader[iDet]->LoadRecPoints("read");
1404         clustersTree = fLoader[iDet]->TreeR();
1405         if (!clustersTree) {
1406           AliError(Form("Can't get the %s clusters tree", 
1407                         fgkDetectorName[iDet]));
1408           if (fStopOnError) return kFALSE;
1409         }
1410       }
1411       if (fRawReader && !reconstructor->HasDigitConversion()) {
1412         reconstructor->FillESD(fRawReader, clustersTree, esd);
1413       } else {
1414         TTree* digitsTree = NULL;
1415         if (fLoader[iDet]) {
1416           fLoader[iDet]->LoadDigits("read");
1417           digitsTree = fLoader[iDet]->TreeD();
1418           if (!digitsTree) {
1419             AliError(Form("Can't get the %s digits tree", 
1420                           fgkDetectorName[iDet]));
1421             if (fStopOnError) return kFALSE;
1422           }
1423         }
1424         reconstructor->FillESD(digitsTree, clustersTree, esd);
1425         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1426       }
1427       if (fLoader[iDet]) {
1428         fLoader[iDet]->UnloadRecPoints();
1429       }
1430
1431       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1432     }
1433   }
1434
1435   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1436     AliError(Form("the following detectors were not found: %s", 
1437                   detStr.Data()));
1438     if (fStopOnError) return kFALSE;
1439   }
1440    AliSysInfo::AddStamp(Form("FillESD%d",eventNr));
1441   eventNr++;
1442   return kTRUE;
1443 }
1444
1445 //_____________________________________________________________________________
1446 Bool_t AliReconstruction::FillTriggerESD(AliESDEvent*& esd)
1447 {
1448   // Reads the trigger decision which is
1449   // stored in Trigger.root file and fills
1450   // the corresponding esd entries
1451
1452   AliCodeTimerAuto("")
1453   
1454   AliInfo("Filling trigger information into the ESD");
1455
1456   if (fRawReader) {
1457     AliCTPRawStream input(fRawReader);
1458     if (!input.Next()) {
1459       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1460       return kFALSE;
1461     }
1462     esd->SetTriggerMask(input.GetClassMask());
1463     esd->SetTriggerCluster(input.GetClusterMask());
1464   }
1465   else {
1466     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1467     if (runloader) {
1468       if (!runloader->LoadTrigger()) {
1469         AliCentralTrigger *aCTP = runloader->GetTrigger();
1470         esd->SetTriggerMask(aCTP->GetClassMask());
1471         esd->SetTriggerCluster(aCTP->GetClusterMask());
1472       }
1473       else {
1474         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1475         return kFALSE;
1476       }
1477     }
1478     else {
1479       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1480       return kFALSE;
1481     }
1482   }
1483
1484   return kTRUE;
1485 }
1486
1487
1488
1489
1490
1491 //_____________________________________________________________________________
1492 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESDEvent*& esd)
1493 {
1494   // 
1495   // Filling information from RawReader Header
1496   // 
1497
1498   AliInfo("Filling information from RawReader Header");
1499   esd->SetBunchCrossNumber(0);
1500   esd->SetOrbitNumber(0);
1501   esd->SetPeriodNumber(0);
1502   esd->SetTimeStamp(0);
1503   esd->SetEventType(0);
1504   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1505   if (eventHeader){
1506
1507     const UInt_t *id = eventHeader->GetP("Id");
1508     esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1509     esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1510     esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1511
1512     esd->SetTimeStamp((eventHeader->Get("Timestamp")));  
1513     esd->SetEventType((eventHeader->Get("Type")));
1514   }
1515
1516   return kTRUE;
1517 }
1518
1519
1520 //_____________________________________________________________________________
1521 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1522 {
1523 // check whether detName is contained in detectors
1524 // if yes, it is removed from detectors
1525
1526   // check if all detectors are selected
1527   if ((detectors.CompareTo("ALL") == 0) ||
1528       detectors.BeginsWith("ALL ") ||
1529       detectors.EndsWith(" ALL") ||
1530       detectors.Contains(" ALL ")) {
1531     detectors = "ALL";
1532     return kTRUE;
1533   }
1534
1535   // search for the given detector
1536   Bool_t result = kFALSE;
1537   if ((detectors.CompareTo(detName) == 0) ||
1538       detectors.BeginsWith(detName+" ") ||
1539       detectors.EndsWith(" "+detName) ||
1540       detectors.Contains(" "+detName+" ")) {
1541     detectors.ReplaceAll(detName, "");
1542     result = kTRUE;
1543   }
1544
1545   // clean up the detectors string
1546   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1547   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1548   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1549
1550   return result;
1551 }
1552
1553 //_____________________________________________________________________________
1554 Bool_t AliReconstruction::InitRunLoader()
1555 {
1556 // get or create the run loader
1557
1558   if (gAlice) delete gAlice;
1559   gAlice = NULL;
1560
1561   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1562     // load all base libraries to get the loader classes
1563     TString libs = gSystem->GetLibraries();
1564     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1565       TString detName = fgkDetectorName[iDet];
1566       if (detName == "HLT") continue;
1567       if (libs.Contains("lib" + detName + "base.so")) continue;
1568       gSystem->Load("lib" + detName + "base.so");
1569     }
1570     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1571     if (!fRunLoader) {
1572       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1573       CleanUp();
1574       return kFALSE;
1575     }
1576     fRunLoader->CdGAFile();
1577     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1578       if (fRunLoader->LoadgAlice() == 0) {
1579         gAlice = fRunLoader->GetAliRun();
1580         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1581       }
1582     }
1583     if (!gAlice && !fRawReader) {
1584       AliError(Form("no gAlice object found in file %s",
1585                     fGAliceFileName.Data()));
1586       CleanUp();
1587       return kFALSE;
1588     }
1589
1590     //PH This is a temporary fix to give access to the kinematics
1591     //PH that is needed for the labels of ITS clusters
1592     fRunLoader->LoadKinematics();
1593
1594   } else {               // galice.root does not exist
1595     if (!fRawReader) {
1596       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1597       CleanUp();
1598       return kFALSE;
1599     }
1600     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1601                                     AliConfig::GetDefaultEventFolderName(),
1602                                     "recreate");
1603     if (!fRunLoader) {
1604       AliError(Form("could not create run loader in file %s", 
1605                     fGAliceFileName.Data()));
1606       CleanUp();
1607       return kFALSE;
1608     }
1609     fRunLoader->MakeTree("E");
1610     Int_t iEvent = 0;
1611     while (fRawReader->NextEvent()) {
1612       fRunLoader->SetEventNumber(iEvent);
1613       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
1614                                      iEvent, iEvent);
1615       fRunLoader->MakeTree("H");
1616       fRunLoader->TreeE()->Fill();
1617       iEvent++;
1618     }
1619     fRawReader->RewindEvents();
1620     if (fNumberOfEventsPerFile > 0)
1621       fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1622     else
1623       fRunLoader->SetNumberOfEventsPerFile(iEvent);
1624     fRunLoader->WriteHeader("OVERWRITE");
1625     fRunLoader->CdGAFile();
1626     fRunLoader->Write(0, TObject::kOverwrite);
1627 //    AliTracker::SetFieldMap(???);
1628   }
1629
1630   return kTRUE;
1631 }
1632
1633 //_____________________________________________________________________________
1634 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1635 {
1636 // get the reconstructor object and the loader for a detector
1637
1638   if (fReconstructor[iDet]) return fReconstructor[iDet];
1639
1640   // load the reconstructor object
1641   TPluginManager* pluginManager = gROOT->GetPluginManager();
1642   TString detName = fgkDetectorName[iDet];
1643   TString recName = "Ali" + detName + "Reconstructor";
1644   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1645
1646   AliReconstructor* reconstructor = NULL;
1647   // first check if a plugin is defined for the reconstructor
1648   TPluginHandler* pluginHandler = 
1649     pluginManager->FindHandler("AliReconstructor", detName);
1650   // if not, add a plugin for it
1651   if (!pluginHandler) {
1652     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1653     TString libs = gSystem->GetLibraries();
1654     if (libs.Contains("lib" + detName + "base.so") ||
1655         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1656       pluginManager->AddHandler("AliReconstructor", detName, 
1657                                 recName, detName + "rec", recName + "()");
1658     } else {
1659       pluginManager->AddHandler("AliReconstructor", detName, 
1660                                 recName, detName, recName + "()");
1661     }
1662     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1663   }
1664   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1665     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1666   }
1667   if (reconstructor) {
1668     TObject* obj = fOptions.FindObject(detName.Data());
1669     if (obj) reconstructor->SetOption(obj->GetTitle());
1670     reconstructor->Init();
1671     fReconstructor[iDet] = reconstructor;
1672   }
1673
1674   // get or create the loader
1675   if (detName != "HLT") {
1676     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1677     if (!fLoader[iDet]) {
1678       AliConfig::Instance()
1679         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
1680                                 detName, detName);
1681       // first check if a plugin is defined for the loader
1682       pluginHandler = 
1683         pluginManager->FindHandler("AliLoader", detName);
1684       // if not, add a plugin for it
1685       if (!pluginHandler) {
1686         TString loaderName = "Ali" + detName + "Loader";
1687         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1688         pluginManager->AddHandler("AliLoader", detName, 
1689                                   loaderName, detName + "base", 
1690                                   loaderName + "(const char*, TFolder*)");
1691         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1692       }
1693       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1694         fLoader[iDet] = 
1695           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1696                                                  fRunLoader->GetEventFolder());
1697       }
1698       if (!fLoader[iDet]) {   // use default loader
1699         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1700       }
1701       if (!fLoader[iDet]) {
1702         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1703         if (fStopOnError) return NULL;
1704       } else {
1705         fRunLoader->AddLoader(fLoader[iDet]);
1706         fRunLoader->CdGAFile();
1707         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1708         fRunLoader->Write(0, TObject::kOverwrite);
1709       }
1710     }
1711   }
1712       
1713   return reconstructor;
1714 }
1715
1716 //_____________________________________________________________________________
1717 Bool_t AliReconstruction::CreateVertexer()
1718 {
1719 // create the vertexer
1720
1721   fVertexer = NULL;
1722   AliReconstructor* itsReconstructor = GetReconstructor(0);
1723   if (itsReconstructor) {
1724     fVertexer = itsReconstructor->CreateVertexer();
1725   }
1726   if (!fVertexer) {
1727     AliWarning("couldn't create a vertexer for ITS");
1728     if (fStopOnError) return kFALSE;
1729   }
1730
1731   return kTRUE;
1732 }
1733
1734 //_____________________________________________________________________________
1735 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1736 {
1737 // create the trackers
1738
1739   TString detStr = detectors;
1740   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1741     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1742     AliReconstructor* reconstructor = GetReconstructor(iDet);
1743     if (!reconstructor) continue;
1744     TString detName = fgkDetectorName[iDet];
1745     if (detName == "HLT") {
1746       fRunHLTTracking = kTRUE;
1747       continue;
1748     }
1749     if (detName == "MUON") {
1750       fRunMuonTracking = kTRUE;
1751       continue;
1752     }
1753
1754
1755     fTracker[iDet] = reconstructor->CreateTracker();
1756     if (!fTracker[iDet] && (iDet < 7)) {
1757       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1758       if (fStopOnError) return kFALSE;
1759     }
1760   }
1761
1762   return kTRUE;
1763 }
1764
1765 //_____________________________________________________________________________
1766 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1767 {
1768 // delete trackers and the run loader and close and delete the file
1769
1770   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1771     delete fReconstructor[iDet];
1772     fReconstructor[iDet] = NULL;
1773     fLoader[iDet] = NULL;
1774     delete fTracker[iDet];
1775     fTracker[iDet] = NULL;
1776     delete fQualAssDataMaker[iDet];
1777     fQualAssDataMaker[iDet] = NULL;
1778   }
1779   delete fVertexer;
1780   fVertexer = NULL;
1781   delete fDiamondProfile;
1782   fDiamondProfile = NULL;
1783
1784   delete fRunLoader;
1785   fRunLoader = NULL;
1786   delete fRawReader;
1787   fRawReader = NULL;
1788
1789   if (file) {
1790     file->Close();
1791     delete file;
1792   }
1793
1794   if (fileOld) {
1795     fileOld->Close();
1796     delete fileOld;
1797     gSystem->Unlink("AliESDs.old.root");
1798   }
1799 }
1800
1801 //_____________________________________________________________________________
1802
1803 Bool_t AliReconstruction::ReadESD(AliESDEvent*& esd, const char* recStep) const
1804 {
1805 // read the ESD event from a file
1806
1807   if (!esd) return kFALSE;
1808   char fileName[256];
1809   sprintf(fileName, "ESD_%d.%d_%s.root", 
1810           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1811   if (gSystem->AccessPathName(fileName)) return kFALSE;
1812
1813   AliInfo(Form("reading ESD from file %s", fileName));
1814   AliDebug(1, Form("reading ESD from file %s", fileName));
1815   TFile* file = TFile::Open(fileName);
1816   if (!file || !file->IsOpen()) {
1817     AliError(Form("opening %s failed", fileName));
1818     delete file;
1819     return kFALSE;
1820   }
1821
1822   gROOT->cd();
1823   delete esd;
1824   esd = (AliESDEvent*) file->Get("ESD");
1825   file->Close();
1826   delete file;
1827   return kTRUE;
1828
1829 }
1830
1831
1832
1833 //_____________________________________________________________________________
1834 void AliReconstruction::WriteESD(AliESDEvent* esd, const char* recStep) const
1835 {
1836 // write the ESD event to a file
1837
1838   if (!esd) return;
1839   char fileName[256];
1840   sprintf(fileName, "ESD_%d.%d_%s.root", 
1841           esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1842
1843   AliDebug(1, Form("writing ESD to file %s", fileName));
1844   TFile* file = TFile::Open(fileName, "recreate");
1845   if (!file || !file->IsOpen()) {
1846     AliError(Form("opening %s failed", fileName));
1847   } else {
1848     esd->Write("ESD");
1849     file->Close();
1850   }
1851   delete file;
1852 }
1853
1854
1855
1856
1857
1858 //_____________________________________________________________________________
1859 void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
1860 {
1861   // write all files from the given esd file to an aod file
1862
1863   // create an AliAOD object 
1864   AliAODEvent *aod = new AliAODEvent();
1865   aod->CreateStdContent();
1866   
1867   // go to the file
1868   aodFile->cd();
1869   
1870   // create the tree
1871   TTree *aodTree = new TTree("aodTree", "AliAOD tree");
1872   aodTree->Branch(aod->GetList());
1873
1874   // connect to ESD
1875   TTree *t = (TTree*) esdFile->Get("esdTree");
1876   AliESDEvent *esd = new AliESDEvent();
1877   esd->ReadFromTree(t);
1878
1879   Int_t nEvents = t->GetEntries();
1880
1881   // set arrays and pointers
1882   Float_t posF[3];
1883   Double_t pos[3];
1884   Double_t p[3];
1885   Double_t covVtx[6];
1886   Double_t covTr[21];
1887   Double_t pid[10];
1888
1889   // loop over events and fill them
1890   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1891     t->GetEntry(iEvent);
1892
1893     // Multiplicity information needed by the header (to be revised!)
1894     Int_t nTracks   = esd->GetNumberOfTracks();
1895     Int_t nPosTracks = 0;
1896     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) 
1897       if (esd->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
1898
1899     // Access the header
1900     AliAODHeader *header = aod->GetHeader();
1901
1902     // fill the header
1903     header->SetRunNumber       (esd->GetRunNumber()       );
1904     header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
1905     header->SetOrbitNumber     (esd->GetOrbitNumber()     );
1906     header->SetPeriodNumber    (esd->GetPeriodNumber()    );
1907     header->SetTriggerMask     (esd->GetTriggerMask()     ); 
1908     header->SetTriggerCluster  (esd->GetTriggerCluster()  );
1909     header->SetEventType       (esd->GetEventType()       );
1910     header->SetMagneticField   (esd->GetMagneticField()   );
1911     header->SetZDCN1Energy     (esd->GetZDCN1Energy()     );
1912     header->SetZDCP1Energy     (esd->GetZDCP1Energy()     );
1913     header->SetZDCN2Energy     (esd->GetZDCN2Energy()     );
1914     header->SetZDCP2Energy     (esd->GetZDCP2Energy()     );
1915     header->SetZDCEMEnergy     (esd->GetZDCEMEnergy()     );
1916     header->SetRefMultiplicity   (nTracks);
1917     header->SetRefMultiplicityPos(nPosTracks);
1918     header->SetRefMultiplicityNeg(nTracks - nPosTracks);
1919     header->SetMuonMagFieldScale(-999.); // FIXME
1920     header->SetCentrality(-999.);        // FIXME
1921
1922     Int_t nV0s      = esd->GetNumberOfV0s();
1923     Int_t nCascades = esd->GetNumberOfCascades();
1924     Int_t nKinks    = esd->GetNumberOfKinks();
1925     Int_t nVertices = nV0s + nCascades + nKinks;
1926     
1927     aod->ResetStd(nTracks, nVertices);
1928     AliAODTrack *aodTrack;
1929     
1930     // Array to take into account the tracks already added to the AOD
1931     Bool_t * usedTrack = NULL;
1932     if (nTracks>0) {
1933       usedTrack = new Bool_t[nTracks];
1934       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
1935     }
1936     // Array to take into account the V0s already added to the AOD
1937     Bool_t * usedV0 = NULL;
1938     if (nV0s>0) {
1939       usedV0 = new Bool_t[nV0s];
1940       for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
1941     }
1942     // Array to take into account the kinks already added to the AOD
1943     Bool_t * usedKink = NULL;
1944     if (nKinks>0) {
1945       usedKink = new Bool_t[nKinks];
1946       for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
1947     }
1948
1949     // Access to the AOD container of vertices
1950     TClonesArray &vertices = *(aod->GetVertices());
1951     Int_t jVertices=0;
1952
1953     // Access to the AOD container of tracks
1954     TClonesArray &tracks = *(aod->GetTracks());
1955     Int_t jTracks=0; 
1956   
1957     // Add primary vertex. The primary tracks will be defined
1958     // after the loops on the composite objects (V0, cascades, kinks)
1959     const AliESDVertex *vtx = esd->GetPrimaryVertex();
1960       
1961     vtx->GetXYZ(pos); // position
1962     vtx->GetCovMatrix(covVtx); //covariance matrix
1963
1964     AliAODVertex * primary = new(vertices[jVertices++])
1965       AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1966          
1967     // Create vertices starting from the most complex objects
1968       
1969     // Cascades
1970     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
1971       AliESDcascade *cascade = esd->GetCascade(nCascade);
1972       
1973       cascade->GetXYZcascade(pos[0], pos[1], pos[2]);  // Bo: bug correction
1974       cascade->GetPosCovXi(covVtx);
1975      
1976       // Add the cascade vertex
1977       AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
1978                                                                         covVtx,
1979                                                                         cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
1980                                                                         primary,
1981                                                                         nCascade,
1982                                                                         AliAODVertex::kCascade);
1983
1984       primary->AddDaughter(vcascade);
1985
1986       // Add the V0 from the cascade. The ESD class have to be optimized...
1987       // Now we have to search for the corresponding V0 in the list of V0s
1988       // using the indeces of the positive and negative tracks
1989
1990       Int_t posFromV0 = cascade->GetPindex();
1991       Int_t negFromV0 = cascade->GetNindex();
1992
1993       AliESDv0 * v0 = 0x0;
1994       Int_t indV0 = -1;
1995
1996       for (Int_t iV0=0; iV0<nV0s; ++iV0) {
1997
1998         v0 = esd->GetV0(iV0);
1999         Int_t posV0 = v0->GetPindex();
2000         Int_t negV0 = v0->GetNindex();
2001
2002         if (posV0==posFromV0 && negV0==negFromV0) {
2003           indV0 = iV0;
2004           break;
2005         }
2006       }
2007
2008       AliAODVertex * vV0FromCascade = 0x0;
2009
2010       if (indV0>-1 && !usedV0[indV0] ) {
2011         
2012         // the V0 exists in the array of V0s and is not used
2013
2014         usedV0[indV0] = kTRUE;
2015         
2016         v0->GetXYZ(pos[0], pos[1], pos[2]);
2017         v0->GetPosCov(covVtx);
2018         
2019         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2020                                                                  covVtx,
2021                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2022                                                                  vcascade,
2023                                                                  indV0,
2024                                                                  AliAODVertex::kV0);
2025       } else {
2026
2027         // the V0 doesn't exist in the array of V0s or was used
2028         cerr << "Error: event " << iEvent << " cascade " << nCascade
2029              << " The V0 " << indV0 
2030              << " doesn't exist in the array of V0s or was used!" << endl;
2031
2032         cascade->GetXYZ(pos[0], pos[1], pos[2]);
2033         cascade->GetPosCov(covVtx);
2034       
2035         vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2036                                                                  covVtx,
2037                                                                  v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2038                                                                  vcascade,
2039                                                                  indV0,
2040                                                                  AliAODVertex::kV0);
2041         vcascade->AddDaughter(vV0FromCascade);
2042       }
2043
2044       // Add the positive tracks from the V0
2045
2046       if (! usedTrack[posFromV0]) {
2047
2048         usedTrack[posFromV0] = kTRUE;
2049
2050         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2051         esdTrack->GetPxPyPz(p);
2052         esdTrack->GetXYZ(pos);
2053         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2054         esdTrack->GetESDpid(pid);
2055         
2056         vV0FromCascade->AddDaughter(aodTrack =
2057                                     new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2058                                            esdTrack->GetLabel(), 
2059                                            p, 
2060                                            kTRUE,
2061                                            pos,
2062                                            kFALSE,
2063                                            covTr, 
2064                                            (Short_t)esdTrack->Charge(),
2065                                            esdTrack->GetITSClusterMap(), 
2066                                            pid,
2067                                            vV0FromCascade,
2068                                            kTRUE,  // check if this is right
2069                                            kFALSE, // check if this is right
2070                                            AliAODTrack::kSecondary)
2071                 );
2072         aodTrack->ConvertAliPIDtoAODPID();
2073       }
2074       else {
2075         cerr << "Error: event " << iEvent << " cascade " << nCascade
2076              << " track " << posFromV0 << " has already been used!" << endl;
2077       }
2078
2079       // Add the negative tracks from the V0
2080
2081       if (!usedTrack[negFromV0]) {
2082         
2083         usedTrack[negFromV0] = kTRUE;
2084         
2085         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2086         esdTrack->GetPxPyPz(p);
2087         esdTrack->GetXYZ(pos);
2088         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2089         esdTrack->GetESDpid(pid);
2090         
2091         vV0FromCascade->AddDaughter(aodTrack =
2092                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2093                                            esdTrack->GetLabel(),
2094                                            p,
2095                                            kTRUE,
2096                                            pos,
2097                                            kFALSE,
2098                                            covTr, 
2099                                            (Short_t)esdTrack->Charge(),
2100                                            esdTrack->GetITSClusterMap(), 
2101                                            pid,
2102                                            vV0FromCascade,
2103                                            kTRUE,  // check if this is right
2104                                            kFALSE, // check if this is right
2105                                            AliAODTrack::kSecondary)
2106                 );
2107         aodTrack->ConvertAliPIDtoAODPID();
2108       }
2109       else {
2110         cerr << "Error: event " << iEvent << " cascade " << nCascade
2111              << " track " << negFromV0 << " has already been used!" << endl;
2112       }
2113
2114       // Add the bachelor track from the cascade
2115
2116       Int_t bachelor = cascade->GetBindex();
2117       
2118       if(!usedTrack[bachelor]) {
2119       
2120         usedTrack[bachelor] = kTRUE;
2121         
2122         AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2123         esdTrack->GetPxPyPz(p);
2124         esdTrack->GetXYZ(pos);
2125         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2126         esdTrack->GetESDpid(pid);
2127
2128         vcascade->AddDaughter(aodTrack =
2129                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2130                                            esdTrack->GetLabel(),
2131                                            p,
2132                                            kTRUE,
2133                                            pos,
2134                                            kFALSE,
2135                                            covTr, 
2136                                            (Short_t)esdTrack->Charge(),
2137                                            esdTrack->GetITSClusterMap(), 
2138                                            pid,
2139                                            vcascade,
2140                                            kTRUE,  // check if this is right
2141                                            kFALSE, // check if this is right
2142                                            AliAODTrack::kSecondary)
2143                 );
2144         aodTrack->ConvertAliPIDtoAODPID();
2145      }
2146       else {
2147         cerr << "Error: event " << iEvent << " cascade " << nCascade
2148              << " track " << bachelor << " has already been used!" << endl;
2149       }
2150
2151       // Add the primary track of the cascade (if any)
2152
2153     } // end of the loop on cascades
2154     
2155     // V0s
2156         
2157     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2158
2159       if (usedV0[nV0]) continue; // skip if aready added to the AOD
2160
2161       AliESDv0 *v0 = esd->GetV0(nV0);
2162       
2163       v0->GetXYZ(pos[0], pos[1], pos[2]);
2164       v0->GetPosCov(covVtx);
2165
2166       AliAODVertex * vV0 = 
2167         new(vertices[jVertices++]) AliAODVertex(pos,
2168                                                 covVtx,
2169                                                 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2170                                                 primary,
2171                                                 nV0,
2172                                                 AliAODVertex::kV0);
2173       primary->AddDaughter(vV0);
2174
2175       Int_t posFromV0 = v0->GetPindex();
2176       Int_t negFromV0 = v0->GetNindex();
2177       
2178       // Add the positive tracks from the V0
2179
2180       if (!usedTrack[posFromV0]) {
2181         
2182         usedTrack[posFromV0] = kTRUE;
2183
2184         AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2185         esdTrack->GetPxPyPz(p);
2186         esdTrack->GetXYZ(pos);
2187         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2188         esdTrack->GetESDpid(pid);
2189         
2190         vV0->AddDaughter(aodTrack =
2191                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2192                                            esdTrack->GetLabel(), 
2193                                            p, 
2194                                            kTRUE,
2195                                            pos,
2196                                            kFALSE,
2197                                            covTr, 
2198                                            (Short_t)esdTrack->Charge(),
2199                                            esdTrack->GetITSClusterMap(), 
2200                                            pid,
2201                                            vV0,
2202                                            kTRUE,  // check if this is right
2203                                            kFALSE, // check if this is right
2204                                            AliAODTrack::kSecondary)
2205                 );
2206         aodTrack->ConvertAliPIDtoAODPID();
2207       }
2208       else {
2209         cerr << "Error: event " << iEvent << " V0 " << nV0
2210              << " track " << posFromV0 << " has already been used!" << endl;
2211       }
2212
2213       // Add the negative tracks from the V0
2214
2215       if (!usedTrack[negFromV0]) {
2216
2217         usedTrack[negFromV0] = kTRUE;
2218
2219         AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2220         esdTrack->GetPxPyPz(p);
2221         esdTrack->GetXYZ(pos);
2222         esdTrack->GetCovarianceXYZPxPyPz(covTr);
2223         esdTrack->GetESDpid(pid);
2224
2225         vV0->AddDaughter(aodTrack =
2226                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2227                                            esdTrack->GetLabel(),
2228                                            p,
2229                                            kTRUE,
2230                                            pos,
2231                                            kFALSE,
2232                                            covTr, 
2233                                            (Short_t)esdTrack->Charge(),
2234                                            esdTrack->GetITSClusterMap(), 
2235                                            pid,
2236                                            vV0,
2237                                            kTRUE,  // check if this is right
2238                                            kFALSE, // check if this is right
2239                                            AliAODTrack::kSecondary)
2240                 );
2241         aodTrack->ConvertAliPIDtoAODPID();
2242       }
2243       else {
2244         cerr << "Error: event " << iEvent << " V0 " << nV0
2245              << " track " << negFromV0 << " has already been used!" << endl;
2246       }
2247
2248     } // end of the loop on V0s
2249     
2250     // Kinks: it is a big mess the access to the information in the kinks
2251     // The loop is on the tracks in order to find the mother and daugther of each kink
2252
2253
2254     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2255
2256
2257       AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2258
2259       Int_t ikink = esdTrack->GetKinkIndex(0);
2260
2261       if (ikink) {
2262         // Negative kink index: mother, positive: daughter
2263
2264         // Search for the second track of the kink
2265
2266         for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2267
2268           AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2269
2270           Int_t jkink = esdTrack1->GetKinkIndex(0);
2271
2272           if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2273
2274             // The two tracks are from the same kink
2275           
2276             if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2277
2278             Int_t imother = -1;
2279             Int_t idaughter = -1;
2280
2281             if (ikink<0 && jkink>0) {
2282
2283               imother = iTrack;
2284               idaughter = jTrack;
2285             }
2286             else if (ikink>0 && jkink<0) {
2287
2288               imother = jTrack;
2289               idaughter = iTrack;
2290             }
2291             else {
2292               cerr << "Error: Wrong combination of kink indexes: "
2293               << ikink << " " << jkink << endl;
2294               continue;
2295             }
2296
2297             // Add the mother track
2298
2299             AliAODTrack * mother = NULL;
2300
2301             if (!usedTrack[imother]) {
2302         
2303               usedTrack[imother] = kTRUE;
2304         
2305               AliESDtrack *esdTrack = esd->GetTrack(imother);
2306               esdTrack->GetPxPyPz(p);
2307               esdTrack->GetXYZ(pos);
2308               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2309               esdTrack->GetESDpid(pid);
2310
2311               mother = 
2312                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2313                                            esdTrack->GetLabel(),
2314                                            p,
2315                                            kTRUE,
2316                                            pos,
2317                                            kFALSE,
2318                                            covTr, 
2319                                            (Short_t)esdTrack->Charge(),
2320                                            esdTrack->GetITSClusterMap(), 
2321                                            pid,
2322                                            primary,
2323                                            kTRUE, // check if this is right
2324                                            kTRUE, // check if this is right
2325                                            AliAODTrack::kPrimary);
2326               primary->AddDaughter(mother);
2327               mother->ConvertAliPIDtoAODPID();
2328             }
2329             else {
2330               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2331               << " track " << imother << " has already been used!" << endl;
2332             }
2333
2334             // Add the kink vertex
2335             AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2336
2337             AliAODVertex * vkink = 
2338             new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2339                                                     NULL,
2340                                                     0.,
2341                                                     mother,
2342                                                     esdTrack->GetID(), // This is the track ID of the mother's track!
2343                                                     AliAODVertex::kKink);
2344             // Add the daughter track
2345
2346             AliAODTrack * daughter = NULL;
2347
2348             if (!usedTrack[idaughter]) {
2349         
2350               usedTrack[idaughter] = kTRUE;
2351         
2352               AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2353               esdTrack->GetPxPyPz(p);
2354               esdTrack->GetXYZ(pos);
2355               esdTrack->GetCovarianceXYZPxPyPz(covTr);
2356               esdTrack->GetESDpid(pid);
2357
2358               daughter = 
2359                 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2360                                            esdTrack->GetLabel(),
2361                                            p,
2362                                            kTRUE,
2363                                            pos,
2364                                            kFALSE,
2365                                            covTr, 
2366                                            (Short_t)esdTrack->Charge(),
2367                                            esdTrack->GetITSClusterMap(), 
2368                                            pid,
2369                                            vkink,
2370                                            kTRUE, // check if this is right
2371                                            kTRUE, // check if this is right
2372                                            AliAODTrack::kPrimary);
2373               vkink->AddDaughter(daughter);
2374               daughter->ConvertAliPIDtoAODPID();
2375             }
2376             else {
2377               cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2378               << " track " << idaughter << " has already been used!" << endl;
2379             }
2380
2381
2382           }
2383         }
2384
2385       }      
2386
2387     }
2388
2389     
2390     // Tracks (primary and orphan)
2391       
2392     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2393         
2394
2395       if (usedTrack[nTrack]) continue;
2396
2397       AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2398       esdTrack->GetPxPyPz(p);
2399       esdTrack->GetXYZ(pos);
2400       esdTrack->GetCovarianceXYZPxPyPz(covTr);
2401       esdTrack->GetESDpid(pid);
2402
2403       Float_t impactXY, impactZ;
2404
2405       esdTrack->GetImpactParameters(impactXY,impactZ);
2406
2407       if (impactXY<3) {
2408         // track inside the beam pipe
2409       
2410         primary->AddDaughter(aodTrack =
2411             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2412                                          esdTrack->GetLabel(),
2413                                          p,
2414                                          kTRUE,
2415                                          pos,
2416                                          kFALSE,
2417                                          covTr, 
2418                                          (Short_t)esdTrack->Charge(),
2419                                          esdTrack->GetITSClusterMap(), 
2420                                          pid,
2421                                          primary,
2422                                          kTRUE, // check if this is right
2423                                          kTRUE, // check if this is right
2424                                          AliAODTrack::kPrimary)
2425             );
2426         aodTrack->ConvertAliPIDtoAODPID();
2427       }
2428       else {
2429         // outside the beam pipe: orphan track
2430             aodTrack =
2431             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2432                                          esdTrack->GetLabel(),
2433                                          p,
2434                                          kTRUE,
2435                                          pos,
2436                                          kFALSE,
2437                                          covTr, 
2438                                          (Short_t)esdTrack->Charge(),
2439                                          esdTrack->GetITSClusterMap(), 
2440                                          pid,
2441                                          NULL,
2442                                          kFALSE, // check if this is right
2443                                          kFALSE, // check if this is right
2444                                          AliAODTrack::kOrphan);
2445             aodTrack->ConvertAliPIDtoAODPID();
2446       } 
2447     } // end of loop on tracks
2448
2449     // muon tracks
2450     Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2451     for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2452       
2453       AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);     
2454       p[0] = esdMuTrack->Px(); 
2455       p[1] = esdMuTrack->Py(); 
2456       p[2] = esdMuTrack->Pz();
2457       pos[0] = primary->GetX(); 
2458       pos[1] = primary->GetY(); 
2459       pos[2] = primary->GetZ();
2460       
2461       // has to be changed once the muon pid is provided by the ESD
2462       for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2463       
2464       primary->AddDaughter(aodTrack =
2465           new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2466                                              0, // no label provided
2467                                              p,
2468                                              kTRUE,
2469                                              pos,
2470                                              kFALSE,
2471                                              NULL, // no covariance matrix provided
2472                                              esdMuTrack->Charge(),
2473                                              0, // ITSClusterMap is set below
2474                                              pid,
2475                                              primary,
2476                                              kFALSE,  // muon tracks are not used to fit the primary vtx
2477                                              kFALSE,  // not used for vertex fit
2478                                              AliAODTrack::kPrimary)
2479           );
2480
2481       aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
2482       Int_t track2Trigger = esdMuTrack->GetMatchTrigger();
2483       aodTrack->SetMatchTrigger(track2Trigger);
2484       if (track2Trigger) 
2485         aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
2486       else 
2487         aodTrack->SetChi2MatchTrigger(0.);
2488     }
2489     
2490     // Access to the AOD container of clusters
2491     TClonesArray &clusters = *(aod->GetClusters());
2492     Int_t jClusters=0;
2493
2494     // Calo Clusters
2495     Int_t nClusters    = esd->GetNumberOfCaloClusters();
2496
2497     for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2498
2499       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2500
2501       Int_t id = cluster->GetID();
2502       Int_t label = -1;
2503       Float_t energy = cluster->E();
2504       cluster->GetPosition(posF);
2505       AliAODVertex *prodVertex = primary;
2506       AliAODTrack *primTrack = NULL;
2507       Char_t ttype=AliAODCluster::kUndef;
2508
2509       if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2510       else if (cluster->IsEMCAL()) {
2511
2512         if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2513           ttype = AliAODCluster::kEMCALPseudoCluster;
2514         else
2515           ttype = AliAODCluster::kEMCALClusterv1;
2516
2517       }
2518
2519       new(clusters[jClusters++]) AliAODCluster(id,
2520                                                label,
2521                                                energy,
2522                                                pos,
2523                                                NULL, // no covariance matrix provided
2524                                                NULL, // no pid for clusters provided
2525                                                prodVertex,
2526                                                primTrack,
2527                                                ttype);
2528
2529     } // end of loop on calo clusters
2530
2531     // tracklets
2532     const AliMultiplicity *mult = esd->GetMultiplicity();
2533     if (mult) {
2534       if (mult->GetNumberOfTracklets()>0) {
2535         aod->GetTracklets()->CreateContainer(mult->GetNumberOfTracklets());
2536
2537         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
2538           aod->GetTracklets()->SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n));
2539         }
2540       }
2541     } else {
2542       Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
2543     }
2544
2545     delete [] usedTrack;
2546     delete [] usedV0;
2547     delete [] usedKink;
2548
2549     // fill the tree for this event
2550     aodTree->Fill();
2551   } // end of event loop
2552
2553   aodTree->GetUserInfo()->Add(aod);
2554
2555   // close ESD file
2556   esdFile->Close();
2557
2558   // write the tree to the specified file
2559   aodFile = aodTree->GetCurrentFile();
2560   aodFile->cd();
2561   aodTree->Write();
2562
2563   return;
2564 }
2565
2566 void AliReconstruction::WriteAlignmentData(AliESDEvent* esd)
2567 {
2568   // Write space-points which are then used in the alignment procedures
2569   // For the moment only ITS, TRD and TPC
2570
2571   // Load TOF clusters
2572   if (fTracker[3]){
2573     fLoader[3]->LoadRecPoints("read");
2574     TTree* tree = fLoader[3]->TreeR();
2575     if (!tree) {
2576       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2577       return;
2578     }
2579     fTracker[3]->LoadClusters(tree);
2580   }
2581   Int_t ntracks = esd->GetNumberOfTracks();
2582   for (Int_t itrack = 0; itrack < ntracks; itrack++)
2583     {
2584       AliESDtrack *track = esd->GetTrack(itrack);
2585       Int_t nsp = 0;
2586       Int_t idx[200];
2587       for (Int_t iDet = 3; iDet >= 0; iDet--)
2588         nsp += track->GetNcls(iDet);
2589       if (nsp) {
2590         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2591         track->SetTrackPointArray(sp);
2592         Int_t isptrack = 0;
2593         for (Int_t iDet = 3; iDet >= 0; iDet--) {
2594           AliTracker *tracker = fTracker[iDet];
2595           if (!tracker) continue;
2596           Int_t nspdet = track->GetNcls(iDet);
2597           if (nspdet <= 0) continue;
2598           track->GetClusters(iDet,idx);
2599           AliTrackPoint p;
2600           Int_t isp = 0;
2601           Int_t isp2 = 0;
2602           while (isp < nspdet) {
2603             Bool_t isvalid;
2604             if(IsSelected(fgkDetectorName[iDet],fUseTrackingErrorsForAlignment)) {
2605               isvalid = tracker->GetTrackPointTrackingError(idx[isp2],p,track);
2606             } else {
2607               isvalid = tracker->GetTrackPoint(idx[isp2],p); 
2608             } 
2609             isp2++;
2610             const Int_t kNTPCmax = 159;
2611             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
2612             if (!isvalid) continue;
2613             sp->AddPoint(isptrack,&p); isptrack++; isp++;
2614           }
2615         }       
2616       }
2617     }
2618   if (fTracker[3]){
2619     fTracker[3]->UnloadClusters();
2620     fLoader[3]->UnloadRecPoints();
2621   }
2622 }
2623
2624 //_____________________________________________________________________________
2625 void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESDEvent* esd)
2626 {
2627   // The method reads the raw-data error log
2628   // accumulated within the rawReader.
2629   // It extracts the raw-data errors related to
2630   // the current event and stores them into
2631   // a TClonesArray inside the esd object.
2632
2633   if (!fRawReader) return;
2634
2635   for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2636
2637     AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2638     if (!log) continue;
2639     if (iEvent != log->GetEventNumber()) continue;
2640
2641     esd->AddRawDataErrorLog(log);
2642   }
2643
2644 }
2645
2646 TNamed* AliReconstruction::CopyFileToTNamed(TString fPath,TString fName){
2647   // Dump a file content into a char in TNamed
2648   ifstream in;
2649   in.open(fPath.Data(),ios::in | ios::binary|ios::ate);
2650   Int_t kBytes = (Int_t)in.tellg();
2651   printf("Size: %d \n",kBytes);
2652   TNamed *fn = 0;
2653   if(in.good()){
2654     char* memblock = new char [kBytes];
2655     in.seekg (0, ios::beg);
2656     in.read (memblock, kBytes);
2657     in.close();
2658     TString fData(memblock,kBytes);
2659     fn = new TNamed(fName,fData);
2660     printf("fData Size: %d \n",fData.Sizeof());
2661     printf("fName Size: %d \n",fName.Sizeof());
2662     printf("fn    Size: %d \n",fn->Sizeof());
2663     delete[] memblock;
2664   }
2665   else{
2666     AliInfo(Form("Could not Open %s\n",fPath.Data()));
2667   }
2668
2669   return fn;
2670 }
2671
2672 void AliReconstruction::TNamedToFile(TTree* fTree, TString fName){
2673   // This is not really needed in AliReconstruction at the moment
2674   // but can serve as a template
2675
2676   TList *fList = fTree->GetUserInfo();
2677   TNamed *fn = (TNamed*)fList->FindObject(fName.Data());
2678   printf("fn Size: %d \n",fn->Sizeof());
2679
2680   TString fTmp(fn->GetName()); // to be 100% sure in principle fName also works
2681   const char* cdata = fn->GetTitle();
2682   printf("fTmp Size %d\n",fTmp.Sizeof());
2683
2684   int size = fn->Sizeof()-fTmp.Sizeof()-sizeof(UChar_t)-sizeof(Int_t); // see dfinition of TString::SizeOf()...
2685   printf("calculated size %d\n",size);
2686   ofstream out(fName.Data(),ios::out | ios::binary);
2687   out.write(cdata,size);
2688   out.close();
2689
2690 }
2691
2692 //_____________________________________________________________________________
2693 AliQualAssDataMaker * AliReconstruction::GetQualAssDataMaker(Int_t iDet)
2694 {
2695 // get the quality assurance data maker object and the loader for a detector
2696
2697   if (fQualAssDataMaker[iDet]) 
2698     return fQualAssDataMaker[iDet];
2699
2700   // load the QA data maker object
2701   TPluginManager* pluginManager = gROOT->GetPluginManager();
2702   TString detName = fgkDetectorName[iDet];
2703   TString qadmName = "Ali" + detName + "QualAssDataMaker";
2704   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) 
2705     return NULL;
2706
2707   AliQualAssDataMaker * qadm = NULL;
2708   // first check if a plugin is defined for the quality assurance data maker
2709   TPluginHandler* pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2710   // if not, add a plugin for it
2711   if (!pluginHandler) {
2712     AliDebug(1, Form("defining plugin for %s", qadmName.Data()));
2713     TString libs = gSystem->GetLibraries();
2714     if (libs.Contains("lib" + detName + "base.so") ||
2715         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
2716       pluginManager->AddHandler("AliQualAssDataMaker", detName, 
2717                                 qadmName, detName + "qadm", qadmName + "()");
2718     } else {
2719       pluginManager->AddHandler("AliQualAssDataMaker", detName, 
2720                                 qadmName, detName, qadmName + "()");
2721     }
2722     pluginHandler = pluginManager->FindHandler("AliQualAssDataMaker", detName);
2723   }
2724   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
2725     qadm = (AliQualAssDataMaker *) pluginHandler->ExecPlugin(0);
2726   }
2727   if (qadm) {
2728     AliInfo(Form("Initializing quality assurance data maker for %s", fgkDetectorName[iDet]));
2729     qadm->Init(AliQualAss::kRECPOINTS, fRunLoader->GetHeader()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2730     qadm->Init(AliQualAss::kESDS, fRunLoader->GetHeader()->GetRun(), GetQACycles(fgkDetectorName[iDet]));
2731     qadm->StartOfCycle(AliQualAss::kRECPOINTS);
2732     qadm->StartOfCycle(AliQualAss::kESDS, "same") ;     
2733     fQualAssDataMaker[iDet] = qadm;
2734   }
2735
2736   return qadm;
2737 }
2738
2739 //_____________________________________________________________________________
2740 Bool_t AliReconstruction::RunQualAss(const char* detectors, AliESDEvent *& esd)
2741 {
2742   // run the Quality Assurance data producer
2743
2744   AliCodeTimerAuto("")
2745   TString detStr = detectors;
2746   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
2747    if (!IsSelected(fgkDetectorName[iDet], detStr)) 
2748      continue;
2749    AliQualAssDataMaker * qadm = GetQualAssDataMaker(iDet);
2750    if (!qadm) 
2751      continue;
2752    AliCodeTimerStart(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2753    AliInfo(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2754     
2755    qadm->Exec(AliQualAss::kESDS, esd) ; 
2756    qadm->Increment() ; 
2757
2758    AliCodeTimerStop(Form("running quality assurance data maker for %s", fgkDetectorName[iDet]));
2759  }
2760  if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
2761    AliError(Form("the following detectors were not found: %s",
2762                  detStr.Data()));
2763    if (fStopOnError) 
2764      return kFALSE;
2765  }
2766  
2767  return kTRUE;
2768   
2769
2770 }
2771
2772 //_____________________________________________________________________________
2773 Int_t AliReconstruction::GetDetIndex(const char* detector)
2774 {
2775   // return the detector index corresponding to detector
2776   Int_t index = -1 ; 
2777   for (index = 0; index < fgkNDetectors ; index++) {
2778     if ( strcmp(detector, fgkDetectorName[index]) == 0 )
2779         break ; 
2780   }     
2781   return index ; 
2782 }