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