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