6333f977d1647ed5f4d51abb4198b45397c7c35c
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // class for running the reconstruction                                      //
21 //                                                                           //
22 // Clusters and tracks are created for all detectors and all events by       //
23 // typing:                                                                   //
24 //                                                                           //
25 //   AliReconstruction rec;                                                  //
26 //   rec.Run();                                                              //
27 //                                                                           //
28 // The Run method returns kTRUE in case of successful execution.             //
29 //                                                                           //
30 // If the input to the reconstruction are not simulated digits but raw data, //
31 // this can be specified by an argument of the Run method or by the method   //
32 //                                                                           //
33 //   rec.SetInput("...");                                                    //
34 //                                                                           //
35 // The input formats and the corresponding argument are:                     //
36 // - DDL raw data files: directory name, ends with "/"                       //
37 // - raw data root file: root file name, extension ".root"                   //
38 // - raw data DATE file: DATE file name, any other non-empty string          //
39 // - MC root files     : empty string, default                               //
40 //                                                                           //
41 // By default all events are reconstructed. The reconstruction can be        //
42 // limited to a range of events by giving the index of the first and the     //
43 // last event as an argument to the Run method or by calling                 //
44 //                                                                           //
45 //   rec.SetEventRange(..., ...);                                            //
46 //                                                                           //
47 // The index -1 (default) can be used for the last event to indicate no      //
48 // upper limit of the event range.                                           //
49 //                                                                           //
50 // The name of the galice file can be changed from the default               //
51 // "galice.root" by passing it as argument to the AliReconstruction          //
52 // constructor or by                                                         //
53 //                                                                           //
54 //   rec.SetGAliceFile("...");                                               //
55 //                                                                           //
56 // The local reconstruction can be switched on or off for individual         //
57 // detectors by                                                              //
58 //                                                                           //
59 //   rec.SetRunLocalReconstruction("...");                                   //
60 //                                                                           //
61 // The argument is a (case sensitive) string with the names of the           //
62 // detectors separated by a space. The special string "ALL" selects all      //
63 // available detectors. This is the default.                                 //
64 //                                                                           //
65 // The reconstruction of the primary vertex position can be switched off by  //
66 //                                                                           //
67 //   rec.SetRunVertexFinder(kFALSE);                                         //
68 //                                                                           //
69 // The tracking and the creation of ESD tracks can be switched on for        //
70 // selected detectors by                                                     //
71 //                                                                           //
72 //   rec.SetRunTracking("...");                                              //
73 //                                                                           //
74 // Uniform/nonuniform field tracking switches (default: uniform field)       //
75 //                                                                           //
76 //   rec.SetUniformFieldTracking();  ( rec.SetNonuniformFieldTracking(); )   //
77 //                                                                           //
78 // The filling of additional ESD information can be steered by               //
79 //                                                                           //
80 //   rec.SetFillESD("...");                                                  //
81 //                                                                           //
82 // Again, for both methods the string specifies the list of detectors.       //
83 // The default is "ALL".                                                     //
84 //                                                                           //
85 // The call of the shortcut method                                           //
86 //                                                                           //
87 //   rec.SetRunReconstruction("...");                                        //
88 //                                                                           //
89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
90 // SetFillESD with the same detector selecting string as argument.           //
91 //                                                                           //
92 // The reconstruction requires digits or raw data as input. For the creation //
93 // of digits and raw data have a look at the class AliSimulation.            //
94 //                                                                           //
95 // For debug purposes the method SetCheckPointLevel can be used. If the      //
96 // argument is greater than 0, files with ESD events will be written after   //
97 // selected steps of the reconstruction for each event:                      //
98 //   level 1: after tracking and after filling of ESD (final)                //
99 //   level 2: in addition after each tracking step                           //
100 //   level 3: in addition after the filling of ESD for each detector         //
101 // If a final check point file exists for an event, this event will be       //
102 // skipped in the reconstruction. The tracking and the filling of ESD for    //
103 // a detector will be skipped as well, if the corresponding check point      //
104 // file exists. The ESD event will then be loaded from the file instead.     //
105 //                                                                           //
106 ///////////////////////////////////////////////////////////////////////////////
107
108 #include <TArrayF.h>
109 #include <TFile.h>
110 #include <TSystem.h>
111 #include <TROOT.h>
112 #include <TPluginManager.h>
113 #include <TStopwatch.h>
114
115 #include "AliReconstruction.h"
116 #include "AliReconstructor.h"
117 #include "AliLog.h"
118 #include "AliRunLoader.h"
119 #include "AliRun.h"
120 #include "AliRawReaderFile.h"
121 #include "AliRawReaderDate.h"
122 #include "AliRawReaderRoot.h"
123 #include "AliESD.h"
124 #include "AliESDVertex.h"
125 #include "AliTracker.h"
126 #include "AliVertexer.h"
127 #include "AliHeader.h"
128 #include "AliGenEventHeader.h"
129 #include "AliPID.h"
130 #include "AliESDpid.h"
131 //#include "AliMagF.h"
132
133
134
135 #include "AliRunTag.h"
136 //#include "AliLHCTag.h"
137 #include "AliDetectorTag.h"
138 #include "AliEventTag.h"
139
140 #include "AliTrackPointArray.h"
141
142 ClassImp(AliReconstruction)
143
144
145 //_____________________________________________________________________________
146 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
147
148 //_____________________________________________________________________________
149 AliReconstruction::AliReconstruction(const char* gAliceFilename,
150                                      const char* name, const char* title) :
151   TNamed(name, title),
152
153   fRunLocalReconstruction("ALL"),
154   fUniformField(kTRUE),
155   fRunVertexFinder(kTRUE),
156   fRunHLTTracking(kFALSE),
157   fRunTracking("ALL"),
158   fFillESD("ALL"),
159   fGAliceFileName(gAliceFilename),
160   fInput(""),
161   fFirstEvent(0),
162   fLastEvent(-1),
163   fStopOnError(kFALSE),
164   fCheckPointLevel(0),
165   fOptions(),
166
167   fRunLoader(NULL),
168   fRawReader(NULL),
169
170   fVertexer(NULL),
171
172   fWriteAlignmentData(kFALSE)
173 {
174 // create reconstruction object with default parameters
175   
176   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
177     fReconstructor[iDet] = NULL;
178     fLoader[iDet] = NULL;
179     fTracker[iDet] = NULL;
180   }
181   AliPID pid;
182 }
183
184 //_____________________________________________________________________________
185 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
186   TNamed(rec),
187
188   fRunLocalReconstruction(rec.fRunLocalReconstruction),
189   fUniformField(rec.fUniformField),
190   fRunVertexFinder(rec.fRunVertexFinder),
191   fRunHLTTracking(rec.fRunHLTTracking),
192   fRunTracking(rec.fRunTracking),
193   fFillESD(rec.fFillESD),
194   fGAliceFileName(rec.fGAliceFileName),
195   fInput(rec.fInput),
196   fFirstEvent(rec.fFirstEvent),
197   fLastEvent(rec.fLastEvent),
198   fStopOnError(rec.fStopOnError),
199   fCheckPointLevel(0),
200   fOptions(),
201
202   fRunLoader(NULL),
203   fRawReader(NULL),
204
205   fVertexer(NULL),
206
207   fWriteAlignmentData(rec.fWriteAlignmentData)
208 {
209 // copy constructor
210
211   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
212     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
213   }
214   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
215     fReconstructor[iDet] = NULL;
216     fLoader[iDet] = NULL;
217     fTracker[iDet] = NULL;
218   }
219 }
220
221 //_____________________________________________________________________________
222 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
223 {
224 // assignment operator
225
226   this->~AliReconstruction();
227   new(this) AliReconstruction(rec);
228   return *this;
229 }
230
231 //_____________________________________________________________________________
232 AliReconstruction::~AliReconstruction()
233 {
234 // clean up
235
236   CleanUp();
237   fOptions.Delete();
238 }
239
240
241 //_____________________________________________________________________________
242 void AliReconstruction::SetGAliceFile(const char* fileName)
243 {
244 // set the name of the galice file
245
246   fGAliceFileName = fileName;
247 }
248
249 //_____________________________________________________________________________
250 void AliReconstruction::SetOption(const char* detector, const char* option)
251 {
252 // set options for the reconstruction of a detector
253
254   TObject* obj = fOptions.FindObject(detector);
255   if (obj) fOptions.Remove(obj);
256   fOptions.Add(new TNamed(detector, option));
257 }
258
259
260 //_____________________________________________________________________________
261 Bool_t AliReconstruction::Run(const char* input,
262                               Int_t firstEvent, Int_t lastEvent)
263 {
264 // run the reconstruction
265
266   // set the input
267   if (!input) input = fInput.Data();
268   TString fileName(input);
269   if (fileName.EndsWith("/")) {
270     fRawReader = new AliRawReaderFile(fileName);
271   } else if (fileName.EndsWith(".root")) {
272     fRawReader = new AliRawReaderRoot(fileName);
273   } else if (!fileName.IsNull()) {
274     fRawReader = new AliRawReaderDate(fileName);
275     fRawReader->SelectEvents(7);
276   }
277
278   // get the run loader
279   if (!InitRunLoader()) return kFALSE;
280
281   // local reconstruction
282   if (!fRunLocalReconstruction.IsNull()) {
283     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
284       if (fStopOnError) {CleanUp(); return kFALSE;}
285     }
286   }
287 //  if (!fRunVertexFinder && fRunTracking.IsNull() && 
288 //      fFillESD.IsNull()) return kTRUE;
289
290   // get vertexer
291   if (fRunVertexFinder && !CreateVertexer()) {
292     if (fStopOnError) {
293       CleanUp(); 
294       return kFALSE;
295     }
296   }
297
298   // get trackers
299   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
300     if (fStopOnError) {
301       CleanUp(); 
302       return kFALSE;
303     }      
304   }
305
306
307   TStopwatch stopwatch;
308   stopwatch.Start();
309
310   // get the possibly already existing ESD file and tree
311   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
312   TFile* fileOld = NULL;
313   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
314   if (!gSystem->AccessPathName("AliESDs.root")){
315     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
316     fileOld = TFile::Open("AliESDs.old.root");
317     if (fileOld && fileOld->IsOpen()) {
318       treeOld = (TTree*) fileOld->Get("esdTree");
319       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
320       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
321       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
322     }
323   }
324
325   // create the ESD output file and tree
326   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
327   if (!file->IsOpen()) {
328     AliError("opening AliESDs.root failed");
329     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
330   }
331   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
332   tree->Branch("ESD", "AliESD", &esd);
333   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
334   hlttree->Branch("ESD", "AliESD", &hltesd);
335   delete esd; delete hltesd;
336   esd = NULL; hltesd = NULL;
337   gROOT->cd();
338
339   // loop over events
340   if (fRawReader) fRawReader->RewindEvents();
341   
342   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
343     if (fRawReader) fRawReader->NextEvent();
344     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
345       // copy old ESD to the new one
346       if (treeOld) {
347         treeOld->SetBranchAddress("ESD", &esd);
348         treeOld->GetEntry(iEvent);
349       }
350       tree->Fill();
351       if (hlttreeOld) {
352         hlttreeOld->SetBranchAddress("ESD", &hltesd);
353         hlttreeOld->GetEntry(iEvent);
354       }
355       hlttree->Fill();
356       continue;
357     }
358
359     AliInfo(Form("processing event %d", iEvent));
360     fRunLoader->GetEvent(iEvent);
361
362     char fileName[256];
363     sprintf(fileName, "ESD_%d.%d_final.root", 
364             fRunLoader->GetHeader()->GetRun(), 
365             fRunLoader->GetHeader()->GetEventNrInRun());
366     if (!gSystem->AccessPathName(fileName)) continue;
367
368     // local reconstruction
369     if (!fRunLocalReconstruction.IsNull()) {
370       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
371         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
372       }
373     }
374
375     esd = new AliESD; hltesd = new AliESD;
376     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
377     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
378     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
379     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
380     if (gAlice) {
381       esd->SetMagneticField(gAlice->Field()->SolenoidField());
382       hltesd->SetMagneticField(gAlice->Field()->SolenoidField());
383     } else {
384       // ???
385     }
386
387     // vertex finder
388     if (fRunVertexFinder) {
389       if (!ReadESD(esd, "vertex")) {
390         if (!RunVertexFinder(esd)) {
391           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
392         }
393         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
394       }
395     }
396
397     // HLT tracking
398     if (!fRunTracking.IsNull()) {
399       if (fRunHLTTracking) {
400         hltesd->SetVertex(esd->GetVertex());
401         if (!RunHLTTracking(hltesd)) {
402           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
403         }
404       }
405     }
406
407     // barrel tracking
408     if (!fRunTracking.IsNull()) {
409       if (!ReadESD(esd, "tracking")) {
410         if (!RunTracking(esd)) {
411           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
412         }
413         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
414       }
415     }
416
417     // fill ESD
418     if (!fFillESD.IsNull()) {
419       if (!FillESD(esd, fFillESD)) {
420         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
421       }
422     }
423
424     // combined PID
425     AliESDpid::MakePID(esd);
426     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
427
428     // write ESD
429     tree->Fill();
430     // write HLT ESD
431     hlttree->Fill();
432
433     if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
434  
435     delete esd; delete hltesd;
436     esd = NULL; hltesd = NULL;
437   }
438
439   AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
440                stopwatch.RealTime(),stopwatch.CpuTime()));
441
442   file->cd();
443   tree->Write();
444   hlttree->Write();
445
446   // Create tags for the events in the ESD tree (the ESD tree is always present)
447   // In case of empty events the tags will contain dummy values
448   CreateTag(file);
449   CleanUp(file, fileOld);
450
451   return kTRUE;
452 }
453
454
455 //_____________________________________________________________________________
456 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
457 {
458 // run the local reconstruction
459
460   TStopwatch stopwatch;
461   stopwatch.Start();
462
463   TString detStr = detectors;
464   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
465     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
466     AliReconstructor* reconstructor = GetReconstructor(iDet);
467     if (!reconstructor) continue;
468     if (reconstructor->HasLocalReconstruction()) continue;
469
470     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
471     TStopwatch stopwatchDet;
472     stopwatchDet.Start();
473     if (fRawReader) {
474       fRawReader->RewindEvents();
475       reconstructor->Reconstruct(fRunLoader, fRawReader);
476     } else {
477       reconstructor->Reconstruct(fRunLoader);
478     }
479     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
480                  fgkDetectorName[iDet],
481                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
482   }
483
484   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
485     AliError(Form("the following detectors were not found: %s",
486                   detStr.Data()));
487     if (fStopOnError) return kFALSE;
488   }
489
490   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
491                stopwatch.RealTime(),stopwatch.CpuTime()));
492
493   return kTRUE;
494 }
495
496 //_____________________________________________________________________________
497 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
498 {
499 // run the local reconstruction
500
501   TStopwatch stopwatch;
502   stopwatch.Start();
503
504   TString detStr = detectors;
505   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
506     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
507     AliReconstructor* reconstructor = GetReconstructor(iDet);
508     if (!reconstructor) continue;
509     AliLoader* loader = fLoader[iDet];
510
511     // conversion of digits
512     if (fRawReader && reconstructor->HasDigitConversion()) {
513       AliInfo(Form("converting raw data digits into root objects for %s", 
514                    fgkDetectorName[iDet]));
515       TStopwatch stopwatchDet;
516       stopwatchDet.Start();
517       loader->LoadDigits("update");
518       loader->CleanDigits();
519       loader->MakeDigitsContainer();
520       TTree* digitsTree = loader->TreeD();
521       reconstructor->ConvertDigits(fRawReader, digitsTree);
522       loader->WriteDigits("OVERWRITE");
523       loader->UnloadDigits();
524       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
525                    fgkDetectorName[iDet],
526                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
527     }
528
529     // local reconstruction
530     if (!reconstructor->HasLocalReconstruction()) continue;
531     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
532     TStopwatch stopwatchDet;
533     stopwatchDet.Start();
534     loader->LoadRecPoints("update");
535     loader->CleanRecPoints();
536     loader->MakeRecPointsContainer();
537     TTree* clustersTree = loader->TreeR();
538     if (fRawReader && !reconstructor->HasDigitConversion()) {
539       reconstructor->Reconstruct(fRawReader, clustersTree);
540     } else {
541       loader->LoadDigits("read");
542       TTree* digitsTree = loader->TreeD();
543       if (!digitsTree) {
544         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
545         if (fStopOnError) return kFALSE;
546       } else {
547         reconstructor->Reconstruct(digitsTree, clustersTree);
548       }
549       loader->UnloadDigits();
550     }
551     loader->WriteRecPoints("OVERWRITE");
552     loader->UnloadRecPoints();
553     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
554                     fgkDetectorName[iDet],
555                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
556   }
557
558   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
559     AliError(Form("the following detectors were not found: %s",
560                   detStr.Data()));
561     if (fStopOnError) return kFALSE;
562   }
563   
564   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
565                stopwatch.RealTime(),stopwatch.CpuTime()));
566
567   return kTRUE;
568 }
569
570 //_____________________________________________________________________________
571 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
572 {
573 // run the barrel tracking
574
575   TStopwatch stopwatch;
576   stopwatch.Start();
577
578   AliESDVertex* vertex = NULL;
579   Double_t vtxPos[3] = {0, 0, 0};
580   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
581   TArrayF mcVertex(3); 
582   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
583     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
584     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
585   }
586
587   if (fVertexer) {
588     AliInfo("running the ITS vertex finder");
589     if (fLoader[0]) fLoader[0]->LoadRecPoints();
590     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
591     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
592     if(!vertex){
593       AliWarning("Vertex not found");
594       vertex = new AliESDVertex();
595     }
596     else {
597       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
598     }
599
600   } else {
601     AliInfo("getting the primary vertex from MC");
602     vertex = new AliESDVertex(vtxPos, vtxErr);
603   }
604
605   if (vertex) {
606     vertex->GetXYZ(vtxPos);
607     vertex->GetSigmaXYZ(vtxErr);
608   } else {
609     AliWarning("no vertex reconstructed");
610     vertex = new AliESDVertex(vtxPos, vtxErr);
611   }
612   esd->SetVertex(vertex);
613   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
614     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
615   }  
616   delete vertex;
617
618   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
619                stopwatch.RealTime(),stopwatch.CpuTime()));
620
621   return kTRUE;
622 }
623
624 //_____________________________________________________________________________
625 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
626 {
627 // run the HLT barrel tracking
628
629   TStopwatch stopwatch;
630   stopwatch.Start();
631
632   if (!fRunLoader) {
633     AliError("Missing runLoader!");
634     return kFALSE;
635   }
636
637   AliInfo("running HLT tracking");
638
639   // Get a pointer to the HLT reconstructor
640   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
641   if (!reconstructor) return kFALSE;
642
643   // TPC + ITS
644   for (Int_t iDet = 1; iDet >= 0; iDet--) {
645     TString detName = fgkDetectorName[iDet];
646     AliDebug(1, Form("%s HLT tracking", detName.Data()));
647     reconstructor->SetOption(detName.Data());
648     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
649     if (!tracker) {
650       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
651       if (fStopOnError) return kFALSE;
652       continue;
653     }
654     Double_t vtxPos[3];
655     Double_t vtxErr[3]={0.005,0.005,0.010};
656     const AliESDVertex *vertex = esd->GetVertex();
657     vertex->GetXYZ(vtxPos);
658     tracker->SetVertex(vtxPos,vtxErr);
659     if(iDet != 1) {
660       fLoader[iDet]->LoadRecPoints("read");
661       TTree* tree = fLoader[iDet]->TreeR();
662       if (!tree) {
663         AliError(Form("Can't get the %s cluster tree", detName.Data()));
664         return kFALSE;
665       }
666       tracker->LoadClusters(tree);
667     }
668     if (tracker->Clusters2Tracks(esd) != 0) {
669       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
670       return kFALSE;
671     }
672     if(iDet != 1) {
673       tracker->UnloadClusters();
674     }
675     delete tracker;
676   }
677
678   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
679                stopwatch.RealTime(),stopwatch.CpuTime()));
680
681   return kTRUE;
682 }
683
684 //_____________________________________________________________________________
685 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
686 {
687 // run the barrel tracking
688
689   TStopwatch stopwatch;
690   stopwatch.Start();
691
692   AliInfo("running tracking");
693
694   // pass 1: TPC + ITS inwards
695   for (Int_t iDet = 1; iDet >= 0; iDet--) {
696     if (!fTracker[iDet]) continue;
697     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
698
699     // load clusters
700     fLoader[iDet]->LoadRecPoints("read");
701     TTree* tree = fLoader[iDet]->TreeR();
702     if (!tree) {
703       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
704       return kFALSE;
705     }
706     fTracker[iDet]->LoadClusters(tree);
707
708     // run tracking
709     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
710       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
711       return kFALSE;
712     }
713     if (fCheckPointLevel > 1) {
714       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
715     }
716     // preliminary PID in TPC needed by the ITS tracker
717     if (iDet == 1) {
718       GetReconstructor(1)->FillESD(fRunLoader, esd);
719       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
720       AliESDpid::MakePID(esd);
721     }
722   }
723
724   // pass 2: ALL backwards
725   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
726     if (!fTracker[iDet]) continue;
727     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
728
729     // load clusters
730     if (iDet > 1) {     // all except ITS, TPC
731       TTree* tree = NULL;
732       fLoader[iDet]->LoadRecPoints("read");
733       tree = fLoader[iDet]->TreeR();
734       if (!tree) {
735         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
736         return kFALSE;
737       }
738       fTracker[iDet]->LoadClusters(tree);
739     }
740
741     // run tracking
742     if (fTracker[iDet]->PropagateBack(esd) != 0) {
743       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
744       return kFALSE;
745     }
746     if (fCheckPointLevel > 1) {
747       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
748     }
749
750     // unload clusters
751     if (iDet > 2) {     // all except ITS, TPC, TRD
752       fTracker[iDet]->UnloadClusters();
753       fLoader[iDet]->UnloadRecPoints();
754     }
755   }
756
757   // write space-points to the ESD in case alignment data output
758   // is switched on
759   if (fWriteAlignmentData)
760     WriteAlignmentData(esd);
761
762   // pass 3: TRD + TPC + ITS refit inwards
763   for (Int_t iDet = 2; iDet >= 0; iDet--) {
764     if (!fTracker[iDet]) continue;
765     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
766
767     // run tracking
768     if (fTracker[iDet]->RefitInward(esd) != 0) {
769       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
770       return kFALSE;
771     }
772     if (fCheckPointLevel > 1) {
773       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
774     }
775
776     // unload clusters
777     fTracker[iDet]->UnloadClusters();
778     fLoader[iDet]->UnloadRecPoints();
779   }
780
781   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
782                stopwatch.RealTime(),stopwatch.CpuTime()));
783
784   return kTRUE;
785 }
786
787 //_____________________________________________________________________________
788 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
789 {
790 // fill the event summary data
791
792   TStopwatch stopwatch;
793   stopwatch.Start();
794   AliInfo("filling ESD");
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
802     if (!ReadESD(esd, fgkDetectorName[iDet])) {
803       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
804       TTree* clustersTree = NULL;
805       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
806         fLoader[iDet]->LoadRecPoints("read");
807         clustersTree = fLoader[iDet]->TreeR();
808         if (!clustersTree) {
809           AliError(Form("Can't get the %s clusters tree", 
810                         fgkDetectorName[iDet]));
811           if (fStopOnError) return kFALSE;
812         }
813       }
814       if (fRawReader && !reconstructor->HasDigitConversion()) {
815         reconstructor->FillESD(fRawReader, clustersTree, esd);
816       } else {
817         TTree* digitsTree = NULL;
818         if (fLoader[iDet]) {
819           fLoader[iDet]->LoadDigits("read");
820           digitsTree = fLoader[iDet]->TreeD();
821           if (!digitsTree) {
822             AliError(Form("Can't get the %s digits tree", 
823                           fgkDetectorName[iDet]));
824             if (fStopOnError) return kFALSE;
825           }
826         }
827         reconstructor->FillESD(digitsTree, clustersTree, esd);
828         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
829       }
830       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
831         fLoader[iDet]->UnloadRecPoints();
832       }
833
834       if (fRawReader) {
835         reconstructor->FillESD(fRunLoader, fRawReader, esd);
836       } else {
837         reconstructor->FillESD(fRunLoader, esd);
838       }
839       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
840     }
841   }
842
843   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
844     AliError(Form("the following detectors were not found: %s", 
845                   detStr.Data()));
846     if (fStopOnError) return kFALSE;
847   }
848
849   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
850                stopwatch.RealTime(),stopwatch.CpuTime()));
851
852   return kTRUE;
853 }
854
855
856 //_____________________________________________________________________________
857 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
858 {
859 // check whether detName is contained in detectors
860 // if yes, it is removed from detectors
861
862   // check if all detectors are selected
863   if ((detectors.CompareTo("ALL") == 0) ||
864       detectors.BeginsWith("ALL ") ||
865       detectors.EndsWith(" ALL") ||
866       detectors.Contains(" ALL ")) {
867     detectors = "ALL";
868     return kTRUE;
869   }
870
871   // search for the given detector
872   Bool_t result = kFALSE;
873   if ((detectors.CompareTo(detName) == 0) ||
874       detectors.BeginsWith(detName+" ") ||
875       detectors.EndsWith(" "+detName) ||
876       detectors.Contains(" "+detName+" ")) {
877     detectors.ReplaceAll(detName, "");
878     result = kTRUE;
879   }
880
881   // clean up the detectors string
882   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
883   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
884   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
885
886   return result;
887 }
888
889 //_____________________________________________________________________________
890 Bool_t AliReconstruction::InitRunLoader()
891 {
892 // get or create the run loader
893
894   if (gAlice) delete gAlice;
895   gAlice = NULL;
896
897   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
898     // load all base libraries to get the loader classes
899     TString libs = gSystem->GetLibraries();
900     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
901       TString detName = fgkDetectorName[iDet];
902       if (detName == "HLT") continue;
903       if (libs.Contains("lib" + detName + "base.so")) continue;
904       gSystem->Load("lib" + detName + "base.so");
905     }
906     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
907     if (!fRunLoader) {
908       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
909       CleanUp();
910       return kFALSE;
911     }
912     fRunLoader->CdGAFile();
913     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
914       if (fRunLoader->LoadgAlice() == 0) {
915         gAlice = fRunLoader->GetAliRun();
916         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
917         AliExternalTrackParam::SetFieldMap(gAlice->Field());
918         if(fUniformField)
919           AliExternalTrackParam::SetUniformFieldTracking();
920         else
921           AliExternalTrackParam::SetNonuniformFieldTracking();
922       }
923     }
924     if (!gAlice && !fRawReader) {
925       AliError(Form("no gAlice object found in file %s",
926                     fGAliceFileName.Data()));
927       CleanUp();
928       return kFALSE;
929     }
930
931   } else {               // galice.root does not exist
932     if (!fRawReader) {
933       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
934       CleanUp();
935       return kFALSE;
936     }
937     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
938                                     AliConfig::GetDefaultEventFolderName(),
939                                     "recreate");
940     if (!fRunLoader) {
941       AliError(Form("could not create run loader in file %s", 
942                     fGAliceFileName.Data()));
943       CleanUp();
944       return kFALSE;
945     }
946     fRunLoader->MakeTree("E");
947     Int_t iEvent = 0;
948     while (fRawReader->NextEvent()) {
949       fRunLoader->SetEventNumber(iEvent);
950       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
951                                      iEvent, iEvent);
952       fRunLoader->MakeTree("H");
953       fRunLoader->TreeE()->Fill();
954       iEvent++;
955     }
956     fRawReader->RewindEvents();
957     fRunLoader->WriteHeader("OVERWRITE");
958     fRunLoader->CdGAFile();
959     fRunLoader->Write(0, TObject::kOverwrite);
960 //    AliTracker::SetFieldMap(???);
961   }
962
963   return kTRUE;
964 }
965
966 //_____________________________________________________________________________
967 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
968 {
969 // get the reconstructor object and the loader for a detector
970
971   if (fReconstructor[iDet]) return fReconstructor[iDet];
972
973   // load the reconstructor object
974   TPluginManager* pluginManager = gROOT->GetPluginManager();
975   TString detName = fgkDetectorName[iDet];
976   TString recName = "Ali" + detName + "Reconstructor";
977   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
978
979   if (detName == "HLT") {
980     if (!gROOT->GetClass("AliLevel3")) {
981       gSystem->Load("libAliL3Src.so");
982       gSystem->Load("libAliL3Misc.so");
983       gSystem->Load("libAliL3Hough.so");
984       gSystem->Load("libAliL3Comp.so");
985     }
986   }
987
988   AliReconstructor* reconstructor = NULL;
989   // first check if a plugin is defined for the reconstructor
990   TPluginHandler* pluginHandler = 
991     pluginManager->FindHandler("AliReconstructor", detName);
992   // if not, add a plugin for it
993   if (!pluginHandler) {
994     AliDebug(1, Form("defining plugin for %s", recName.Data()));
995     TString libs = gSystem->GetLibraries();
996     if (libs.Contains("lib" + detName + "base.so") ||
997         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
998       pluginManager->AddHandler("AliReconstructor", detName, 
999                                 recName, detName + "rec", recName + "()");
1000     } else {
1001       pluginManager->AddHandler("AliReconstructor", detName, 
1002                                 recName, detName, recName + "()");
1003     }
1004     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1005   }
1006   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1007     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1008   }
1009   if (reconstructor) {
1010     TObject* obj = fOptions.FindObject(detName.Data());
1011     if (obj) reconstructor->SetOption(obj->GetTitle());
1012     reconstructor->Init(fRunLoader);
1013     fReconstructor[iDet] = reconstructor;
1014   }
1015
1016   // get or create the loader
1017   if (detName != "HLT") {
1018     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1019     if (!fLoader[iDet]) {
1020       AliConfig::Instance()
1021         ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
1022                                 detName, detName);
1023       // first check if a plugin is defined for the loader
1024       TPluginHandler* pluginHandler = 
1025         pluginManager->FindHandler("AliLoader", detName);
1026       // if not, add a plugin for it
1027       if (!pluginHandler) {
1028         TString loaderName = "Ali" + detName + "Loader";
1029         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1030         pluginManager->AddHandler("AliLoader", detName, 
1031                                   loaderName, detName + "base", 
1032                                   loaderName + "(const char*, TFolder*)");
1033         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1034       }
1035       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1036         fLoader[iDet] = 
1037           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
1038                                                  fRunLoader->GetEventFolder());
1039       }
1040       if (!fLoader[iDet]) {   // use default loader
1041         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1042       }
1043       if (!fLoader[iDet]) {
1044         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1045         if (fStopOnError) return NULL;
1046       } else {
1047         fRunLoader->AddLoader(fLoader[iDet]);
1048         fRunLoader->CdGAFile();
1049         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1050         fRunLoader->Write(0, TObject::kOverwrite);
1051       }
1052     }
1053   }
1054       
1055   return reconstructor;
1056 }
1057
1058 //_____________________________________________________________________________
1059 Bool_t AliReconstruction::CreateVertexer()
1060 {
1061 // create the vertexer
1062
1063   fVertexer = NULL;
1064   AliReconstructor* itsReconstructor = GetReconstructor(0);
1065   if (itsReconstructor) {
1066     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1067   }
1068   if (!fVertexer) {
1069     AliWarning("couldn't create a vertexer for ITS");
1070     if (fStopOnError) return kFALSE;
1071   }
1072
1073   return kTRUE;
1074 }
1075
1076 //_____________________________________________________________________________
1077 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1078 {
1079 // create the trackers
1080
1081   TString detStr = detectors;
1082   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1083     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1084     AliReconstructor* reconstructor = GetReconstructor(iDet);
1085     if (!reconstructor) continue;
1086     TString detName = fgkDetectorName[iDet];
1087     if (detName == "HLT") {
1088       fRunHLTTracking = kTRUE;
1089       continue;
1090     }
1091
1092     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1093     if (!fTracker[iDet] && (iDet < 7)) {
1094       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1095       if (fStopOnError) return kFALSE;
1096     }
1097   }
1098
1099   return kTRUE;
1100 }
1101
1102 //_____________________________________________________________________________
1103 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1104 {
1105 // delete trackers and the run loader and close and delete the file
1106
1107   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1108     delete fReconstructor[iDet];
1109     fReconstructor[iDet] = NULL;
1110     fLoader[iDet] = NULL;
1111     delete fTracker[iDet];
1112     fTracker[iDet] = NULL;
1113   }
1114   delete fVertexer;
1115   fVertexer = NULL;
1116
1117   delete fRunLoader;
1118   fRunLoader = NULL;
1119   delete fRawReader;
1120   fRawReader = NULL;
1121
1122   if (file) {
1123     file->Close();
1124     delete file;
1125   }
1126
1127   if (fileOld) {
1128     fileOld->Close();
1129     delete fileOld;
1130     gSystem->Unlink("AliESDs.old.root");
1131   }
1132 }
1133
1134
1135 //_____________________________________________________________________________
1136 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1137 {
1138 // read the ESD event from a file
1139
1140   if (!esd) return kFALSE;
1141   char fileName[256];
1142   sprintf(fileName, "ESD_%d.%d_%s.root", 
1143           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1144   if (gSystem->AccessPathName(fileName)) return kFALSE;
1145
1146   AliInfo(Form("reading ESD from file %s", fileName));
1147   AliDebug(1, Form("reading ESD from file %s", fileName));
1148   TFile* file = TFile::Open(fileName);
1149   if (!file || !file->IsOpen()) {
1150     AliError(Form("opening %s failed", fileName));
1151     delete file;
1152     return kFALSE;
1153   }
1154
1155   gROOT->cd();
1156   delete esd;
1157   esd = (AliESD*) file->Get("ESD");
1158   file->Close();
1159   delete file;
1160   return kTRUE;
1161 }
1162
1163 //_____________________________________________________________________________
1164 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1165 {
1166 // write the ESD event to a file
1167
1168   if (!esd) return;
1169   char fileName[256];
1170   sprintf(fileName, "ESD_%d.%d_%s.root", 
1171           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1172
1173   AliDebug(1, Form("writing ESD to file %s", fileName));
1174   TFile* file = TFile::Open(fileName, "recreate");
1175   if (!file || !file->IsOpen()) {
1176     AliError(Form("opening %s failed", fileName));
1177   } else {
1178     esd->Write("ESD");
1179     file->Close();
1180   }
1181   delete file;
1182 }
1183
1184
1185
1186
1187 //_____________________________________________________________________________
1188 void AliReconstruction::CreateTag(TFile* file)
1189 {
1190   // Creates the tags for all the events in a given ESD file
1191   Int_t ntrack;
1192   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1193   Int_t nPos, nNeg, nNeutr;
1194   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1195   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1196   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1197   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1198   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1199
1200   AliRunTag *tag = new AliRunTag();
1201   AliDetectorTag *detTag = new AliDetectorTag();
1202   AliEventTag *evTag = new AliEventTag();
1203   TTree ttag("T","A Tree with event tags");
1204   TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
1205   btag->SetCompressionLevel(9);
1206
1207   AliInfo(Form("Creating the tags......."));    
1208   
1209   if (!file || !file->IsOpen()) {
1210     AliError(Form("opening failed"));
1211     delete file;
1212     return ;
1213   }
1214
1215   TTree *t = (TTree*) file->Get("esdTree");
1216   TBranch * b = t->GetBranch("ESD");
1217   AliESD *esd = 0;
1218   b->SetAddress(&esd);
1219
1220   tag->SetRunId(esd->GetRunNumber());
1221
1222   Int_t firstEvent = 0,lastEvent = 0;
1223   Int_t iNumberOfEvents = b->GetEntries();
1224   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++)
1225     {
1226       ntrack = 0;
1227       nPos = 0;
1228       nNeg = 0;
1229       nNeutr =0;
1230       nK0s = 0;
1231       nNeutrons = 0;
1232       nPi0s = 0;
1233       nGamas = 0;
1234       nProtons = 0;
1235       nKaons = 0;
1236       nPions = 0;
1237       nMuons = 0;
1238       nElectrons = 0;     
1239       nCh1GeV = 0;
1240       nCh3GeV = 0;
1241       nCh10GeV = 0;
1242       nMu1GeV = 0;
1243       nMu3GeV = 0;
1244       nMu10GeV = 0;
1245       nEl1GeV = 0;
1246       nEl3GeV = 0;
1247       nEl10GeV = 0;
1248       maxPt = .0;
1249       meanPt = .0;
1250       totalP = .0;
1251
1252       b->GetEntry(iEventNumber);
1253       const AliESDVertex * vertexIn = esd->GetVertex();
1254
1255       for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++)
1256         {
1257           AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1258           UInt_t status = esdTrack->GetStatus();
1259           
1260           //select only tracks with ITS refit
1261           if ((status&AliESDtrack::kITSrefit)==0) continue;
1262           
1263           //select only tracks with TPC refit-->remove extremely high Pt tracks
1264           if ((status&AliESDtrack::kTPCrefit)==0) continue;
1265           
1266           //select only tracks with the "combined PID"
1267           if ((status&AliESDtrack::kESDpid)==0) continue;
1268                   Double_t p[3];
1269           esdTrack->GetPxPyPz(p);
1270           Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1271           Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1272           totalP += momentum;
1273           meanPt += fPt;
1274           if(fPt > maxPt)
1275             maxPt = fPt;
1276           
1277           if(esdTrack->GetSign() > 0)
1278             {
1279               nPos++;
1280               if(fPt > 1.0)
1281                 nCh1GeV++;
1282               if(fPt > 3.0)
1283                 nCh3GeV++;
1284               if(fPt > 10.0)
1285                 nCh10GeV++;
1286             }
1287           if(esdTrack->GetSign() < 0)
1288             {
1289               nNeg++;
1290               if(fPt > 1.0)
1291                 nCh1GeV++;
1292               if(fPt > 3.0)
1293                 nCh3GeV++;
1294               if(fPt > 10.0)
1295                 nCh10GeV++;
1296             }
1297           if(esdTrack->GetSign() == 0)
1298             nNeutr++;
1299           
1300           //PID
1301           Double_t prob[10];
1302           esdTrack->GetESDpid(prob);
1303                     
1304           //K0s
1305           if ((prob[8]>prob[7])&&(prob[8]>prob[6])&&(prob[8]>prob[5])&&(prob[8]>prob[4])&&(prob[8]>prob[3])&&(prob[8]>prob[2])&&(prob[8]>prob[1])&&(prob[8]>prob[0]))
1306             nK0s++;
1307           //neutrons
1308           if ((prob[7]>prob[8])&&(prob[7]>prob[6])&&(prob[7]>prob[5])&&(prob[7]>prob[4])&&(prob[7]>prob[3])&&(prob[7]>prob[2])&&(prob[7]>prob[1])&&(prob[7]>prob[0]))
1309             nNeutrons++; 
1310           //pi0s
1311           if ((prob[6]>prob[8])&&(prob[6]>prob[7])&&(prob[6]>prob[5])&&(prob[6]>prob[4])&&(prob[6]>prob[3])&&(prob[6]>prob[2])&&(prob[6]>prob[1])&&(prob[6]>prob[0]))
1312             nPi0s++;
1313           //gamas
1314           if ((prob[5]>prob[8])&&(prob[5]>prob[7])&&(prob[5]>prob[6])&&(prob[5]>prob[4])&&(prob[5]>prob[3])&&(prob[5]>prob[2])&&(prob[5]>prob[1])&&(prob[5]>prob[0]))
1315             nGamas++;
1316           //protons
1317           if ((prob[4]>prob[8])&&(prob[4]>prob[7])&&(prob[4]>prob[6])&&(prob[4]>prob[5])&&(prob[4]>prob[3])&&(prob[4]>prob[2])&&(prob[4]>prob[1])&&(prob[4]>prob[0]))
1318             nProtons++;
1319           //kaons
1320           if ((prob[3]>prob[8])&&(prob[3]>prob[7])&&(prob[3]>prob[6])&&(prob[3]>prob[5])&&(prob[3]>prob[4])&&(prob[3]>prob[2])&&(prob[3]>prob[1])&&(prob[3]>prob[0]))
1321             nKaons++;
1322           //kaons
1323           if ((prob[2]>prob[8])&&(prob[2]>prob[7])&&(prob[2]>prob[6])&&(prob[2]>prob[5])&&(prob[2]>prob[4])&&(prob[2]>prob[3])&&(prob[2]>prob[1])&&(prob[2]>prob[0]))
1324             nPions++; 
1325           //muons
1326           if ((prob[1]>prob[8])&&(prob[1]>prob[7])&&(prob[1]>prob[6])&&(prob[1]>prob[5])&&(prob[1]>prob[4])&&(prob[1]>prob[3])&&(prob[1]>prob[2])&&(prob[1]>prob[0]))
1327             {
1328               nMuons++;
1329               if(fPt > 1.0)
1330                 nMu1GeV++;
1331               if(fPt > 3.0)
1332                 nMu3GeV++;
1333               if(fPt > 10.0)
1334                 nMu10GeV++;
1335             }
1336           //electrons
1337           if ((prob[0]>prob[8])&&(prob[0]>prob[7])&&(prob[0]>prob[6])&&(prob[0]>prob[5])&&(prob[0]>prob[4])&&(prob[0]>prob[3])&&(prob[0]>prob[2])&&(prob[0]>prob[1]))
1338             {
1339               nElectrons++;
1340               if(fPt > 1.0)
1341                 nEl1GeV++;
1342               if(fPt > 3.0)
1343                 nEl3GeV++;
1344               if(fPt > 10.0)
1345                 nEl10GeV++;
1346             }
1347           
1348           
1349           
1350           ntrack++;
1351         }//track loop
1352       // Fill the event tags 
1353       if(ntrack != 0)
1354         meanPt = meanPt/ntrack;
1355       
1356       evTag->SetEventId(iEventNumber+1);
1357       evTag->SetVertexX(vertexIn->GetXv());
1358       evTag->SetVertexY(vertexIn->GetYv());
1359       evTag->SetVertexZ(vertexIn->GetZv());
1360       
1361       evTag->SetT0VertexZ(esd->GetT0zVertex());
1362       
1363       evTag->SetTrigger(esd->GetTrigger());
1364       
1365       evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
1366       evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
1367       evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1368       evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1369       
1370       
1371       evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1372       evTag->SetNumOfPosTracks(nPos);
1373       evTag->SetNumOfNegTracks(nNeg);
1374       evTag->SetNumOfNeutrTracks(nNeutr);
1375       
1376       evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1377       evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1378       evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1379       evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1380       
1381       evTag->SetNumOfProtons(nProtons);
1382       evTag->SetNumOfKaons(nKaons);
1383       evTag->SetNumOfPions(nPions);
1384       evTag->SetNumOfMuons(nMuons);
1385       evTag->SetNumOfElectrons(nElectrons);
1386       evTag->SetNumOfPhotons(nGamas);
1387       evTag->SetNumOfPi0s(nPi0s);
1388       evTag->SetNumOfNeutrons(nNeutrons);
1389       evTag->SetNumOfKaon0s(nK0s);
1390       
1391       evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1392       evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1393       evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1394       evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1395       evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1396       evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1397       evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1398       evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1399       evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1400       
1401       evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
1402       evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
1403       
1404       evTag->SetTotalMomentum(totalP);
1405       evTag->SetMeanPt(meanPt);
1406       evTag->SetMaxPt(maxPt);
1407   
1408       tag->AddEventTag(*evTag);
1409     }
1410   lastEvent = iNumberOfEvents;
1411         
1412   ttag.Fill();
1413   tag->Clear();
1414
1415   char fileName[256];
1416   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
1417           tag->GetRunId(),firstEvent,lastEvent );
1418   AliInfo(Form("writing tags to file %s", fileName));
1419   AliDebug(1, Form("writing tags to file %s", fileName));
1420  
1421   TFile* ftag = TFile::Open(fileName, "recreate");
1422   ftag->cd();
1423   ttag.Write();
1424   ftag->Close();
1425   file->cd();
1426   delete tag;
1427   delete detTag;
1428   delete evTag;
1429 }
1430
1431 void AliReconstruction::WriteAlignmentData(AliESD* esd)
1432 {
1433   // Write space-points which are then used in the alignment procedures
1434   // For the moment only ITS, TRD and TPC
1435
1436   // Load TOF clusters
1437   fLoader[3]->LoadRecPoints("read");
1438   TTree* tree = fLoader[3]->TreeR();
1439   if (!tree) {
1440     AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1441     return;
1442   }
1443   fTracker[3]->LoadClusters(tree);
1444
1445   Int_t ntracks = esd->GetNumberOfTracks();
1446   for (Int_t itrack = 0; itrack < ntracks; itrack++)
1447     {
1448       AliESDtrack *track = esd->GetTrack(itrack);
1449       Int_t nsp = 0;
1450       UInt_t idx[200];
1451       for (Int_t iDet = 3; iDet >= 0; iDet--)
1452         nsp += track->GetNcls(iDet);
1453       if (nsp) {
1454         AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1455         track->SetTrackPointArray(sp);
1456         Int_t isptrack = 0;
1457         for (Int_t iDet = 3; iDet >= 0; iDet--) {
1458           AliTracker *tracker = fTracker[iDet];
1459           if (!tracker) continue;
1460           Int_t nspdet = track->GetNcls(iDet);
1461           cout<<iDet<<" "<<nspdet<<endl;
1462           if (nspdet <= 0) continue;
1463           track->GetClusters(iDet,idx);
1464           AliTrackPoint p;
1465           Int_t isp = 0;
1466           Int_t isp2 = 0;
1467           while (isp < nspdet) {
1468             Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
1469             if (!isvalid) continue;
1470             sp->AddPoint(isptrack,&p); isptrack++; isp++;
1471           }
1472           //      for (Int_t isp = 0; isp < nspdet; isp++) {
1473             //      AliCluster *cl = tracker->GetCluster(idx[isp]);
1474             //      UShort_t volid = tracker->GetVolumeID(idx[isp]);
1475           //        tracker->GetTrackPoint(idx[isp],p);
1476           //        sp->AddPoint(isptrack,&p); isptrack++;
1477           //      }
1478         }       
1479       }
1480     }
1481   fTracker[3]->UnloadClusters();
1482   fLoader[3]->UnloadRecPoints();
1483 }