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