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