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