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