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