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