Calling the V0 and cascade finders from AliReconstruction, and the possibility 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 // The name of the galice file can be changed from the default               //
51 // "galice.root" by passing it as argument to the AliReconstruction          //
52 // constructor or by                                                         //
53 //                                                                           //
54 //   rec.SetGAliceFile("...");                                               //
55 //                                                                           //
56 // The local reconstruction can be switched on or off for individual         //
57 // detectors by                                                              //
58 //                                                                           //
59 //   rec.SetRunLocalReconstruction("...");                                   //
60 //                                                                           //
61 // The argument is a (case sensitive) string with the names of the           //
62 // detectors separated by a space. The special string "ALL" selects all      //
63 // available detectors. This is the default.                                 //
64 //                                                                           //
65 // The reconstruction of the primary vertex position can be switched off by  //
66 //                                                                           //
67 //   rec.SetRunVertexFinder(kFALSE);                                         //
68 //                                                                           //
69 // The tracking and the creation of ESD tracks can be switched on for        //
70 // selected detectors by                                                     //
71 //                                                                           //
72 //   rec.SetRunTracking("...");                                              //
73 //                                                                           //
74 // Uniform/nonuniform field tracking switches (default: uniform field)       //
75 //                                                                           //
76 //   rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
77 //                                                                           //
78 // The filling of additional ESD information can be steered by               //
79 //                                                                           //
80 //   rec.SetFillESD("...");                                                  //
81 //                                                                           //
82 // Again, for both methods the string specifies the list of detectors.       //
83 // The default is "ALL".                                                     //
84 //                                                                           //
85 // The call of the shortcut method                                           //
86 //                                                                           //
87 //   rec.SetRunReconstruction("...");                                        //
88 //                                                                           //
89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
90 // SetFillESD with the same detector selecting string as argument.           //
91 //                                                                           //
92 // The reconstruction requires digits or raw data as input. For the creation //
93 // of digits and raw data have a look at the class AliSimulation.            //
94 //                                                                           //
95 // For debug purposes the method SetCheckPointLevel can be used. If the      //
96 // argument is greater than 0, files with ESD events will be written after   //
97 // selected steps of the reconstruction for each event:                      //
98 //   level 1: after tracking and after filling of ESD (final)                //
99 //   level 2: in addition after each tracking step                           //
100 //   level 3: in addition after the filling of ESD for each detector         //
101 // If a final check point file exists for an event, this event will be       //
102 // skipped in the reconstruction. The tracking and the filling of ESD for    //
103 // a detector will be skipped as well, if the corresponding check point      //
104 // file exists. The ESD event will then be loaded from the file instead.     //
105 //                                                                           //
106 ///////////////////////////////////////////////////////////////////////////////
107
108 #include <TArrayF.h>
109 #include <TFile.h>
110 #include <TSystem.h>
111 #include <TROOT.h>
112 #include <TPluginManager.h>
113 #include <TStopwatch.h>
114 #include <TGeoManager.h>
115 #include <TLorentzVector.h>
116
117 #include "AliReconstruction.h"
118 #include "AliReconstructor.h"
119 #include "AliLog.h"
120 #include "AliRunLoader.h"
121 #include "AliRun.h"
122 #include "AliRawReaderFile.h"
123 #include "AliRawReaderDate.h"
124 #include "AliRawReaderRoot.h"
125 #include "AliRawEventHeaderBase.h"
126 #include "AliESD.h"
127 #include "AliESDfriend.h"
128 #include "AliESDVertex.h"
129 #include "AliMultiplicity.h"
130 #include "AliTracker.h"
131 #include "AliVertexer.h"
132 #include "AliVertexerTracks.h"
133 #include "AliV0vertexer.h"
134 #include "AliCascadeVertexer.h"
135 #include "AliHeader.h"
136 #include "AliGenEventHeader.h"
137 #include "AliPID.h"
138 #include "AliESDpid.h"
139 #include "AliESDtrack.h"
140
141 #include "AliRunTag.h"
142 #include "AliDetectorTag.h"
143 #include "AliEventTag.h"
144
145 #include "AliTrackPointArray.h"
146 #include "AliCDBManager.h"
147 #include "AliCDBEntry.h"
148 #include "AliAlignObj.h"
149
150 #include "AliCentralTrigger.h"
151 #include "AliCTPRawStream.h"
152
153 ClassImp(AliReconstruction)
154
155
156 //_____________________________________________________________________________
157 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
158
159 //_____________________________________________________________________________
160 AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
161                                      const char* name, const char* title) :
162   TNamed(name, title),
163
164   fUniformField(kTRUE),
165   fRunVertexFinder(kTRUE),
166   fRunHLTTracking(kFALSE),
167   fStopOnError(kFALSE),
168   fWriteAlignmentData(kFALSE),
169   fWriteESDfriend(kFALSE),
170   fFillTriggerESD(kTRUE),
171
172   fRunLocalReconstruction("ALL"),
173   fRunTracking("ALL"),
174   fFillESD("ALL"),
175   fGAliceFileName(gAliceFilename),
176   fInput(""),
177   fEquipIdMap(""),
178   fFirstEvent(0),
179   fLastEvent(-1),
180   fCheckPointLevel(0),
181   fOptions(),
182   fLoadAlignFromCDB(kTRUE),
183   fLoadAlignData("ALL"),
184
185   fRunLoader(NULL),
186   fRawReader(NULL),
187
188   fVertexer(NULL),
189   fDiamondProfile(NULL),
190
191   fAlignObjArray(NULL),
192   fCDBUri(cdbUri),
193   fSpecCDBUri()
194 {
195 // create reconstruction object with default parameters
196   
197   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
198     fReconstructor[iDet] = NULL;
199     fLoader[iDet] = NULL;
200     fTracker[iDet] = NULL;
201   }
202   AliPID pid;
203 }
204
205 //_____________________________________________________________________________
206 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
207   TNamed(rec),
208
209   fUniformField(rec.fUniformField),
210   fRunVertexFinder(rec.fRunVertexFinder),
211   fRunHLTTracking(rec.fRunHLTTracking),
212   fStopOnError(rec.fStopOnError),
213   fWriteAlignmentData(rec.fWriteAlignmentData),
214   fWriteESDfriend(rec.fWriteESDfriend),
215   fFillTriggerESD(rec.fFillTriggerESD),
216
217   fRunLocalReconstruction(rec.fRunLocalReconstruction),
218   fRunTracking(rec.fRunTracking),
219   fFillESD(rec.fFillESD),
220   fGAliceFileName(rec.fGAliceFileName),
221   fInput(rec.fInput),
222   fEquipIdMap(rec.fEquipIdMap),
223   fFirstEvent(rec.fFirstEvent),
224   fLastEvent(rec.fLastEvent),
225   fCheckPointLevel(0),
226   fOptions(),
227   fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
228   fLoadAlignData(rec.fLoadAlignData),
229
230   fRunLoader(NULL),
231   fRawReader(NULL),
232
233   fVertexer(NULL),
234   fDiamondProfile(NULL),
235
236   fAlignObjArray(rec.fAlignObjArray),
237   fCDBUri(rec.fCDBUri),
238   fSpecCDBUri()
239 {
240 // copy constructor
241
242   for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
243     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
244   }
245   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
246     fReconstructor[iDet] = NULL;
247     fLoader[iDet] = NULL;
248     fTracker[iDet] = NULL;
249   }
250   for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
251     if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
252   }
253 }
254
255 //_____________________________________________________________________________
256 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
257 {
258 // assignment operator
259
260   this->~AliReconstruction();
261   new(this) AliReconstruction(rec);
262   return *this;
263 }
264
265 //_____________________________________________________________________________
266 AliReconstruction::~AliReconstruction()
267 {
268 // clean up
269
270   CleanUp();
271   fOptions.Delete();
272   fSpecCDBUri.Delete();
273 }
274
275 //_____________________________________________________________________________
276 void AliReconstruction::InitCDBStorage()
277 {
278 // activate a default CDB storage
279 // First check if we have any CDB storage set, because it is used 
280 // to retrieve the calibration and alignment constants
281
282   AliCDBManager* man = AliCDBManager::Instance();
283   if (man->IsDefaultStorageSet())
284   {
285     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
286     AliWarning("Default CDB storage has been already set !");
287     AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
288     AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
289     fCDBUri = "";
290   }
291   else {
292     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
293     AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
294     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
295     man->SetDefaultStorage(fCDBUri);
296   }
297
298   // Now activate the detector specific CDB storage locations
299   for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
300     TObject* obj = fSpecCDBUri[i];
301     if (!obj) continue;
302     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
303     AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
304     AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
305     man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
306   }
307   man->Print();
308 }
309
310 //_____________________________________________________________________________
311 void AliReconstruction::SetDefaultStorage(const char* uri) {
312 // Store the desired default CDB storage location
313 // Activate it later within the Run() method
314
315   fCDBUri = uri;
316
317 }
318
319 //_____________________________________________________________________________
320 void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
321 // Store a detector-specific CDB storage location
322 // Activate it later within the Run() method
323
324   AliCDBPath aPath(calibType);
325   if(!aPath.IsValid()){
326         // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
327         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
328                 if(!strcmp(calibType, fgkDetectorName[iDet])) {
329                         aPath.SetPath(Form("%s/*", calibType));
330                         AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
331                         break;
332                 }
333         }
334         if(!aPath.IsValid()){
335                 AliError(Form("Not a valid path or detector: %s", calibType));
336                 return;
337         }
338   }
339
340   // check that calibType refers to a "valid" detector name
341   Bool_t isDetector = kFALSE;
342   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
343     TString detName = fgkDetectorName[iDet];
344     if(aPath.GetLevel0() == detName) {
345         isDetector = kTRUE;
346         break;
347     }
348   }
349
350   if(!isDetector) {
351         AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
352         return;
353   }
354
355   TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
356   if (obj) fSpecCDBUri.Remove(obj);
357   fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
358
359 }
360
361 //_____________________________________________________________________________
362 Bool_t AliReconstruction::SetRunNumber()
363 {
364   // The method is called in Run() in order
365   // to set a correct run number.
366   // In case of raw data reconstruction the
367   // run number is taken from the raw data header
368
369   if(AliCDBManager::Instance()->GetRun() < 0) {
370     if (!fRunLoader) {
371       AliError("No run loader is found !"); 
372       return kFALSE;
373     }
374     // read run number from gAlice
375     if(fRunLoader->GetAliRun())
376       AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
377     else {
378       if(fRawReader) {
379         if(fRawReader->NextEvent()) {
380           AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
381           fRawReader->RewindEvents();
382         }
383         else {
384           AliError("No raw-data events found !");
385           return kFALSE;
386         }
387       }
388       else {
389         AliError("Neither gAlice nor RawReader objects are found !");
390         return kFALSE;
391       }
392     }
393     AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
394   }
395   return kTRUE;
396 }
397
398 //_____________________________________________________________________________
399 Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
400 {
401   // Read collection of alignment objects (AliAlignObj derived) saved
402   // in the TClonesArray ClArrayName and apply them to the geometry
403   // manager singleton.
404   //
405   alObjArray->Sort();
406   Int_t nvols = alObjArray->GetEntriesFast();
407
408   Bool_t flag = kTRUE;
409
410   for(Int_t j=0; j<nvols; j++)
411     {
412       AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
413       if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
414     }
415
416   if (AliDebugLevelClass() >= 1) {
417     gGeoManager->GetTopNode()->CheckOverlaps(20);
418     TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
419     if(ovexlist->GetEntriesFast()){  
420       AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
421    }
422   }
423
424   return flag;
425
426 }
427
428 //_____________________________________________________________________________
429 Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
430 {
431   // Fills array of single detector's alignable objects from CDB
432   
433   AliDebug(2, Form("Loading alignment data for detector: %s",detName));
434   
435   AliCDBEntry *entry;
436         
437   AliCDBPath path(detName,"Align","Data");
438         
439   entry=AliCDBManager::Instance()->Get(path.GetPath());
440   if(!entry){ 
441         AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
442         return kFALSE;
443   }
444   entry->SetOwner(1);
445   TClonesArray *alignArray = (TClonesArray*) entry->GetObject();        
446   alignArray->SetOwner(0);
447   AliDebug(2,Form("Found %d alignment objects for %s",
448                         alignArray->GetEntries(),detName));
449
450   AliAlignObj *alignObj=0;
451   TIter iter(alignArray);
452         
453   // loop over align objects in detector
454   while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
455         fAlignObjArray->Add(alignObj);
456   }
457   // delete entry --- Don't delete, it is cached!
458         
459   AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
460   return kTRUE;
461
462 }
463
464 //_____________________________________________________________________________
465 Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
466 {
467   // Read the alignment objects from CDB.
468   // Each detector is supposed to have the
469   // alignment objects in DET/Align/Data CDB path.
470   // All the detector objects are then collected,
471   // sorted by geometry level (starting from ALIC) and
472   // then applied to the TGeo geometry.
473   // Finally an overlaps check is performed.
474
475   // Load alignment data from CDB and fill fAlignObjArray 
476   if(fLoadAlignFromCDB){
477         if(!fAlignObjArray) fAlignObjArray = new TObjArray();
478         
479         //fAlignObjArray->RemoveAll(); 
480         fAlignObjArray->Clear();        
481         fAlignObjArray->SetOwner(0);
482  
483         TString detStr = detectors;
484         TString dataNotLoaded="";
485         TString dataLoaded="";
486   
487         for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
488           if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
489           if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
490             dataNotLoaded += fgkDetectorName[iDet];
491             dataNotLoaded += " ";
492           } else {
493             dataLoaded += fgkDetectorName[iDet];
494             dataLoaded += " ";
495           }
496         } // end loop over detectors
497   
498         if ((detStr.CompareTo("ALL") == 0)) detStr = "";
499         dataNotLoaded += detStr;
500         if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
501                           dataLoaded.Data()));
502         if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
503                           dataNotLoaded.Data()));
504   } // fLoadAlignFromCDB flag
505  
506   // Check if the array with alignment objects was
507   // provided by the user. If yes, apply the objects
508   // to the present TGeo geometry
509   if (fAlignObjArray) {
510     if (gGeoManager && gGeoManager->IsClosed()) {
511       if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
512         AliError("The misalignment of one or more volumes failed!"
513                  "Compare the list of simulated detectors and the list of detector alignment data!");
514         return kFALSE;
515       }
516     }
517     else {
518       AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
519       return kFALSE;
520     }
521   }
522
523   delete fAlignObjArray; fAlignObjArray=0;
524   return kTRUE;
525 }
526
527 //_____________________________________________________________________________
528 void AliReconstruction::SetGAliceFile(const char* fileName)
529 {
530 // set the name of the galice file
531
532   fGAliceFileName = fileName;
533 }
534
535 //_____________________________________________________________________________
536 void AliReconstruction::SetOption(const char* detector, const char* option)
537 {
538 // set options for the reconstruction of a detector
539
540   TObject* obj = fOptions.FindObject(detector);
541   if (obj) fOptions.Remove(obj);
542   fOptions.Add(new TNamed(detector, option));
543 }
544
545
546 //_____________________________________________________________________________
547 Bool_t AliReconstruction::Run(const char* input)
548 {
549 // run the reconstruction
550
551   // set the input
552   if (!input) input = fInput.Data();
553   TString fileName(input);
554   if (fileName.EndsWith("/")) {
555     fRawReader = new AliRawReaderFile(fileName);
556   } else if (fileName.EndsWith(".root")) {
557     fRawReader = new AliRawReaderRoot(fileName);
558   } else if (!fileName.IsNull()) {
559     fRawReader = new AliRawReaderDate(fileName);
560     fRawReader->SelectEvents(7);
561   }
562   if (!fEquipIdMap.IsNull() && fRawReader)
563     fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
564
565
566   // get the run loader
567   if (!InitRunLoader()) return kFALSE;
568
569   // Initialize the CDB storage
570   InitCDBStorage();
571
572   // Set run number in CDBManager (if it is not already set by the user)
573   if (!SetRunNumber()) if (fStopOnError) return kFALSE;
574
575   // Import ideal TGeo geometry and apply misalignment
576   if (!gGeoManager) {
577     TString geom(gSystem->DirName(fGAliceFileName));
578     geom += "/geometry.root";
579     TGeoManager::Import(geom.Data());
580     if (!gGeoManager) if (fStopOnError) return kFALSE;
581   }
582
583   AliCDBManager* man = AliCDBManager::Instance();
584   if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
585
586   // local reconstruction
587   if (!fRunLocalReconstruction.IsNull()) {
588     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
589       if (fStopOnError) {CleanUp(); return kFALSE;}
590     }
591   }
592 //  if (!fRunVertexFinder && fRunTracking.IsNull() && 
593 //      fFillESD.IsNull()) return kTRUE;
594
595   // get vertexer
596   if (fRunVertexFinder && !CreateVertexer()) {
597     if (fStopOnError) {
598       CleanUp(); 
599       return kFALSE;
600     }
601   }
602
603   // get trackers
604   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
605     if (fStopOnError) {
606       CleanUp(); 
607       return kFALSE;
608     }      
609   }
610
611
612   TStopwatch stopwatch;
613   stopwatch.Start();
614
615   // get the possibly already existing ESD file and tree
616   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
617   TFile* fileOld = NULL;
618   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
619   if (!gSystem->AccessPathName("AliESDs.root")){
620     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
621     fileOld = TFile::Open("AliESDs.old.root");
622     if (fileOld && fileOld->IsOpen()) {
623       treeOld = (TTree*) fileOld->Get("esdTree");
624       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
625       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
626       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
627     }
628   }
629
630   // create the ESD output file and tree
631   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
632   if (!file->IsOpen()) {
633     AliError("opening AliESDs.root failed");
634     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
635   }
636   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
637   tree->Branch("ESD", "AliESD", &esd);
638   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
639   hlttree->Branch("ESD", "AliESD", &hltesd);
640   delete esd; delete hltesd;
641   esd = NULL; hltesd = NULL;
642
643   // create the branch with ESD additions
644   AliESDfriend *esdf=0;
645   if (fWriteESDfriend) {
646      TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
647      br->SetFile("AliESDfriends.root");
648   }
649
650   AliVertexerTracks tVertexer;
651   if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
652
653   // loop over events
654   if (fRawReader) fRawReader->RewindEvents();
655   
656   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
657     if (fRawReader) fRawReader->NextEvent();
658     if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
659       // copy old ESD to the new one
660       if (treeOld) {
661         treeOld->SetBranchAddress("ESD", &esd);
662         treeOld->GetEntry(iEvent);
663       }
664       tree->Fill();
665       if (hlttreeOld) {
666         hlttreeOld->SetBranchAddress("ESD", &hltesd);
667         hlttreeOld->GetEntry(iEvent);
668       }
669       hlttree->Fill();
670       continue;
671     }
672
673     AliInfo(Form("processing event %d", iEvent));
674     fRunLoader->GetEvent(iEvent);
675
676     char fileName[256];
677     sprintf(fileName, "ESD_%d.%d_final.root", 
678             fRunLoader->GetHeader()->GetRun(), 
679             fRunLoader->GetHeader()->GetEventNrInRun());
680     if (!gSystem->AccessPathName(fileName)) continue;
681
682     // local reconstruction
683     if (!fRunLocalReconstruction.IsNull()) {
684       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
685         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
686       }
687     }
688
689     esd = new AliESD; hltesd = new AliESD;
690     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
691     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
692     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
693     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
694
695     // Set magnetic field from the tracker
696     esd->SetMagneticField(AliTracker::GetBz());
697     hltesd->SetMagneticField(AliTracker::GetBz());
698
699     // vertex finder
700     if (fRunVertexFinder) {
701       if (!ReadESD(esd, "vertex")) {
702         if (!RunVertexFinder(esd)) {
703           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
704         }
705         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
706       }
707     }
708
709     // HLT tracking
710     if (!fRunTracking.IsNull()) {
711       if (fRunHLTTracking) {
712         hltesd->SetVertex(esd->GetVertex());
713         if (!RunHLTTracking(hltesd)) {
714           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
715         }
716       }
717     }
718
719     // barrel tracking
720     if (!fRunTracking.IsNull()) {
721       if (!ReadESD(esd, "tracking")) {
722         if (!RunTracking(esd)) {
723           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
724         }
725         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
726       }
727     }
728
729     // fill ESD
730     if (!fFillESD.IsNull()) {
731       if (!FillESD(esd, fFillESD)) {
732         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
733       }
734     }
735     // fill Event header information from the RawEventHeader
736     if (fRawReader){FillRawEventHeaderESD(esd);}
737
738     // combined PID
739     AliESDpid::MakePID(esd);
740     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
741
742     if (fFillTriggerESD) {
743       if (!ReadESD(esd, "trigger")) {
744         if (!FillTriggerESD(esd)) {
745           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
746         }
747         if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
748       }
749     }
750
751     esd->SetPrimaryVertex(tVertexer.FindPrimaryVertex(esd));
752
753     {
754     // V0 finding
755     AliV0vertexer vtxer;
756     vtxer.Tracks2V0vertices(esd);
757
758     // Cascade finding
759     AliCascadeVertexer cvtxer;
760     cvtxer.V0sTracks2CascadeVertices(esd);
761     }
762  
763     // write ESD
764     if (fWriteESDfriend) {
765        esdf=new AliESDfriend();
766        esd->GetESDfriend(esdf);
767     }
768     tree->Fill();
769
770     // write HLT ESD
771     hlttree->Fill();
772
773     if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
774  
775     delete esd; delete esdf; delete hltesd;
776     esd = NULL; esdf=NULL; hltesd = NULL;
777   }
778
779   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
780                stopwatch.RealTime(),stopwatch.CpuTime()));
781
782   file->cd();
783   if (fWriteESDfriend)
784     tree->SetBranchStatus("ESDfriend*",0);
785   tree->Write();
786   hlttree->Write();
787
788   // Create tags for the events in the ESD tree (the ESD tree is always present)
789   // In case of empty events the tags will contain dummy values
790   CreateTag(file);
791   CleanUp(file, fileOld);
792
793   return kTRUE;
794 }
795
796
797 //_____________________________________________________________________________
798 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
799 {
800 // run the local reconstruction
801
802   TStopwatch stopwatch;
803   stopwatch.Start();
804
805   AliCDBManager* man = AliCDBManager::Instance();
806   Bool_t origCache = man->GetCacheFlag();
807
808   TString detStr = detectors;
809   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
810     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
811     AliReconstructor* reconstructor = GetReconstructor(iDet);
812     if (!reconstructor) continue;
813     if (reconstructor->HasLocalReconstruction()) continue;
814
815     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
816     TStopwatch stopwatchDet;
817     stopwatchDet.Start();
818
819     AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
820
821     man->SetCacheFlag(kTRUE);
822     TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
823     man->GetAll(calibPath); // entries are cached!
824
825     if (fRawReader) {
826       fRawReader->RewindEvents();
827       reconstructor->Reconstruct(fRunLoader, fRawReader);
828     } else {
829       reconstructor->Reconstruct(fRunLoader);
830     }
831     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
832                  fgkDetectorName[iDet],
833                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
834
835     // unload calibration data
836     man->ClearCache();
837   }
838
839   man->SetCacheFlag(origCache);
840
841   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
842     AliError(Form("the following detectors were not found: %s",
843                   detStr.Data()));
844     if (fStopOnError) return kFALSE;
845   }
846
847   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
848                stopwatch.RealTime(),stopwatch.CpuTime()));
849
850   return kTRUE;
851 }
852
853 //_____________________________________________________________________________
854 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
855 {
856 // run the local reconstruction
857
858   TStopwatch stopwatch;
859   stopwatch.Start();
860
861   TString detStr = detectors;
862   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
863     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
864     AliReconstructor* reconstructor = GetReconstructor(iDet);
865     if (!reconstructor) continue;
866     AliLoader* loader = fLoader[iDet];
867
868     // conversion of digits
869     if (fRawReader && reconstructor->HasDigitConversion()) {
870       AliInfo(Form("converting raw data digits into root objects for %s", 
871                    fgkDetectorName[iDet]));
872       TStopwatch stopwatchDet;
873       stopwatchDet.Start();
874       loader->LoadDigits("update");
875       loader->CleanDigits();
876       loader->MakeDigitsContainer();
877       TTree* digitsTree = loader->TreeD();
878       reconstructor->ConvertDigits(fRawReader, digitsTree);
879       loader->WriteDigits("OVERWRITE");
880       loader->UnloadDigits();
881       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
882                    fgkDetectorName[iDet],
883                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
884     }
885
886     // local reconstruction
887     if (!reconstructor->HasLocalReconstruction()) continue;
888     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
889     TStopwatch stopwatchDet;
890     stopwatchDet.Start();
891     loader->LoadRecPoints("update");
892     loader->CleanRecPoints();
893     loader->MakeRecPointsContainer();
894     TTree* clustersTree = loader->TreeR();
895     if (fRawReader && !reconstructor->HasDigitConversion()) {
896       reconstructor->Reconstruct(fRawReader, clustersTree);
897     } else {
898       loader->LoadDigits("read");
899       TTree* digitsTree = loader->TreeD();
900       if (!digitsTree) {
901         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
902         if (fStopOnError) return kFALSE;
903       } else {
904         reconstructor->Reconstruct(digitsTree, clustersTree);
905       }
906       loader->UnloadDigits();
907     }
908     loader->WriteRecPoints("OVERWRITE");
909     loader->UnloadRecPoints();
910     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
911                     fgkDetectorName[iDet],
912                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
913   }
914
915   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
916     AliError(Form("the following detectors were not found: %s",
917                   detStr.Data()));
918     if (fStopOnError) return kFALSE;
919   }
920   
921   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
922                stopwatch.RealTime(),stopwatch.CpuTime()));
923
924   return kTRUE;
925 }
926
927 //_____________________________________________________________________________
928 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
929 {
930 // run the barrel tracking
931
932   TStopwatch stopwatch;
933   stopwatch.Start();
934
935   AliESDVertex* vertex = NULL;
936   Double_t vtxPos[3] = {0, 0, 0};
937   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
938   TArrayF mcVertex(3); 
939   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
940     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
941     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
942   }
943
944   if (fVertexer) {
945     AliInfo("running the ITS vertex finder");
946     if (fLoader[0]) fLoader[0]->LoadRecPoints();
947     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
948     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
949     if(!vertex){
950       AliWarning("Vertex not found");
951       vertex = new AliESDVertex();
952       vertex->SetName("default");
953     }
954     else {
955       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
956       vertex->SetName("reconstructed");
957     }
958
959   } else {
960     AliInfo("getting the primary vertex from MC");
961     vertex = new AliESDVertex(vtxPos, vtxErr);
962   }
963
964   if (vertex) {
965     vertex->GetXYZ(vtxPos);
966     vertex->GetSigmaXYZ(vtxErr);
967   } else {
968     AliWarning("no vertex reconstructed");
969     vertex = new AliESDVertex(vtxPos, vtxErr);
970   }
971   esd->SetVertex(vertex);
972   // if SPD multiplicity has been determined, it is stored in the ESD
973   AliMultiplicity *mult= fVertexer->GetMultiplicity();
974   if(mult)esd->SetMultiplicity(mult);
975
976   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
977     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
978   }  
979   delete vertex;
980
981   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
982                stopwatch.RealTime(),stopwatch.CpuTime()));
983
984   return kTRUE;
985 }
986
987 //_____________________________________________________________________________
988 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
989 {
990 // run the HLT barrel tracking
991
992   TStopwatch stopwatch;
993   stopwatch.Start();
994
995   if (!fRunLoader) {
996     AliError("Missing runLoader!");
997     return kFALSE;
998   }
999
1000   AliInfo("running HLT tracking");
1001
1002   // Get a pointer to the HLT reconstructor
1003   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1004   if (!reconstructor) return kFALSE;
1005
1006   // TPC + ITS
1007   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1008     TString detName = fgkDetectorName[iDet];
1009     AliDebug(1, Form("%s HLT tracking", detName.Data()));
1010     reconstructor->SetOption(detName.Data());
1011     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1012     if (!tracker) {
1013       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1014       if (fStopOnError) return kFALSE;
1015       continue;
1016     }
1017     Double_t vtxPos[3];
1018     Double_t vtxErr[3]={0.005,0.005,0.010};
1019     const AliESDVertex *vertex = esd->GetVertex();
1020     vertex->GetXYZ(vtxPos);
1021     tracker->SetVertex(vtxPos,vtxErr);
1022     if(iDet != 1) {
1023       fLoader[iDet]->LoadRecPoints("read");
1024       TTree* tree = fLoader[iDet]->TreeR();
1025       if (!tree) {
1026         AliError(Form("Can't get the %s cluster tree", detName.Data()));
1027         return kFALSE;
1028       }
1029       tracker->LoadClusters(tree);
1030     }
1031     if (tracker->Clusters2Tracks(esd) != 0) {
1032       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1033       return kFALSE;
1034     }
1035     if(iDet != 1) {
1036       tracker->UnloadClusters();
1037     }
1038     delete tracker;
1039   }
1040
1041   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1042                stopwatch.RealTime(),stopwatch.CpuTime()));
1043
1044   return kTRUE;
1045 }
1046
1047 //_____________________________________________________________________________
1048 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1049 {
1050 // run the barrel tracking
1051
1052   TStopwatch stopwatch;
1053   stopwatch.Start();
1054
1055   AliInfo("running tracking");
1056
1057   // pass 1: TPC + ITS inwards
1058   for (Int_t iDet = 1; iDet >= 0; iDet--) {
1059     if (!fTracker[iDet]) continue;
1060     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1061
1062     // load clusters
1063     fLoader[iDet]->LoadRecPoints("read");
1064     TTree* tree = fLoader[iDet]->TreeR();
1065     if (!tree) {
1066       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1067       return kFALSE;
1068     }
1069     fTracker[iDet]->LoadClusters(tree);
1070
1071     // run tracking
1072     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1073       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1074       return kFALSE;
1075     }
1076     if (fCheckPointLevel > 1) {
1077       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1078     }
1079     // preliminary PID in TPC needed by the ITS tracker
1080     if (iDet == 1) {
1081       GetReconstructor(1)->FillESD(fRunLoader, esd);
1082       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1083       AliESDpid::MakePID(esd);
1084     }
1085   }
1086
1087   // pass 2: ALL backwards
1088   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1089     if (!fTracker[iDet]) continue;
1090     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1091
1092     // load clusters
1093     if (iDet > 1) {     // all except ITS, TPC
1094       TTree* tree = NULL;
1095       fLoader[iDet]->LoadRecPoints("read");
1096       tree = fLoader[iDet]->TreeR();
1097       if (!tree) {
1098         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1099         return kFALSE;
1100       }
1101       fTracker[iDet]->LoadClusters(tree);
1102     }
1103
1104     // run tracking
1105     if (fTracker[iDet]->PropagateBack(esd) != 0) {
1106       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1107       return kFALSE;
1108     }
1109     if (fCheckPointLevel > 1) {
1110       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1111     }
1112
1113     // unload clusters
1114     if (iDet > 2) {     // all except ITS, TPC, TRD
1115       fTracker[iDet]->UnloadClusters();
1116       fLoader[iDet]->UnloadRecPoints();
1117     }
1118     // updated PID in TPC needed by the ITS tracker -MI
1119     if (iDet == 1) {
1120       GetReconstructor(1)->FillESD(fRunLoader, esd);
1121       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1122       AliESDpid::MakePID(esd);
1123     }
1124   }
1125
1126   // write space-points to the ESD in case alignment data output
1127   // is switched on
1128   if (fWriteAlignmentData)
1129     WriteAlignmentData(esd);
1130
1131   // pass 3: TRD + TPC + ITS refit inwards
1132   for (Int_t iDet = 2; iDet >= 0; iDet--) {
1133     if (!fTracker[iDet]) continue;
1134     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1135
1136     // run tracking
1137     if (fTracker[iDet]->RefitInward(esd) != 0) {
1138       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1139       return kFALSE;
1140     }
1141     if (fCheckPointLevel > 1) {
1142       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1143     }
1144
1145     // unload clusters
1146     fTracker[iDet]->UnloadClusters();
1147     fLoader[iDet]->UnloadRecPoints();
1148   }
1149   //
1150   // Propagate track to the vertex - if not done by ITS
1151   //
1152   Int_t ntracks = esd->GetNumberOfTracks();
1153   for (Int_t itrack=0; itrack<ntracks; itrack++){
1154     const Double_t kRadius  = 3;   // beam pipe radius
1155     const Double_t kMaxStep = 5;   // max step
1156     const Double_t kMaxD    = 123456;  // max distance to prim vertex
1157     Double_t       fieldZ   = AliTracker::GetBz();  //
1158     AliESDtrack * track = esd->GetTrack(itrack);
1159     if (!track) continue;
1160     if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1161     track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1162     track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1163   }
1164   
1165   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1166                stopwatch.RealTime(),stopwatch.CpuTime()));
1167
1168   return kTRUE;
1169 }
1170
1171 //_____________________________________________________________________________
1172 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1173 {
1174 // fill the event summary data
1175
1176   TStopwatch stopwatch;
1177   stopwatch.Start();
1178   AliInfo("filling ESD");
1179
1180   TString detStr = detectors;
1181   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1182     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1183     AliReconstructor* reconstructor = GetReconstructor(iDet);
1184     if (!reconstructor) continue;
1185
1186     if (!ReadESD(esd, fgkDetectorName[iDet])) {
1187       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1188       TTree* clustersTree = NULL;
1189       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1190         fLoader[iDet]->LoadRecPoints("read");
1191         clustersTree = fLoader[iDet]->TreeR();
1192         if (!clustersTree) {
1193           AliError(Form("Can't get the %s clusters tree", 
1194                         fgkDetectorName[iDet]));
1195           if (fStopOnError) return kFALSE;
1196         }
1197       }
1198       if (fRawReader && !reconstructor->HasDigitConversion()) {
1199         reconstructor->FillESD(fRawReader, clustersTree, esd);
1200       } else {
1201         TTree* digitsTree = NULL;
1202         if (fLoader[iDet]) {
1203           fLoader[iDet]->LoadDigits("read");
1204           digitsTree = fLoader[iDet]->TreeD();
1205           if (!digitsTree) {
1206             AliError(Form("Can't get the %s digits tree", 
1207                           fgkDetectorName[iDet]));
1208             if (fStopOnError) return kFALSE;
1209           }
1210         }
1211         reconstructor->FillESD(digitsTree, clustersTree, esd);
1212         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1213       }
1214       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1215         fLoader[iDet]->UnloadRecPoints();
1216       }
1217
1218       if (fRawReader) {
1219         reconstructor->FillESD(fRunLoader, fRawReader, esd);
1220       } else {
1221         reconstructor->FillESD(fRunLoader, esd);
1222       }
1223       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1224     }
1225   }
1226
1227   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1228     AliError(Form("the following detectors were not found: %s", 
1229                   detStr.Data()));
1230     if (fStopOnError) return kFALSE;
1231   }
1232
1233   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1234                stopwatch.RealTime(),stopwatch.CpuTime()));
1235
1236   return kTRUE;
1237 }
1238
1239 //_____________________________________________________________________________
1240 Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1241 {
1242   // Reads the trigger decision which is
1243   // stored in Trigger.root file and fills
1244   // the corresponding esd entries
1245
1246   AliInfo("Filling trigger information into the ESD");
1247
1248   if (fRawReader) {
1249     AliCTPRawStream input(fRawReader);
1250     if (!input.Next()) {
1251       AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1252       return kFALSE;
1253     }
1254     esd->SetTriggerMask(input.GetClassMask());
1255     esd->SetTriggerCluster(input.GetClusterMask());
1256   }
1257   else {
1258     AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1259     if (runloader) {
1260       if (!runloader->LoadTrigger()) {
1261         AliCentralTrigger *aCTP = runloader->GetTrigger();
1262         esd->SetTriggerMask(aCTP->GetClassMask());
1263         esd->SetTriggerCluster(aCTP->GetClusterMask());
1264       }
1265       else {
1266         AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1267         return kFALSE;
1268       }
1269     }
1270     else {
1271       AliError("No run loader is available! The trigger information is not stored in the ESD !");
1272       return kFALSE;
1273     }
1274   }
1275
1276   return kTRUE;
1277 }
1278
1279
1280
1281
1282
1283 //_____________________________________________________________________________
1284 Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1285 {
1286   // 
1287   // Filling information from RawReader Header
1288   // 
1289
1290   AliInfo("Filling information from RawReader Header");
1291   esd->SetTimeStamp(0);
1292   esd->SetEventType(0);
1293   const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1294   if (eventHeader){
1295     esd->SetTimeStamp((eventHeader->Get("Timestamp")));  
1296     esd->SetEventType((eventHeader->Get("Type")));  
1297   }
1298
1299   return kTRUE;
1300 }
1301
1302
1303 //_____________________________________________________________________________
1304 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1305 {
1306 // check whether detName is contained in detectors
1307 // if yes, it is removed from detectors
1308
1309   // check if all detectors are selected
1310   if ((detectors.CompareTo("ALL") == 0) ||
1311       detectors.BeginsWith("ALL ") ||
1312       detectors.EndsWith(" ALL") ||
1313       detectors.Contains(" ALL ")) {
1314     detectors = "ALL";
1315     return kTRUE;
1316   }
1317
1318   // search for the given detector
1319   Bool_t result = kFALSE;
1320   if ((detectors.CompareTo(detName) == 0) ||
1321       detectors.BeginsWith(detName+" ") ||
1322       detectors.EndsWith(" "+detName) ||
1323       detectors.Contains(" "+detName+" ")) {
1324     detectors.ReplaceAll(detName, "");
1325     result = kTRUE;
1326   }
1327
1328   // clean up the detectors string
1329   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
1330   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1331   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1332
1333   return result;
1334 }
1335
1336 //_____________________________________________________________________________
1337 Bool_t AliReconstruction::InitRunLoader()
1338 {
1339 // get or create the run loader
1340
1341   if (gAlice) delete gAlice;
1342   gAlice = NULL;
1343
1344   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1345     // load all base libraries to get the loader classes
1346     TString libs = gSystem->GetLibraries();
1347     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1348       TString detName = fgkDetectorName[iDet];
1349       if (detName == "HLT") continue;
1350       if (libs.Contains("lib" + detName + "base.so")) continue;
1351       gSystem->Load("lib" + detName + "base.so");
1352     }
1353     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1354     if (!fRunLoader) {
1355       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1356       CleanUp();
1357       return kFALSE;
1358     }
1359     fRunLoader->CdGAFile();
1360     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1361       if (fRunLoader->LoadgAlice() == 0) {
1362         gAlice = fRunLoader->GetAliRun();
1363         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1364       }
1365     }
1366     if (!gAlice && !fRawReader) {
1367       AliError(Form("no gAlice object found in file %s",
1368                     fGAliceFileName.Data()));
1369       CleanUp();
1370       return kFALSE;
1371     }
1372
1373   } else {               // galice.root does not exist
1374     if (!fRawReader) {
1375       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1376       CleanUp();
1377       return kFALSE;
1378     }
1379     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1380                                     AliConfig::GetDefaultEventFolderName(),
1381                                     "recreate");
1382     if (!fRunLoader) {
1383       AliError(Form("could not create run loader in file %s", 
1384                     fGAliceFileName.Data()));
1385       CleanUp();
1386       return kFALSE;
1387     }
1388     fRunLoader->MakeTree("E");
1389     Int_t iEvent = 0;
1390     while (fRawReader->NextEvent()) {
1391       fRunLoader->SetEventNumber(iEvent);
1392       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
1393                                      iEvent, iEvent);
1394       fRunLoader->MakeTree("H");
1395       fRunLoader->TreeE()->Fill();
1396       iEvent++;
1397     }
1398     fRawReader->RewindEvents();
1399     fRunLoader->WriteHeader("OVERWRITE");
1400     fRunLoader->CdGAFile();
1401     fRunLoader->Write(0, TObject::kOverwrite);
1402 //    AliTracker::SetFieldMap(???);
1403   }
1404
1405   return kTRUE;
1406 }
1407
1408 //_____________________________________________________________________________
1409 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1410 {
1411 // get the reconstructor object and the loader for a detector
1412
1413   if (fReconstructor[iDet]) return fReconstructor[iDet];
1414
1415   // load the reconstructor object
1416   TPluginManager* pluginManager = gROOT->GetPluginManager();
1417   TString detName = fgkDetectorName[iDet];
1418   TString recName = "Ali" + detName + "Reconstructor";
1419   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1420
1421   if (detName == "HLT") {
1422     if (!gROOT->GetClass("AliLevel3")) {
1423       gSystem->Load("libAliHLTSrc.so");
1424       gSystem->Load("libAliHLTMisc.so");
1425       gSystem->Load("libAliHLTHough.so");
1426       gSystem->Load("libAliHLTComp.so");
1427     }
1428   }
1429
1430   AliReconstructor* reconstructor = NULL;
1431   // first check if a plugin is defined for the reconstructor
1432   TPluginHandler* pluginHandler = 
1433     pluginManager->FindHandler("AliReconstructor", detName);
1434   // if not, add a plugin for it
1435   if (!pluginHandler) {
1436     AliDebug(1, Form("defining plugin for %s", recName.Data()));
1437     TString libs = gSystem->GetLibraries();
1438     if (libs.Contains("lib" + detName + "base.so") ||
1439         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1440       pluginManager->AddHandler("AliReconstructor", detName, 
1441                                 recName, detName + "rec", recName + "()");
1442     } else {
1443       pluginManager->AddHandler("AliReconstructor", detName, 
1444                                 recName, detName, recName + "()");
1445     }
1446     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1447   }
1448   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1449     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1450   }
1451   if (reconstructor) {
1452     TObject* obj = fOptions.FindObject(detName.Data());
1453     if (obj) reconstructor->SetOption(obj->GetTitle());
1454     reconstructor->Init(fRunLoader);
1455     fReconstructor[iDet] = reconstructor;
1456   }
1457
1458   // get or create the loader
1459   if (detName != "HLT") {
1460     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1461     if (!fLoader[iDet]) {
1462       AliConfig::Instance()
1463         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
1464                                 detName, detName);
1465       // first check if a plugin is defined for the loader
1466       TPluginHandler* pluginHandler = 
1467         pluginManager->FindHandler("AliLoader", detName);
1468       // if not, add a plugin for it
1469       if (!pluginHandler) {
1470         TString loaderName = "Ali" + detName + "Loader";
1471         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1472         pluginManager->AddHandler("AliLoader", detName, 
1473                                   loaderName, detName + "base", 
1474                                   loaderName + "(const char*, TFolder*)");
1475         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1476       }
1477       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1478         fLoader[iDet] = 
1479           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1480                                                  fRunLoader->GetEventFolder());
1481       }
1482       if (!fLoader[iDet]) {   // use default loader
1483         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1484       }
1485       if (!fLoader[iDet]) {
1486         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1487         if (fStopOnError) return NULL;
1488       } else {
1489         fRunLoader->AddLoader(fLoader[iDet]);
1490         fRunLoader->CdGAFile();
1491         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1492         fRunLoader->Write(0, TObject::kOverwrite);
1493       }
1494     }
1495   }
1496       
1497   return reconstructor;
1498 }
1499
1500 //_____________________________________________________________________________
1501 Bool_t AliReconstruction::CreateVertexer()
1502 {
1503 // create the vertexer
1504
1505   fVertexer = NULL;
1506   AliReconstructor* itsReconstructor = GetReconstructor(0);
1507   if (itsReconstructor) {
1508     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1509   }
1510   if (!fVertexer) {
1511     AliWarning("couldn't create a vertexer for ITS");
1512     if (fStopOnError) return kFALSE;
1513   }
1514
1515   return kTRUE;
1516 }
1517
1518 //_____________________________________________________________________________
1519 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1520 {
1521 // create the trackers
1522
1523   TString detStr = detectors;
1524   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1525     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1526     AliReconstructor* reconstructor = GetReconstructor(iDet);
1527     if (!reconstructor) continue;
1528     TString detName = fgkDetectorName[iDet];
1529     if (detName == "HLT") {
1530       fRunHLTTracking = kTRUE;
1531       continue;
1532     }
1533
1534     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1535     if (!fTracker[iDet] && (iDet < 7)) {
1536       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1537       if (fStopOnError) return kFALSE;
1538     }
1539   }
1540
1541   return kTRUE;
1542 }
1543
1544 //_____________________________________________________________________________
1545 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1546 {
1547 // delete trackers and the run loader and close and delete the file
1548
1549   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1550     delete fReconstructor[iDet];
1551     fReconstructor[iDet] = NULL;
1552     fLoader[iDet] = NULL;
1553     delete fTracker[iDet];
1554     fTracker[iDet] = NULL;
1555   }
1556   delete fVertexer;
1557   fVertexer = NULL;
1558   delete fDiamondProfile;
1559   fDiamondProfile = NULL;
1560
1561   delete fRunLoader;
1562   fRunLoader = NULL;
1563   delete fRawReader;
1564   fRawReader = NULL;
1565
1566   if (file) {
1567     file->Close();
1568     delete file;
1569   }
1570
1571   if (fileOld) {
1572     fileOld->Close();
1573     delete fileOld;
1574     gSystem->Unlink("AliESDs.old.root");
1575   }
1576 }
1577
1578
1579 //_____________________________________________________________________________
1580 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1581 {
1582 // read the ESD event from a file
1583
1584   if (!esd) return kFALSE;
1585   char fileName[256];
1586   sprintf(fileName, "ESD_%d.%d_%s.root", 
1587           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1588   if (gSystem->AccessPathName(fileName)) return kFALSE;
1589
1590   AliInfo(Form("reading ESD from file %s", fileName));
1591   AliDebug(1, Form("reading ESD from file %s", fileName));
1592   TFile* file = TFile::Open(fileName);
1593   if (!file || !file->IsOpen()) {
1594     AliError(Form("opening %s failed", fileName));
1595     delete file;
1596     return kFALSE;
1597   }
1598
1599   gROOT->cd();
1600   delete esd;
1601   esd = (AliESD*) file->Get("ESD");
1602   file->Close();
1603   delete file;
1604   return kTRUE;
1605 }
1606
1607 //_____________________________________________________________________________
1608 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1609 {
1610 // write the ESD event to a file
1611
1612   if (!esd) return;
1613   char fileName[256];
1614   sprintf(fileName, "ESD_%d.%d_%s.root", 
1615           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1616
1617   AliDebug(1, Form("writing ESD to file %s", fileName));
1618   TFile* file = TFile::Open(fileName, "recreate");
1619   if (!file || !file->IsOpen()) {
1620     AliError(Form("opening %s failed", fileName));
1621   } else {
1622     esd->Write("ESD");
1623     file->Close();
1624   }
1625   delete file;
1626 }
1627
1628
1629
1630
1631 //_____________________________________________________________________________
1632 void AliReconstruction::CreateTag(TFile* file)
1633 {
1634   /////////////
1635   //muon code//
1636   ////////////
1637   Double_t fMUONMASS = 0.105658369;
1638   //Variables
1639   Double_t fX,fY,fZ ;
1640   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1641   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1642   Int_t fCharge;
1643   TLorentzVector fEPvector;
1644
1645   Float_t fZVertexCut = 10.0; 
1646   Float_t fRhoVertexCut = 2.0; 
1647
1648   Float_t fLowPtCut = 1.0;
1649   Float_t fHighPtCut = 3.0;
1650   Float_t fVeryHighPtCut = 10.0;
1651   ////////////
1652
1653   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1654
1655   // Creates the tags for all the events in a given ESD file
1656   Int_t ntrack;
1657   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1658   Int_t nPos, nNeg, nNeutr;
1659   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1660   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1661   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1662   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1663   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1664   Int_t fVertexflag;
1665   Int_t iRunNumber = 0;
1666   TString fVertexName("default");
1667
1668   AliRunTag *tag = new AliRunTag();
1669   AliEventTag *evTag = new AliEventTag();
1670   TTree ttag("T","A Tree with event tags");
1671   TBranch * btag = ttag.Branch("AliTAG", &tag);
1672   btag->SetCompressionLevel(9);
1673   
1674   AliInfo(Form("Creating the tags......."));    
1675   
1676   if (!file || !file->IsOpen()) {
1677     AliError(Form("opening failed"));
1678     delete file;
1679     return ;
1680   }  
1681   Int_t lastEvent = 0;
1682   TTree *t = (TTree*) file->Get("esdTree");
1683   TBranch * b = t->GetBranch("ESD");
1684   AliESD *esd = 0;
1685   b->SetAddress(&esd);
1686
1687   b->GetEntry(fFirstEvent);
1688   Int_t iInitRunNumber = esd->GetRunNumber();
1689   
1690   Int_t iNumberOfEvents = b->GetEntries();
1691   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1692   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1693     ntrack = 0;
1694     nPos = 0;
1695     nNeg = 0;
1696     nNeutr =0;
1697     nK0s = 0;
1698     nNeutrons = 0;
1699     nPi0s = 0;
1700     nGamas = 0;
1701     nProtons = 0;
1702     nKaons = 0;
1703     nPions = 0;
1704     nMuons = 0;
1705     nElectrons = 0;       
1706     nCh1GeV = 0;
1707     nCh3GeV = 0;
1708     nCh10GeV = 0;
1709     nMu1GeV = 0;
1710     nMu3GeV = 0;
1711     nMu10GeV = 0;
1712     nEl1GeV = 0;
1713     nEl3GeV = 0;
1714     nEl10GeV = 0;
1715     maxPt = .0;
1716     meanPt = .0;
1717     totalP = .0;
1718     fVertexflag = 0;
1719
1720     b->GetEntry(iEventNumber);
1721     iRunNumber = esd->GetRunNumber();
1722     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1723     const AliESDVertex * vertexIn = esd->GetVertex();
1724     if (!vertexIn) AliError("ESD has not defined vertex.");
1725     if (vertexIn) fVertexName = vertexIn->GetName();
1726     if(fVertexName != "default") fVertexflag = 1;
1727     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1728       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1729       UInt_t status = esdTrack->GetStatus();
1730       
1731       //select only tracks with ITS refit
1732       if ((status&AliESDtrack::kITSrefit)==0) continue;
1733       //select only tracks with TPC refit
1734       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1735       
1736       //select only tracks with the "combined PID"
1737       if ((status&AliESDtrack::kESDpid)==0) continue;
1738       Double_t p[3];
1739       esdTrack->GetPxPyPz(p);
1740       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1741       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1742       totalP += momentum;
1743       meanPt += fPt;
1744       if(fPt > maxPt) maxPt = fPt;
1745       
1746       if(esdTrack->GetSign() > 0) {
1747         nPos++;
1748         if(fPt > fLowPtCut) nCh1GeV++;
1749         if(fPt > fHighPtCut) nCh3GeV++;
1750         if(fPt > fVeryHighPtCut) nCh10GeV++;
1751       }
1752       if(esdTrack->GetSign() < 0) {
1753         nNeg++;
1754         if(fPt > fLowPtCut) nCh1GeV++;
1755         if(fPt > fHighPtCut) nCh3GeV++;
1756         if(fPt > fVeryHighPtCut) nCh10GeV++;
1757       }
1758       if(esdTrack->GetSign() == 0) nNeutr++;
1759       
1760       //PID
1761       Double_t prob[5];
1762       esdTrack->GetESDpid(prob);
1763       
1764       Double_t rcc = 0.0;
1765       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1766       if(rcc == 0.0) continue;
1767       //Bayes' formula
1768       Double_t w[5];
1769       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1770       
1771       //protons
1772       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1773       //kaons
1774       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1775       //pions
1776       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
1777       //electrons
1778       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1779         nElectrons++;
1780         if(fPt > fLowPtCut) nEl1GeV++;
1781         if(fPt > fHighPtCut) nEl3GeV++;
1782         if(fPt > fVeryHighPtCut) nEl10GeV++;
1783       }   
1784       ntrack++;
1785     }//track loop
1786     
1787     /////////////
1788     //muon code//
1789     ////////////
1790     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1791     // loop over all reconstructed tracks (also first track of combination)
1792     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1793       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1794       if (muonTrack == 0x0) continue;
1795       
1796       // Coordinates at vertex
1797       fZ = muonTrack->GetZ(); 
1798       fY = muonTrack->GetBendingCoor();
1799       fX = muonTrack->GetNonBendingCoor(); 
1800       
1801       fThetaX = muonTrack->GetThetaX();
1802       fThetaY = muonTrack->GetThetaY();
1803       
1804       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1805       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1806       fPxRec = fPzRec * TMath::Tan(fThetaX);
1807       fPyRec = fPzRec * TMath::Tan(fThetaY);
1808       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1809       
1810       //ChiSquare of the track if needed
1811       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1812       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1813       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1814       
1815       // total number of muons inside a vertex cut 
1816       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1817         nMuons++;
1818         if(fEPvector.Pt() > fLowPtCut) {
1819           nMu1GeV++; 
1820           if(fEPvector.Pt() > fHighPtCut) {
1821             nMu3GeV++; 
1822             if (fEPvector.Pt() > fVeryHighPtCut) {
1823               nMu10GeV++;
1824             }
1825           }
1826         }
1827       }
1828     }//muon track loop
1829     
1830     // Fill the event tags 
1831     if(ntrack != 0)
1832       meanPt = meanPt/ntrack;
1833     
1834     evTag->SetEventId(iEventNumber+1);
1835     if (vertexIn) {
1836       evTag->SetVertexX(vertexIn->GetXv());
1837       evTag->SetVertexY(vertexIn->GetYv());
1838       evTag->SetVertexZ(vertexIn->GetZv());
1839       evTag->SetVertexZError(vertexIn->GetZRes());
1840     }  
1841     evTag->SetVertexFlag(fVertexflag);
1842
1843     evTag->SetT0VertexZ(esd->GetT0zVertex());
1844     
1845     evTag->SetTriggerMask(esd->GetTriggerMask());
1846     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1847     
1848     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1849     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1850     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1851     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1852     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1853     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1854     
1855     
1856     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1857     evTag->SetNumOfPosTracks(nPos);
1858     evTag->SetNumOfNegTracks(nNeg);
1859     evTag->SetNumOfNeutrTracks(nNeutr);
1860     
1861     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1862     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1863     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1864     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1865     
1866     evTag->SetNumOfProtons(nProtons);
1867     evTag->SetNumOfKaons(nKaons);
1868     evTag->SetNumOfPions(nPions);
1869     evTag->SetNumOfMuons(nMuons);
1870     evTag->SetNumOfElectrons(nElectrons);
1871     evTag->SetNumOfPhotons(nGamas);
1872     evTag->SetNumOfPi0s(nPi0s);
1873     evTag->SetNumOfNeutrons(nNeutrons);
1874     evTag->SetNumOfKaon0s(nK0s);
1875     
1876     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1877     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1878     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1879     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1880     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1881     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1882     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1883     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1884     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1885     
1886     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1887     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1888     
1889     evTag->SetTotalMomentum(totalP);
1890     evTag->SetMeanPt(meanPt);
1891     evTag->SetMaxPt(maxPt);
1892     
1893     tag->SetRunId(iInitRunNumber);
1894     tag->AddEventTag(*evTag);
1895   }
1896   if(fLastEvent == -1) lastEvent = b->GetEntries();
1897   else lastEvent = fLastEvent;
1898         
1899   ttag.Fill();
1900   tag->Clear();
1901
1902   char fileName[256];
1903   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
1904           tag->GetRunId(),fFirstEvent,lastEvent );
1905   AliInfo(Form("writing tags to file %s", fileName));
1906   AliDebug(1, Form("writing tags to file %s", fileName));
1907  
1908   TFile* ftag = TFile::Open(fileName, "recreate");
1909   ftag->cd();
1910   ttag.Write();
1911   ftag->Close();
1912   file->cd();
1913   delete tag;
1914   delete evTag;
1915 }
1916
1917 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1918 {
1919   // Write space-points which are then used in the alignment procedures
1920   // For the moment only ITS, TRD and TPC
1921
1922   // Load TOF clusters
1923   if (fTracker[3]){
1924     fLoader[3]->LoadRecPoints("read");
1925     TTree* tree = fLoader[3]->TreeR();
1926     if (!tree) {
1927       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1928       return;
1929     }
1930     fTracker[3]->LoadClusters(tree);
1931   }
1932   Int_t ntracks = esd->GetNumberOfTracks();
1933   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1934     {
1935       AliESDtrack *track = esd->GetTrack(itrack);
1936       Int_t nsp = 0;
1937       Int_t idx[200];
1938       for (Int_t iDet = 3; iDet >= 0; iDet--)
1939         nsp += track->GetNcls(iDet);
1940       if (nsp) {
1941         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1942         track->SetTrackPointArray(sp);
1943         Int_t isptrack = 0;
1944         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1945           AliTracker *tracker = fTracker[iDet];
1946           if (!tracker) continue;
1947           Int_t nspdet = track->GetNcls(iDet);
1948           if (nspdet <= 0) continue;
1949           track->GetClusters(iDet,idx);
1950           AliTrackPoint p;
1951           Int_t isp = 0;
1952           Int_t isp2 = 0;
1953           while (isp < nspdet) {
1954             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1955             const Int_t kNTPCmax = 159;
1956             if (iDet==1 && isp2>kNTPCmax) break;   // to be fixed
1957             if (!isvalid) continue;
1958             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1959           }
1960         }       
1961       }
1962     }
1963   if (fTracker[3]){
1964     fTracker[3]->UnloadClusters();
1965     fLoader[3]->UnloadRecPoints();
1966   }
1967 }