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