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