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