]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
removed warnings
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
CommitLineData
596a855f 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. //
c71de921 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// //
b26c3770 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// //
973388c2 50// In case of raw-data reconstruction the user can modify the default //
51// number of events per digits/clusters/tracks file. In case the option //
52// is not used the number is set 1. In case the user provides 0, than //
53// the number of events is equal to the number of events inside the //
54// raw-data file (i.e. one digits/clusters/tracks file): //
55// //
56// rec.SetNumberOfEventsPerFile(...); //
57// //
58// //
596a855f 59// The name of the galice file can be changed from the default //
e583c30d 60// "galice.root" by passing it as argument to the AliReconstruction //
61// constructor or by //
596a855f 62// //
63// rec.SetGAliceFile("..."); //
64// //
59697224 65// The local reconstruction can be switched on or off for individual //
66// detectors by //
596a855f 67// //
59697224 68// rec.SetRunLocalReconstruction("..."); //
596a855f 69// //
70// The argument is a (case sensitive) string with the names of the //
71// detectors separated by a space. The special string "ALL" selects all //
72// available detectors. This is the default. //
73// //
c71de921 74// The reconstruction of the primary vertex position can be switched off by //
75// //
76// rec.SetRunVertexFinder(kFALSE); //
77// //
b8cd5251 78// The tracking and the creation of ESD tracks can be switched on for //
79// selected detectors by //
596a855f 80// //
b8cd5251 81// rec.SetRunTracking("..."); //
596a855f 82// //
c84a5e9e 83// Uniform/nonuniform field tracking switches (default: uniform field) //
84// //
1d99986f 85// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
c84a5e9e 86// //
596a855f 87// The filling of additional ESD information can be steered by //
88// //
89// rec.SetFillESD("..."); //
90// //
b8cd5251 91// Again, for both methods the string specifies the list of detectors. //
92// The default is "ALL". //
93// //
94// The call of the shortcut method //
95// //
96// rec.SetRunReconstruction("..."); //
97// //
98// is equivalent to calling SetRunLocalReconstruction, SetRunTracking and //
99// SetFillESD with the same detector selecting string as argument. //
596a855f 100// //
c71de921 101// The reconstruction requires digits or raw data as input. For the creation //
102// of digits and raw data have a look at the class AliSimulation. //
596a855f 103// //
24f7a148 104// For debug purposes the method SetCheckPointLevel can be used. If the //
105// argument is greater than 0, files with ESD events will be written after //
106// selected steps of the reconstruction for each event: //
107// level 1: after tracking and after filling of ESD (final) //
108// level 2: in addition after each tracking step //
109// level 3: in addition after the filling of ESD for each detector //
110// If a final check point file exists for an event, this event will be //
111// skipped in the reconstruction. The tracking and the filling of ESD for //
112// a detector will be skipped as well, if the corresponding check point //
113// file exists. The ESD event will then be loaded from the file instead. //
114// //
596a855f 115///////////////////////////////////////////////////////////////////////////////
116
024a7e64 117#include <TArrayF.h>
118#include <TFile.h>
119#include <TSystem.h>
120#include <TROOT.h>
121#include <TPluginManager.h>
fd46e2d2 122#include <TStopwatch.h>
3103d196 123#include <TGeoManager.h>
2bdb9d38 124#include <TLorentzVector.h>
596a855f 125
126#include "AliReconstruction.h"
b8cd5251 127#include "AliReconstructor.h"
815c2b38 128#include "AliLog.h"
596a855f 129#include "AliRunLoader.h"
130#include "AliRun.h"
b649205a 131#include "AliRawReaderFile.h"
132#include "AliRawReaderDate.h"
133#include "AliRawReaderRoot.h"
001397cd 134#include "AliRawEventHeaderBase.h"
596a855f 135#include "AliESD.h"
1d99986f 136#include "AliESDfriend.h"
2257f27e 137#include "AliESDVertex.h"
32e449be 138#include "AliMultiplicity.h"
c84a5e9e 139#include "AliTracker.h"
2257f27e 140#include "AliVertexer.h"
c5e3e5d1 141#include "AliVertexerTracks.h"
5e4ff34d 142#include "AliV0vertexer.h"
143#include "AliCascadeVertexer.h"
596a855f 144#include "AliHeader.h"
145#include "AliGenEventHeader.h"
b26c3770 146#include "AliPID.h"
596a855f 147#include "AliESDpid.h"
ff8bb5ae 148#include "AliESDtrack.h"
f3a97c86 149
150#include "AliRunTag.h"
f3a97c86 151#include "AliDetectorTag.h"
152#include "AliEventTag.h"
153
98937d93 154#include "AliTrackPointArray.h"
b0314964 155#include "AliCDBManager.h"
6bae477a 156#include "AliCDBEntry.h"
157#include "AliAlignObj.h"
f3a97c86 158
b647652d 159#include "AliCentralTrigger.h"
160#include "AliCTPRawStream.h"
161
596a855f 162ClassImp(AliReconstruction)
163
164
c757bafd 165//_____________________________________________________________________________
b384f8a4 166const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
c757bafd 167
596a855f 168//_____________________________________________________________________________
024cf675 169AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
e583c30d 170 const char* name, const char* title) :
171 TNamed(name, title),
172
c84a5e9e 173 fUniformField(kTRUE),
2257f27e 174 fRunVertexFinder(kTRUE),
1f46a9ae 175 fRunHLTTracking(kFALSE),
1d99986f 176 fStopOnError(kFALSE),
177 fWriteAlignmentData(kFALSE),
178 fWriteESDfriend(kFALSE),
b647652d 179 fFillTriggerESD(kTRUE),
1d99986f 180
181 fRunLocalReconstruction("ALL"),
b8cd5251 182 fRunTracking("ALL"),
e583c30d 183 fFillESD("ALL"),
184 fGAliceFileName(gAliceFilename),
b649205a 185 fInput(""),
35042093 186 fEquipIdMap(""),
b26c3770 187 fFirstEvent(0),
188 fLastEvent(-1),
973388c2 189 fNumberOfEventsPerFile(1),
24f7a148 190 fCheckPointLevel(0),
b8cd5251 191 fOptions(),
6bae477a 192 fLoadAlignFromCDB(kTRUE),
193 fLoadAlignData("ALL"),
e583c30d 194
195 fRunLoader(NULL),
b649205a 196 fRawReader(NULL),
b8cd5251 197
98937d93 198 fVertexer(NULL),
9178838a 199 fDiamondProfile(NULL),
98937d93 200
6bae477a 201 fAlignObjArray(NULL),
ec92bee0 202 fCDBUri(cdbUri),
203 fSpecCDBUri()
596a855f 204{
205// create reconstruction object with default parameters
b8cd5251 206
207 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
208 fReconstructor[iDet] = NULL;
209 fLoader[iDet] = NULL;
210 fTracker[iDet] = NULL;
211 }
e47c4c2e 212 AliPID pid;
596a855f 213}
214
215//_____________________________________________________________________________
216AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
e583c30d 217 TNamed(rec),
218
c84a5e9e 219 fUniformField(rec.fUniformField),
2257f27e 220 fRunVertexFinder(rec.fRunVertexFinder),
1f46a9ae 221 fRunHLTTracking(rec.fRunHLTTracking),
1d99986f 222 fStopOnError(rec.fStopOnError),
223 fWriteAlignmentData(rec.fWriteAlignmentData),
224 fWriteESDfriend(rec.fWriteESDfriend),
b647652d 225 fFillTriggerESD(rec.fFillTriggerESD),
1d99986f 226
227 fRunLocalReconstruction(rec.fRunLocalReconstruction),
e583c30d 228 fRunTracking(rec.fRunTracking),
229 fFillESD(rec.fFillESD),
230 fGAliceFileName(rec.fGAliceFileName),
b649205a 231 fInput(rec.fInput),
35042093 232 fEquipIdMap(rec.fEquipIdMap),
b26c3770 233 fFirstEvent(rec.fFirstEvent),
234 fLastEvent(rec.fLastEvent),
973388c2 235 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
24f7a148 236 fCheckPointLevel(0),
b8cd5251 237 fOptions(),
6bae477a 238 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
239 fLoadAlignData(rec.fLoadAlignData),
e583c30d 240
241 fRunLoader(NULL),
b649205a 242 fRawReader(NULL),
b8cd5251 243
98937d93 244 fVertexer(NULL),
9178838a 245 fDiamondProfile(NULL),
98937d93 246
6bae477a 247 fAlignObjArray(rec.fAlignObjArray),
ec92bee0 248 fCDBUri(rec.fCDBUri),
249 fSpecCDBUri()
596a855f 250{
251// copy constructor
252
ec92bee0 253 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
efd2085e 254 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
255 }
b8cd5251 256 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
257 fReconstructor[iDet] = NULL;
258 fLoader[iDet] = NULL;
259 fTracker[iDet] = NULL;
260 }
ec92bee0 261 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
262 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
263 }
596a855f 264}
265
266//_____________________________________________________________________________
267AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
268{
269// assignment operator
270
271 this->~AliReconstruction();
272 new(this) AliReconstruction(rec);
273 return *this;
274}
275
276//_____________________________________________________________________________
277AliReconstruction::~AliReconstruction()
278{
279// clean up
280
e583c30d 281 CleanUp();
efd2085e 282 fOptions.Delete();
ec92bee0 283 fSpecCDBUri.Delete();
596a855f 284}
285
024cf675 286//_____________________________________________________________________________
287void AliReconstruction::InitCDBStorage()
288{
289// activate a default CDB storage
290// First check if we have any CDB storage set, because it is used
291// to retrieve the calibration and alignment constants
292
293 AliCDBManager* man = AliCDBManager::Instance();
ec92bee0 294 if (man->IsDefaultStorageSet())
024cf675 295 {
296 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 297 AliWarning("Default CDB storage has been already set !");
298 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
024cf675 299 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 300 fCDBUri = "";
301 }
302 else {
b8ec52f6 303 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
304 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
305 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
ec92bee0 306 man->SetDefaultStorage(fCDBUri);
307 }
308
309 // Now activate the detector specific CDB storage locations
c3a7b59a 310 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
311 TObject* obj = fSpecCDBUri[i];
312 if (!obj) continue;
b8ec52f6 313 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
314 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
315 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
c3a7b59a 316 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
ec92bee0 317 }
727d0766 318 man->Print();
024cf675 319}
320
321//_____________________________________________________________________________
322void AliReconstruction::SetDefaultStorage(const char* uri) {
ec92bee0 323// Store the desired default CDB storage location
324// Activate it later within the Run() method
024cf675 325
ec92bee0 326 fCDBUri = uri;
024cf675 327
328}
329
330//_____________________________________________________________________________
c3a7b59a 331void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
ec92bee0 332// Store a detector-specific CDB storage location
333// Activate it later within the Run() method
024cf675 334
c3a7b59a 335 AliCDBPath aPath(calibType);
336 if(!aPath.IsValid()){
337 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
338 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
339 if(!strcmp(calibType, fgkDetectorName[iDet])) {
340 aPath.SetPath(Form("%s/*", calibType));
341 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
342 break;
343 }
344 }
345 if(!aPath.IsValid()){
346 AliError(Form("Not a valid path or detector: %s", calibType));
347 return;
348 }
349 }
350
351 // check that calibType refers to a "valid" detector name
352 Bool_t isDetector = kFALSE;
353 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
354 TString detName = fgkDetectorName[iDet];
355 if(aPath.GetLevel0() == detName) {
356 isDetector = kTRUE;
357 break;
358 }
359 }
360
361 if(!isDetector) {
362 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
363 return;
364 }
365
366 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
ec92bee0 367 if (obj) fSpecCDBUri.Remove(obj);
c3a7b59a 368 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
024cf675 369
370}
371
6bae477a 372//_____________________________________________________________________________
373Bool_t AliReconstruction::SetRunNumber()
374{
375 // The method is called in Run() in order
376 // to set a correct run number.
377 // In case of raw data reconstruction the
378 // run number is taken from the raw data header
379
380 if(AliCDBManager::Instance()->GetRun() < 0) {
381 if (!fRunLoader) {
382 AliError("No run loader is found !");
383 return kFALSE;
384 }
385 // read run number from gAlice
ec92bee0 386 if(fRunLoader->GetAliRun())
387 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
388 else {
389 if(fRawReader) {
390 if(fRawReader->NextEvent()) {
391 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
392 fRawReader->RewindEvents();
393 }
394 else {
395 AliError("No raw-data events found !");
396 return kFALSE;
397 }
398 }
399 else {
400 AliError("Neither gAlice nor RawReader objects are found !");
401 return kFALSE;
402 }
403 }
404 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
6bae477a 405 }
406 return kTRUE;
407}
408
409//_____________________________________________________________________________
410Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
411{
412 // Read collection of alignment objects (AliAlignObj derived) saved
413 // in the TClonesArray ClArrayName and apply them to the geometry
414 // manager singleton.
415 //
416 alObjArray->Sort();
417 Int_t nvols = alObjArray->GetEntriesFast();
418
dc0984f8 419 Bool_t flag = kTRUE;
420
6bae477a 421 for(Int_t j=0; j<nvols; j++)
422 {
423 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
dc0984f8 424 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
6bae477a 425 }
426
427 if (AliDebugLevelClass() >= 1) {
e0625db4 428 gGeoManager->GetTopNode()->CheckOverlaps(20);
6bae477a 429 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
430 if(ovexlist->GetEntriesFast()){
431 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
432 }
433 }
434
dc0984f8 435 return flag;
6bae477a 436
437}
438
439//_____________________________________________________________________________
440Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
441{
442 // Fills array of single detector's alignable objects from CDB
443
444 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
445
446 AliCDBEntry *entry;
447
448 AliCDBPath path(detName,"Align","Data");
449
450 entry=AliCDBManager::Instance()->Get(path.GetPath());
451 if(!entry){
452 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
453 return kFALSE;
454 }
455 entry->SetOwner(1);
456 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
457 alignArray->SetOwner(0);
458 AliDebug(2,Form("Found %d alignment objects for %s",
459 alignArray->GetEntries(),detName));
460
461 AliAlignObj *alignObj=0;
462 TIter iter(alignArray);
463
464 // loop over align objects in detector
465 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
466 fAlignObjArray->Add(alignObj);
467 }
468 // delete entry --- Don't delete, it is cached!
469
470 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
471 return kTRUE;
472
473}
474
475//_____________________________________________________________________________
476Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
477{
478 // Read the alignment objects from CDB.
479 // Each detector is supposed to have the
480 // alignment objects in DET/Align/Data CDB path.
481 // All the detector objects are then collected,
482 // sorted by geometry level (starting from ALIC) and
483 // then applied to the TGeo geometry.
484 // Finally an overlaps check is performed.
485
486 // Load alignment data from CDB and fill fAlignObjArray
487 if(fLoadAlignFromCDB){
488 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
489
490 //fAlignObjArray->RemoveAll();
491 fAlignObjArray->Clear();
492 fAlignObjArray->SetOwner(0);
493
494 TString detStr = detectors;
495 TString dataNotLoaded="";
496 TString dataLoaded="";
497
498 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
499 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
500 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
501 dataNotLoaded += fgkDetectorName[iDet];
502 dataNotLoaded += " ";
503 } else {
504 dataLoaded += fgkDetectorName[iDet];
505 dataLoaded += " ";
506 }
507 } // end loop over detectors
508
509 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
510 dataNotLoaded += detStr;
7250c343 511 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
6bae477a 512 dataLoaded.Data()));
7250c343 513 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
6bae477a 514 dataNotLoaded.Data()));
515 } // fLoadAlignFromCDB flag
516
517 // Check if the array with alignment objects was
518 // provided by the user. If yes, apply the objects
519 // to the present TGeo geometry
520 if (fAlignObjArray) {
521 if (gGeoManager && gGeoManager->IsClosed()) {
522 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
dc0984f8 523 AliError("The misalignment of one or more volumes failed!"
524 "Compare the list of simulated detectors and the list of detector alignment data!");
6bae477a 525 return kFALSE;
526 }
527 }
528 else {
529 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
530 return kFALSE;
531 }
532 }
533
8e245d15 534 delete fAlignObjArray; fAlignObjArray=0;
6bae477a 535 return kTRUE;
536}
596a855f 537
538//_____________________________________________________________________________
539void AliReconstruction::SetGAliceFile(const char* fileName)
540{
541// set the name of the galice file
542
543 fGAliceFileName = fileName;
544}
545
efd2085e 546//_____________________________________________________________________________
547void AliReconstruction::SetOption(const char* detector, const char* option)
548{
549// set options for the reconstruction of a detector
550
551 TObject* obj = fOptions.FindObject(detector);
552 if (obj) fOptions.Remove(obj);
553 fOptions.Add(new TNamed(detector, option));
554}
555
596a855f 556
557//_____________________________________________________________________________
4a33489c 558Bool_t AliReconstruction::Run(const char* input)
596a855f 559{
560// run the reconstruction
561
b649205a 562 // set the input
563 if (!input) input = fInput.Data();
564 TString fileName(input);
565 if (fileName.EndsWith("/")) {
566 fRawReader = new AliRawReaderFile(fileName);
567 } else if (fileName.EndsWith(".root")) {
568 fRawReader = new AliRawReaderRoot(fileName);
569 } else if (!fileName.IsNull()) {
570 fRawReader = new AliRawReaderDate(fileName);
571 fRawReader->SelectEvents(7);
572 }
35042093 573 if (!fEquipIdMap.IsNull() && fRawReader)
574 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
575
b649205a 576
f08fc9f5 577 // get the run loader
578 if (!InitRunLoader()) return kFALSE;
596a855f 579
ec92bee0 580 // Initialize the CDB storage
581 InitCDBStorage();
582
6bae477a 583 // Set run number in CDBManager (if it is not already set by the user)
584 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
585
586 // Import ideal TGeo geometry and apply misalignment
587 if (!gGeoManager) {
588 TString geom(gSystem->DirName(fGAliceFileName));
589 geom += "/geometry.root";
590 TGeoManager::Import(geom.Data());
591 if (!gGeoManager) if (fStopOnError) return kFALSE;
592 }
8e245d15 593
594 AliCDBManager* man = AliCDBManager::Instance();
6bae477a 595 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
596
596a855f 597 // local reconstruction
59697224 598 if (!fRunLocalReconstruction.IsNull()) {
599 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
e583c30d 600 if (fStopOnError) {CleanUp(); return kFALSE;}
596a855f 601 }
602 }
b26c3770 603// if (!fRunVertexFinder && fRunTracking.IsNull() &&
604// fFillESD.IsNull()) return kTRUE;
2257f27e 605
606 // get vertexer
607 if (fRunVertexFinder && !CreateVertexer()) {
608 if (fStopOnError) {
609 CleanUp();
610 return kFALSE;
611 }
612 }
596a855f 613
f08fc9f5 614 // get trackers
b8cd5251 615 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
24f7a148 616 if (fStopOnError) {
617 CleanUp();
618 return kFALSE;
619 }
596a855f 620 }
24f7a148 621
9db3a215 622
623 TStopwatch stopwatch;
624 stopwatch.Start();
625
b26c3770 626 // get the possibly already existing ESD file and tree
1f46a9ae 627 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
b26c3770 628 TFile* fileOld = NULL;
1f46a9ae 629 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
b26c3770 630 if (!gSystem->AccessPathName("AliESDs.root")){
631 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
632 fileOld = TFile::Open("AliESDs.old.root");
633 if (fileOld && fileOld->IsOpen()) {
634 treeOld = (TTree*) fileOld->Get("esdTree");
635 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
1f46a9ae 636 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
637 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
b26c3770 638 }
639 }
640
36711aa4 641 // create the ESD output file and tree
596a855f 642 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
643 if (!file->IsOpen()) {
815c2b38 644 AliError("opening AliESDs.root failed");
b26c3770 645 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 646 }
36711aa4 647 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
648 tree->Branch("ESD", "AliESD", &esd);
1f46a9ae 649 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
650 hlttree->Branch("ESD", "AliESD", &hltesd);
651 delete esd; delete hltesd;
652 esd = NULL; hltesd = NULL;
1d99986f 653
500d54ab 654 // create the branch with ESD additions
655 AliESDfriend *esdf=0;
1d99986f 656 if (fWriteESDfriend) {
500d54ab 657 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
658 br->SetFile("AliESDfriends.root");
1d99986f 659 }
660
c5e3e5d1 661 AliVertexerTracks tVertexer;
9178838a 662 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
c5e3e5d1 663
596a855f 664 // loop over events
b649205a 665 if (fRawReader) fRawReader->RewindEvents();
f08fc9f5 666
596a855f 667 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
b26c3770 668 if (fRawReader) fRawReader->NextEvent();
4a33489c 669 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
b26c3770 670 // copy old ESD to the new one
671 if (treeOld) {
672 treeOld->SetBranchAddress("ESD", &esd);
673 treeOld->GetEntry(iEvent);
674 }
675 tree->Fill();
1f46a9ae 676 if (hlttreeOld) {
677 hlttreeOld->SetBranchAddress("ESD", &hltesd);
678 hlttreeOld->GetEntry(iEvent);
679 }
680 hlttree->Fill();
b26c3770 681 continue;
682 }
683
815c2b38 684 AliInfo(Form("processing event %d", iEvent));
596a855f 685 fRunLoader->GetEvent(iEvent);
24f7a148 686
bb0901a4 687 char aFileName[256];
688 sprintf(aFileName, "ESD_%d.%d_final.root",
f08fc9f5 689 fRunLoader->GetHeader()->GetRun(),
690 fRunLoader->GetHeader()->GetEventNrInRun());
bb0901a4 691 if (!gSystem->AccessPathName(aFileName)) continue;
24f7a148 692
b26c3770 693 // local reconstruction
694 if (!fRunLocalReconstruction.IsNull()) {
695 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
696 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
697 }
698 }
699
1f46a9ae 700 esd = new AliESD; hltesd = new AliESD;
f08fc9f5 701 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
1f46a9ae 702 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
f08fc9f5 703 esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
1f46a9ae 704 hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
d6ee376f 705
706 // Set magnetic field from the tracker
707 esd->SetMagneticField(AliTracker::GetBz());
708 hltesd->SetMagneticField(AliTracker::GetBz());
596a855f 709
2257f27e 710 // vertex finder
711 if (fRunVertexFinder) {
712 if (!ReadESD(esd, "vertex")) {
713 if (!RunVertexFinder(esd)) {
b26c3770 714 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
2257f27e 715 }
716 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
717 }
718 }
719
1f46a9ae 720 // HLT tracking
721 if (!fRunTracking.IsNull()) {
722 if (fRunHLTTracking) {
723 hltesd->SetVertex(esd->GetVertex());
724 if (!RunHLTTracking(hltesd)) {
725 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
726 }
727 }
728 }
729
596a855f 730 // barrel tracking
b8cd5251 731 if (!fRunTracking.IsNull()) {
24f7a148 732 if (!ReadESD(esd, "tracking")) {
733 if (!RunTracking(esd)) {
b26c3770 734 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
24f7a148 735 }
736 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
596a855f 737 }
738 }
739
740 // fill ESD
741 if (!fFillESD.IsNull()) {
742 if (!FillESD(esd, fFillESD)) {
b26c3770 743 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
596a855f 744 }
745 }
001397cd 746 // fill Event header information from the RawEventHeader
747 if (fRawReader){FillRawEventHeaderESD(esd);}
596a855f 748
749 // combined PID
750 AliESDpid::MakePID(esd);
24f7a148 751 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
596a855f 752
b647652d 753 if (fFillTriggerESD) {
754 if (!ReadESD(esd, "trigger")) {
755 if (!FillTriggerESD(esd)) {
756 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
757 }
758 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
759 }
760 }
761
a03cd2e0 762 //Try to improve the reconstructed primary vertex position using the tracks
763 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
764 if (pvtx)
765 if (pvtx->GetStatus()) {
766 // Store the improved primary vertex
767 esd->SetPrimaryVertex(pvtx);
768 // Propagate the tracks to the DCA to the improved primary vertex
769 Double_t somethingbig = 777.;
770 Double_t bz = esd->GetMagneticField();
771 Int_t nt=esd->GetNumberOfTracks();
772 while (nt--) {
773 AliESDtrack *t = esd->GetTrack(nt);
774 t->RelateToVertex(pvtx, bz, somethingbig);
775 }
776 }
c5e3e5d1 777
5e4ff34d 778 {
779 // V0 finding
780 AliV0vertexer vtxer;
781 vtxer.Tracks2V0vertices(esd);
782
783 // Cascade finding
784 AliCascadeVertexer cvtxer;
785 cvtxer.V0sTracks2CascadeVertices(esd);
786 }
787
596a855f 788 // write ESD
1d99986f 789 if (fWriteESDfriend) {
d75007f6 790 esdf=new AliESDfriend();
791 esd->GetESDfriend(esdf);
1d99986f 792 }
500d54ab 793 tree->Fill();
794
795 // write HLT ESD
796 hlttree->Fill();
1d99986f 797
f3a97c86 798 if (fCheckPointLevel > 0) WriteESD(esd, "final");
799
500d54ab 800 delete esd; delete esdf; delete hltesd;
801 esd = NULL; esdf=NULL; hltesd = NULL;
596a855f 802 }
803
9db3a215 804 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
805 stopwatch.RealTime(),stopwatch.CpuTime()));
806
36711aa4 807 file->cd();
a9c0e6db 808 if (fWriteESDfriend)
809 tree->SetBranchStatus("ESDfriend*",0);
36711aa4 810 tree->Write();
1f46a9ae 811 hlttree->Write();
f3a97c86 812
813 // Create tags for the events in the ESD tree (the ESD tree is always present)
814 // In case of empty events the tags will contain dummy values
815 CreateTag(file);
b26c3770 816 CleanUp(file, fileOld);
596a855f 817
818 return kTRUE;
819}
820
821
822//_____________________________________________________________________________
59697224 823Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
596a855f 824{
59697224 825// run the local reconstruction
596a855f 826
030b532d 827 TStopwatch stopwatch;
828 stopwatch.Start();
829
8e245d15 830 AliCDBManager* man = AliCDBManager::Instance();
831 Bool_t origCache = man->GetCacheFlag();
832
596a855f 833 TString detStr = detectors;
b8cd5251 834 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
835 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
836 AliReconstructor* reconstructor = GetReconstructor(iDet);
837 if (!reconstructor) continue;
b26c3770 838 if (reconstructor->HasLocalReconstruction()) continue;
b8cd5251 839
840 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
841 TStopwatch stopwatchDet;
842 stopwatchDet.Start();
8e245d15 843
844 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
845
846 man->SetCacheFlag(kTRUE);
847 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
848 man->GetAll(calibPath); // entries are cached!
849
b8cd5251 850 if (fRawReader) {
851 fRawReader->RewindEvents();
852 reconstructor->Reconstruct(fRunLoader, fRawReader);
853 } else {
854 reconstructor->Reconstruct(fRunLoader);
596a855f 855 }
5f8272e1 856 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
857 fgkDetectorName[iDet],
858 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
8e245d15 859
860 // unload calibration data
861 man->ClearCache();
596a855f 862 }
863
8e245d15 864 man->SetCacheFlag(origCache);
865
596a855f 866 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 867 AliError(Form("the following detectors were not found: %s",
868 detStr.Data()));
596a855f 869 if (fStopOnError) return kFALSE;
870 }
871
5f8272e1 872 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
873 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 874
596a855f 875 return kTRUE;
876}
877
b26c3770 878//_____________________________________________________________________________
879Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
880{
881// run the local reconstruction
882
883 TStopwatch stopwatch;
884 stopwatch.Start();
885
886 TString detStr = detectors;
887 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
888 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
889 AliReconstructor* reconstructor = GetReconstructor(iDet);
890 if (!reconstructor) continue;
891 AliLoader* loader = fLoader[iDet];
892
893 // conversion of digits
894 if (fRawReader && reconstructor->HasDigitConversion()) {
895 AliInfo(Form("converting raw data digits into root objects for %s",
896 fgkDetectorName[iDet]));
897 TStopwatch stopwatchDet;
898 stopwatchDet.Start();
899 loader->LoadDigits("update");
900 loader->CleanDigits();
901 loader->MakeDigitsContainer();
902 TTree* digitsTree = loader->TreeD();
903 reconstructor->ConvertDigits(fRawReader, digitsTree);
904 loader->WriteDigits("OVERWRITE");
905 loader->UnloadDigits();
5f8272e1 906 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
907 fgkDetectorName[iDet],
908 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 909 }
910
911 // local reconstruction
912 if (!reconstructor->HasLocalReconstruction()) continue;
913 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
914 TStopwatch stopwatchDet;
915 stopwatchDet.Start();
916 loader->LoadRecPoints("update");
917 loader->CleanRecPoints();
918 loader->MakeRecPointsContainer();
919 TTree* clustersTree = loader->TreeR();
920 if (fRawReader && !reconstructor->HasDigitConversion()) {
921 reconstructor->Reconstruct(fRawReader, clustersTree);
922 } else {
923 loader->LoadDigits("read");
924 TTree* digitsTree = loader->TreeD();
925 if (!digitsTree) {
926 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
927 if (fStopOnError) return kFALSE;
928 } else {
929 reconstructor->Reconstruct(digitsTree, clustersTree);
930 }
931 loader->UnloadDigits();
932 }
933 loader->WriteRecPoints("OVERWRITE");
934 loader->UnloadRecPoints();
5f8272e1 935 AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
936 fgkDetectorName[iDet],
937 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
b26c3770 938 }
939
940 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
941 AliError(Form("the following detectors were not found: %s",
942 detStr.Data()));
943 if (fStopOnError) return kFALSE;
944 }
5f8272e1 945
946 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
947 stopwatch.RealTime(),stopwatch.CpuTime()));
b26c3770 948
949 return kTRUE;
950}
951
596a855f 952//_____________________________________________________________________________
2257f27e 953Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
596a855f 954{
955// run the barrel tracking
956
030b532d 957 TStopwatch stopwatch;
958 stopwatch.Start();
959
2257f27e 960 AliESDVertex* vertex = NULL;
961 Double_t vtxPos[3] = {0, 0, 0};
962 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
963 TArrayF mcVertex(3);
a6b0b91b 964 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
965 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
966 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
967 }
2257f27e 968
b8cd5251 969 if (fVertexer) {
815c2b38 970 AliInfo("running the ITS vertex finder");
b26c3770 971 if (fLoader[0]) fLoader[0]->LoadRecPoints();
b8cd5251 972 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
b26c3770 973 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
2257f27e 974 if(!vertex){
815c2b38 975 AliWarning("Vertex not found");
c710f220 976 vertex = new AliESDVertex();
d1a50cb5 977 vertex->SetName("default");
2257f27e 978 }
979 else {
980 vertex->SetTruePos(vtxPos); // store also the vertex from MC
d1a50cb5 981 vertex->SetName("reconstructed");
2257f27e 982 }
983
984 } else {
815c2b38 985 AliInfo("getting the primary vertex from MC");
2257f27e 986 vertex = new AliESDVertex(vtxPos, vtxErr);
987 }
988
989 if (vertex) {
990 vertex->GetXYZ(vtxPos);
991 vertex->GetSigmaXYZ(vtxErr);
992 } else {
815c2b38 993 AliWarning("no vertex reconstructed");
2257f27e 994 vertex = new AliESDVertex(vtxPos, vtxErr);
995 }
996 esd->SetVertex(vertex);
32e449be 997 // if SPD multiplicity has been determined, it is stored in the ESD
998 AliMultiplicity *mult= fVertexer->GetMultiplicity();
999 if(mult)esd->SetMultiplicity(mult);
1000
b8cd5251 1001 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1002 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1003 }
2257f27e 1004 delete vertex;
1005
5f8272e1 1006 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1007 stopwatch.RealTime(),stopwatch.CpuTime()));
2257f27e 1008
1009 return kTRUE;
1010}
1011
1f46a9ae 1012//_____________________________________________________________________________
1013Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1014{
1015// run the HLT barrel tracking
1016
1017 TStopwatch stopwatch;
1018 stopwatch.Start();
1019
1020 if (!fRunLoader) {
1021 AliError("Missing runLoader!");
1022 return kFALSE;
1023 }
1024
1025 AliInfo("running HLT tracking");
1026
1027 // Get a pointer to the HLT reconstructor
1028 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1029 if (!reconstructor) return kFALSE;
1030
1031 // TPC + ITS
1032 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1033 TString detName = fgkDetectorName[iDet];
1034 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1035 reconstructor->SetOption(detName.Data());
1036 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1037 if (!tracker) {
1038 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1039 if (fStopOnError) return kFALSE;
9dcc06e1 1040 continue;
1f46a9ae 1041 }
1042 Double_t vtxPos[3];
1043 Double_t vtxErr[3]={0.005,0.005,0.010};
1044 const AliESDVertex *vertex = esd->GetVertex();
1045 vertex->GetXYZ(vtxPos);
1046 tracker->SetVertex(vtxPos,vtxErr);
1047 if(iDet != 1) {
1048 fLoader[iDet]->LoadRecPoints("read");
1049 TTree* tree = fLoader[iDet]->TreeR();
1050 if (!tree) {
1051 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1052 return kFALSE;
1053 }
1054 tracker->LoadClusters(tree);
1055 }
1056 if (tracker->Clusters2Tracks(esd) != 0) {
1057 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1058 return kFALSE;
1059 }
1060 if(iDet != 1) {
1061 tracker->UnloadClusters();
1062 }
1063 delete tracker;
1064 }
1065
5f8272e1 1066 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1067 stopwatch.RealTime(),stopwatch.CpuTime()));
1f46a9ae 1068
1069 return kTRUE;
1070}
1071
2257f27e 1072//_____________________________________________________________________________
1073Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1074{
1075// run the barrel tracking
1076
1077 TStopwatch stopwatch;
1078 stopwatch.Start();
24f7a148 1079
815c2b38 1080 AliInfo("running tracking");
596a855f 1081
b8cd5251 1082 // pass 1: TPC + ITS inwards
1083 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1084 if (!fTracker[iDet]) continue;
1085 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
24f7a148 1086
b8cd5251 1087 // load clusters
1088 fLoader[iDet]->LoadRecPoints("read");
1089 TTree* tree = fLoader[iDet]->TreeR();
1090 if (!tree) {
1091 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1092 return kFALSE;
1093 }
b8cd5251 1094 fTracker[iDet]->LoadClusters(tree);
1095
1096 // run tracking
1097 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1098 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
24f7a148 1099 return kFALSE;
1100 }
b8cd5251 1101 if (fCheckPointLevel > 1) {
1102 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1103 }
878e1fe1 1104 // preliminary PID in TPC needed by the ITS tracker
1105 if (iDet == 1) {
1106 GetReconstructor(1)->FillESD(fRunLoader, esd);
b26c3770 1107 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
878e1fe1 1108 AliESDpid::MakePID(esd);
1109 }
b8cd5251 1110 }
596a855f 1111
b8cd5251 1112 // pass 2: ALL backwards
1113 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1114 if (!fTracker[iDet]) continue;
1115 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1116
1117 // load clusters
1118 if (iDet > 1) { // all except ITS, TPC
1119 TTree* tree = NULL;
7b61cd9c 1120 fLoader[iDet]->LoadRecPoints("read");
1121 tree = fLoader[iDet]->TreeR();
b8cd5251 1122 if (!tree) {
1123 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
24f7a148 1124 return kFALSE;
1125 }
b8cd5251 1126 fTracker[iDet]->LoadClusters(tree);
1127 }
24f7a148 1128
b8cd5251 1129 // run tracking
1130 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1131 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1132 return kFALSE;
1133 }
1134 if (fCheckPointLevel > 1) {
1135 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1136 }
24f7a148 1137
b8cd5251 1138 // unload clusters
1139 if (iDet > 2) { // all except ITS, TPC, TRD
1140 fTracker[iDet]->UnloadClusters();
7b61cd9c 1141 fLoader[iDet]->UnloadRecPoints();
b8cd5251 1142 }
8f37df88 1143 // updated PID in TPC needed by the ITS tracker -MI
1144 if (iDet == 1) {
1145 GetReconstructor(1)->FillESD(fRunLoader, esd);
1146 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1147 AliESDpid::MakePID(esd);
1148 }
b8cd5251 1149 }
596a855f 1150
98937d93 1151 // write space-points to the ESD in case alignment data output
1152 // is switched on
1153 if (fWriteAlignmentData)
1154 WriteAlignmentData(esd);
1155
b8cd5251 1156 // pass 3: TRD + TPC + ITS refit inwards
1157 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1158 if (!fTracker[iDet]) continue;
1159 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
596a855f 1160
b8cd5251 1161 // run tracking
1162 if (fTracker[iDet]->RefitInward(esd) != 0) {
1163 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1164 return kFALSE;
1165 }
1166 if (fCheckPointLevel > 1) {
1167 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1168 }
596a855f 1169
b8cd5251 1170 // unload clusters
1171 fTracker[iDet]->UnloadClusters();
1172 fLoader[iDet]->UnloadRecPoints();
1173 }
ff8bb5ae 1174 //
1175 // Propagate track to the vertex - if not done by ITS
1176 //
1177 Int_t ntracks = esd->GetNumberOfTracks();
1178 for (Int_t itrack=0; itrack<ntracks; itrack++){
1179 const Double_t kRadius = 3; // beam pipe radius
1180 const Double_t kMaxStep = 5; // max step
1181 const Double_t kMaxD = 123456; // max distance to prim vertex
1182 Double_t fieldZ = AliTracker::GetBz(); //
1183 AliESDtrack * track = esd->GetTrack(itrack);
1184 if (!track) continue;
1185 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
a28195f7 1186 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
ff8bb5ae 1187 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1188 }
1189
5f8272e1 1190 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1191 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1192
596a855f 1193 return kTRUE;
1194}
1195
1196//_____________________________________________________________________________
24f7a148 1197Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
596a855f 1198{
1199// fill the event summary data
1200
030b532d 1201 TStopwatch stopwatch;
1202 stopwatch.Start();
815c2b38 1203 AliInfo("filling ESD");
030b532d 1204
596a855f 1205 TString detStr = detectors;
b8cd5251 1206 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1207 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1208 AliReconstructor* reconstructor = GetReconstructor(iDet);
1209 if (!reconstructor) continue;
1210
1211 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1212 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
b26c3770 1213 TTree* clustersTree = NULL;
1214 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1215 fLoader[iDet]->LoadRecPoints("read");
1216 clustersTree = fLoader[iDet]->TreeR();
1217 if (!clustersTree) {
1218 AliError(Form("Can't get the %s clusters tree",
1219 fgkDetectorName[iDet]));
1220 if (fStopOnError) return kFALSE;
1221 }
1222 }
1223 if (fRawReader && !reconstructor->HasDigitConversion()) {
1224 reconstructor->FillESD(fRawReader, clustersTree, esd);
1225 } else {
1226 TTree* digitsTree = NULL;
1227 if (fLoader[iDet]) {
1228 fLoader[iDet]->LoadDigits("read");
1229 digitsTree = fLoader[iDet]->TreeD();
1230 if (!digitsTree) {
1231 AliError(Form("Can't get the %s digits tree",
1232 fgkDetectorName[iDet]));
1233 if (fStopOnError) return kFALSE;
1234 }
1235 }
1236 reconstructor->FillESD(digitsTree, clustersTree, esd);
1237 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1238 }
1239 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1240 fLoader[iDet]->UnloadRecPoints();
1241 }
1242
b8cd5251 1243 if (fRawReader) {
1244 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1245 } else {
1246 reconstructor->FillESD(fRunLoader, esd);
24f7a148 1247 }
b8cd5251 1248 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
596a855f 1249 }
1250 }
1251
1252 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
815c2b38 1253 AliError(Form("the following detectors were not found: %s",
1254 detStr.Data()));
596a855f 1255 if (fStopOnError) return kFALSE;
1256 }
1257
5f8272e1 1258 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1259 stopwatch.RealTime(),stopwatch.CpuTime()));
030b532d 1260
596a855f 1261 return kTRUE;
1262}
1263
b647652d 1264//_____________________________________________________________________________
1265Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1266{
1267 // Reads the trigger decision which is
1268 // stored in Trigger.root file and fills
1269 // the corresponding esd entries
1270
1271 AliInfo("Filling trigger information into the ESD");
1272
1273 if (fRawReader) {
1274 AliCTPRawStream input(fRawReader);
1275 if (!input.Next()) {
1276 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1277 return kFALSE;
1278 }
1279 esd->SetTriggerMask(input.GetClassMask());
1280 esd->SetTriggerCluster(input.GetClusterMask());
1281 }
1282 else {
1283 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1284 if (runloader) {
1285 if (!runloader->LoadTrigger()) {
1286 AliCentralTrigger *aCTP = runloader->GetTrigger();
1287 esd->SetTriggerMask(aCTP->GetClassMask());
1288 esd->SetTriggerCluster(aCTP->GetClusterMask());
1289 }
1290 else {
1291 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1292 return kFALSE;
1293 }
1294 }
1295 else {
1296 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1297 return kFALSE;
1298 }
1299 }
1300
1301 return kTRUE;
1302}
596a855f 1303
001397cd 1304
1305
1306
1307
1308//_____________________________________________________________________________
1309Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1310{
1311 //
1312 // Filling information from RawReader Header
1313 //
1314
1315 AliInfo("Filling information from RawReader Header");
1316 esd->SetTimeStamp(0);
1317 esd->SetEventType(0);
1318 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1319 if (eventHeader){
1320 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1321 esd->SetEventType((eventHeader->Get("Type")));
1322 }
1323
1324 return kTRUE;
1325}
1326
1327
596a855f 1328//_____________________________________________________________________________
1329Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1330{
1331// check whether detName is contained in detectors
1332// if yes, it is removed from detectors
1333
1334 // check if all detectors are selected
1335 if ((detectors.CompareTo("ALL") == 0) ||
1336 detectors.BeginsWith("ALL ") ||
1337 detectors.EndsWith(" ALL") ||
1338 detectors.Contains(" ALL ")) {
1339 detectors = "ALL";
1340 return kTRUE;
1341 }
1342
1343 // search for the given detector
1344 Bool_t result = kFALSE;
1345 if ((detectors.CompareTo(detName) == 0) ||
1346 detectors.BeginsWith(detName+" ") ||
1347 detectors.EndsWith(" "+detName) ||
1348 detectors.Contains(" "+detName+" ")) {
1349 detectors.ReplaceAll(detName, "");
1350 result = kTRUE;
1351 }
1352
1353 // clean up the detectors string
1354 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1355 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1356 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1357
1358 return result;
1359}
e583c30d 1360
f08fc9f5 1361//_____________________________________________________________________________
1362Bool_t AliReconstruction::InitRunLoader()
1363{
1364// get or create the run loader
1365
1366 if (gAlice) delete gAlice;
1367 gAlice = NULL;
1368
b26c3770 1369 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1370 // load all base libraries to get the loader classes
1371 TString libs = gSystem->GetLibraries();
1372 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1373 TString detName = fgkDetectorName[iDet];
1374 if (detName == "HLT") continue;
1375 if (libs.Contains("lib" + detName + "base.so")) continue;
1376 gSystem->Load("lib" + detName + "base.so");
1377 }
f08fc9f5 1378 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1379 if (!fRunLoader) {
1380 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1381 CleanUp();
1382 return kFALSE;
1383 }
b26c3770 1384 fRunLoader->CdGAFile();
1385 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1386 if (fRunLoader->LoadgAlice() == 0) {
1387 gAlice = fRunLoader->GetAliRun();
c84a5e9e 1388 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
b26c3770 1389 }
f08fc9f5 1390 }
1391 if (!gAlice && !fRawReader) {
1392 AliError(Form("no gAlice object found in file %s",
1393 fGAliceFileName.Data()));
1394 CleanUp();
1395 return kFALSE;
1396 }
1397
1398 } else { // galice.root does not exist
1399 if (!fRawReader) {
1400 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1401 CleanUp();
1402 return kFALSE;
1403 }
1404 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1405 AliConfig::GetDefaultEventFolderName(),
1406 "recreate");
1407 if (!fRunLoader) {
1408 AliError(Form("could not create run loader in file %s",
1409 fGAliceFileName.Data()));
1410 CleanUp();
1411 return kFALSE;
1412 }
1413 fRunLoader->MakeTree("E");
1414 Int_t iEvent = 0;
1415 while (fRawReader->NextEvent()) {
1416 fRunLoader->SetEventNumber(iEvent);
1417 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1418 iEvent, iEvent);
1419 fRunLoader->MakeTree("H");
1420 fRunLoader->TreeE()->Fill();
1421 iEvent++;
1422 }
1423 fRawReader->RewindEvents();
973388c2 1424 if (fNumberOfEventsPerFile > 0)
1425 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1426 else
1427 fRunLoader->SetNumberOfEventsPerFile(iEvent);
f08fc9f5 1428 fRunLoader->WriteHeader("OVERWRITE");
1429 fRunLoader->CdGAFile();
1430 fRunLoader->Write(0, TObject::kOverwrite);
1431// AliTracker::SetFieldMap(???);
1432 }
1433
1434 return kTRUE;
1435}
1436
c757bafd 1437//_____________________________________________________________________________
b8cd5251 1438AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
c757bafd 1439{
f08fc9f5 1440// get the reconstructor object and the loader for a detector
c757bafd 1441
b8cd5251 1442 if (fReconstructor[iDet]) return fReconstructor[iDet];
1443
1444 // load the reconstructor object
1445 TPluginManager* pluginManager = gROOT->GetPluginManager();
1446 TString detName = fgkDetectorName[iDet];
1447 TString recName = "Ali" + detName + "Reconstructor";
f08fc9f5 1448 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
b8cd5251 1449
1450 if (detName == "HLT") {
1451 if (!gROOT->GetClass("AliLevel3")) {
4aa41877 1452 gSystem->Load("libAliHLTSrc.so");
1453 gSystem->Load("libAliHLTMisc.so");
1454 gSystem->Load("libAliHLTHough.so");
1455 gSystem->Load("libAliHLTComp.so");
b8cd5251 1456 }
1457 }
1458
1459 AliReconstructor* reconstructor = NULL;
1460 // first check if a plugin is defined for the reconstructor
1461 TPluginHandler* pluginHandler =
1462 pluginManager->FindHandler("AliReconstructor", detName);
f08fc9f5 1463 // if not, add a plugin for it
1464 if (!pluginHandler) {
b8cd5251 1465 AliDebug(1, Form("defining plugin for %s", recName.Data()));
b26c3770 1466 TString libs = gSystem->GetLibraries();
1467 if (libs.Contains("lib" + detName + "base.so") ||
1468 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
b8cd5251 1469 pluginManager->AddHandler("AliReconstructor", detName,
1470 recName, detName + "rec", recName + "()");
1471 } else {
1472 pluginManager->AddHandler("AliReconstructor", detName,
1473 recName, detName, recName + "()");
c757bafd 1474 }
b8cd5251 1475 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1476 }
1477 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1478 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
c757bafd 1479 }
b8cd5251 1480 if (reconstructor) {
1481 TObject* obj = fOptions.FindObject(detName.Data());
1482 if (obj) reconstructor->SetOption(obj->GetTitle());
b26c3770 1483 reconstructor->Init(fRunLoader);
b8cd5251 1484 fReconstructor[iDet] = reconstructor;
1485 }
1486
f08fc9f5 1487 // get or create the loader
1488 if (detName != "HLT") {
1489 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1490 if (!fLoader[iDet]) {
1491 AliConfig::Instance()
1492 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1493 detName, detName);
1494 // first check if a plugin is defined for the loader
bb0901a4 1495 pluginHandler =
f08fc9f5 1496 pluginManager->FindHandler("AliLoader", detName);
1497 // if not, add a plugin for it
1498 if (!pluginHandler) {
1499 TString loaderName = "Ali" + detName + "Loader";
1500 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1501 pluginManager->AddHandler("AliLoader", detName,
1502 loaderName, detName + "base",
1503 loaderName + "(const char*, TFolder*)");
1504 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1505 }
1506 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1507 fLoader[iDet] =
1508 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1509 fRunLoader->GetEventFolder());
1510 }
1511 if (!fLoader[iDet]) { // use default loader
1512 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1513 }
1514 if (!fLoader[iDet]) {
1515 AliWarning(Form("couldn't get loader for %s", detName.Data()));
6667b602 1516 if (fStopOnError) return NULL;
f08fc9f5 1517 } else {
1518 fRunLoader->AddLoader(fLoader[iDet]);
1519 fRunLoader->CdGAFile();
1520 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1521 fRunLoader->Write(0, TObject::kOverwrite);
1522 }
1523 }
1524 }
1525
b8cd5251 1526 return reconstructor;
c757bafd 1527}
1528
2257f27e 1529//_____________________________________________________________________________
1530Bool_t AliReconstruction::CreateVertexer()
1531{
1532// create the vertexer
1533
b8cd5251 1534 fVertexer = NULL;
1535 AliReconstructor* itsReconstructor = GetReconstructor(0);
59697224 1536 if (itsReconstructor) {
b8cd5251 1537 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
2257f27e 1538 }
b8cd5251 1539 if (!fVertexer) {
815c2b38 1540 AliWarning("couldn't create a vertexer for ITS");
2257f27e 1541 if (fStopOnError) return kFALSE;
1542 }
1543
1544 return kTRUE;
1545}
1546
24f7a148 1547//_____________________________________________________________________________
b8cd5251 1548Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
24f7a148 1549{
f08fc9f5 1550// create the trackers
24f7a148 1551
b8cd5251 1552 TString detStr = detectors;
1553 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1554 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1555 AliReconstructor* reconstructor = GetReconstructor(iDet);
1556 if (!reconstructor) continue;
1557 TString detName = fgkDetectorName[iDet];
1f46a9ae 1558 if (detName == "HLT") {
1559 fRunHLTTracking = kTRUE;
1560 continue;
1561 }
f08fc9f5 1562
1563 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1564 if (!fTracker[iDet] && (iDet < 7)) {
1565 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
8250d5f5 1566 if (fStopOnError) return kFALSE;
1567 }
1568 }
1569
24f7a148 1570 return kTRUE;
1571}
1572
e583c30d 1573//_____________________________________________________________________________
b26c3770 1574void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
e583c30d 1575{
1576// delete trackers and the run loader and close and delete the file
1577
b8cd5251 1578 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1579 delete fReconstructor[iDet];
1580 fReconstructor[iDet] = NULL;
1581 fLoader[iDet] = NULL;
1582 delete fTracker[iDet];
1583 fTracker[iDet] = NULL;
1584 }
1585 delete fVertexer;
1586 fVertexer = NULL;
9178838a 1587 delete fDiamondProfile;
1588 fDiamondProfile = NULL;
e583c30d 1589
1590 delete fRunLoader;
1591 fRunLoader = NULL;
b649205a 1592 delete fRawReader;
1593 fRawReader = NULL;
e583c30d 1594
1595 if (file) {
1596 file->Close();
1597 delete file;
1598 }
b26c3770 1599
1600 if (fileOld) {
1601 fileOld->Close();
1602 delete fileOld;
1603 gSystem->Unlink("AliESDs.old.root");
1604 }
e583c30d 1605}
24f7a148 1606
1607
1608//_____________________________________________________________________________
1609Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1610{
1611// read the ESD event from a file
1612
1613 if (!esd) return kFALSE;
1614 char fileName[256];
1615 sprintf(fileName, "ESD_%d.%d_%s.root",
1616 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1617 if (gSystem->AccessPathName(fileName)) return kFALSE;
1618
f3a97c86 1619 AliInfo(Form("reading ESD from file %s", fileName));
815c2b38 1620 AliDebug(1, Form("reading ESD from file %s", fileName));
24f7a148 1621 TFile* file = TFile::Open(fileName);
1622 if (!file || !file->IsOpen()) {
815c2b38 1623 AliError(Form("opening %s failed", fileName));
24f7a148 1624 delete file;
1625 return kFALSE;
1626 }
1627
1628 gROOT->cd();
1629 delete esd;
1630 esd = (AliESD*) file->Get("ESD");
1631 file->Close();
1632 delete file;
1633 return kTRUE;
1634}
1635
1636//_____________________________________________________________________________
1637void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1638{
1639// write the ESD event to a file
1640
1641 if (!esd) return;
1642 char fileName[256];
1643 sprintf(fileName, "ESD_%d.%d_%s.root",
1644 esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1645
815c2b38 1646 AliDebug(1, Form("writing ESD to file %s", fileName));
24f7a148 1647 TFile* file = TFile::Open(fileName, "recreate");
1648 if (!file || !file->IsOpen()) {
815c2b38 1649 AliError(Form("opening %s failed", fileName));
24f7a148 1650 } else {
1651 esd->Write("ESD");
1652 file->Close();
1653 }
1654 delete file;
1655}
f3a97c86 1656
1657
1658
1659
1660//_____________________________________________________________________________
1661void AliReconstruction::CreateTag(TFile* file)
1662{
2bdb9d38 1663 /////////////
1664 //muon code//
1665 ////////////
56982dd1 1666 Double_t fMUONMASS = 0.105658369;
2bdb9d38 1667 //Variables
56982dd1 1668 Double_t fX,fY,fZ ;
1669 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1670 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1671 Int_t fCharge;
1672 TLorentzVector fEPvector;
1673
1674 Float_t fZVertexCut = 10.0;
1675 Float_t fRhoVertexCut = 2.0;
1676
1677 Float_t fLowPtCut = 1.0;
1678 Float_t fHighPtCut = 3.0;
1679 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 1680 ////////////
1681
1682 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1683
089bf903 1684 // Creates the tags for all the events in a given ESD file
4302e20f 1685 Int_t ntrack;
089bf903 1686 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1687 Int_t nPos, nNeg, nNeutr;
1688 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1689 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1690 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1691 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1692 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 1693 Int_t fVertexflag;
1387d165 1694 Int_t iRunNumber = 0;
547d75a6 1695 TString fVertexName("default");
4302e20f 1696
f3a97c86 1697 AliRunTag *tag = new AliRunTag();
f3a97c86 1698 AliEventTag *evTag = new AliEventTag();
1699 TTree ttag("T","A Tree with event tags");
2bdb9d38 1700 TBranch * btag = ttag.Branch("AliTAG", &tag);
f3a97c86 1701 btag->SetCompressionLevel(9);
2bdb9d38 1702
f3a97c86 1703 AliInfo(Form("Creating the tags......."));
1704
1705 if (!file || !file->IsOpen()) {
1706 AliError(Form("opening failed"));
1707 delete file;
1708 return ;
2bdb9d38 1709 }
d71aa170 1710 Int_t lastEvent = 0;
f3a97c86 1711 TTree *t = (TTree*) file->Get("esdTree");
1712 TBranch * b = t->GetBranch("ESD");
1713 AliESD *esd = 0;
1714 b->SetAddress(&esd);
1387d165 1715
d71aa170 1716 b->GetEntry(fFirstEvent);
1387d165 1717 Int_t iInitRunNumber = esd->GetRunNumber();
2bdb9d38 1718
bb0901a4 1719 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
d71aa170 1720 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1721 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
2bdb9d38 1722 ntrack = 0;
1723 nPos = 0;
1724 nNeg = 0;
1725 nNeutr =0;
1726 nK0s = 0;
1727 nNeutrons = 0;
1728 nPi0s = 0;
1729 nGamas = 0;
1730 nProtons = 0;
1731 nKaons = 0;
1732 nPions = 0;
1733 nMuons = 0;
1734 nElectrons = 0;
1735 nCh1GeV = 0;
1736 nCh3GeV = 0;
1737 nCh10GeV = 0;
1738 nMu1GeV = 0;
1739 nMu3GeV = 0;
1740 nMu10GeV = 0;
1741 nEl1GeV = 0;
1742 nEl3GeV = 0;
1743 nEl10GeV = 0;
1744 maxPt = .0;
1745 meanPt = .0;
1746 totalP = .0;
547d75a6 1747 fVertexflag = 0;
d1a50cb5 1748
2bdb9d38 1749 b->GetEntry(iEventNumber);
1387d165 1750 iRunNumber = esd->GetRunNumber();
1751 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 1752 const AliESDVertex * vertexIn = esd->GetVertex();
547d75a6 1753 if (!vertexIn) AliError("ESD has not defined vertex.");
1754 if (vertexIn) fVertexName = vertexIn->GetName();
1755 if(fVertexName != "default") fVertexflag = 1;
2bdb9d38 1756 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1757 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1758 UInt_t status = esdTrack->GetStatus();
f3a97c86 1759
2bdb9d38 1760 //select only tracks with ITS refit
1761 if ((status&AliESDtrack::kITSrefit)==0) continue;
1762 //select only tracks with TPC refit
1763 if ((status&AliESDtrack::kTPCrefit)==0) continue;
4302e20f 1764
2bdb9d38 1765 //select only tracks with the "combined PID"
1766 if ((status&AliESDtrack::kESDpid)==0) continue;
1767 Double_t p[3];
1768 esdTrack->GetPxPyPz(p);
1769 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1770 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1771 totalP += momentum;
1772 meanPt += fPt;
1773 if(fPt > maxPt) maxPt = fPt;
4302e20f 1774
2bdb9d38 1775 if(esdTrack->GetSign() > 0) {
1776 nPos++;
56982dd1 1777 if(fPt > fLowPtCut) nCh1GeV++;
1778 if(fPt > fHighPtCut) nCh3GeV++;
1779 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1780 }
1781 if(esdTrack->GetSign() < 0) {
1782 nNeg++;
56982dd1 1783 if(fPt > fLowPtCut) nCh1GeV++;
1784 if(fPt > fHighPtCut) nCh3GeV++;
1785 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 1786 }
1787 if(esdTrack->GetSign() == 0) nNeutr++;
4302e20f 1788
2bdb9d38 1789 //PID
1790 Double_t prob[5];
1791 esdTrack->GetESDpid(prob);
4302e20f 1792
2bdb9d38 1793 Double_t rcc = 0.0;
1794 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1795 if(rcc == 0.0) continue;
1796 //Bayes' formula
1797 Double_t w[5];
1798 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
4302e20f 1799
2bdb9d38 1800 //protons
1801 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1802 //kaons
1803 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1804 //pions
1805 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1806 //electrons
1807 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1808 nElectrons++;
56982dd1 1809 if(fPt > fLowPtCut) nEl1GeV++;
1810 if(fPt > fHighPtCut) nEl3GeV++;
1811 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 1812 }
1813 ntrack++;
1814 }//track loop
1815
1816 /////////////
1817 //muon code//
1818 ////////////
1819 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1820 // loop over all reconstructed tracks (also first track of combination)
1821 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1822 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1823 if (muonTrack == 0x0) continue;
4302e20f 1824
2bdb9d38 1825 // Coordinates at vertex
56982dd1 1826 fZ = muonTrack->GetZ();
1827 fY = muonTrack->GetBendingCoor();
1828 fX = muonTrack->GetNonBendingCoor();
4302e20f 1829
56982dd1 1830 fThetaX = muonTrack->GetThetaX();
1831 fThetaY = muonTrack->GetThetaY();
4302e20f 1832
56982dd1 1833 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1834 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1835 fPxRec = fPzRec * TMath::Tan(fThetaX);
1836 fPyRec = fPzRec * TMath::Tan(fThetaY);
1837 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
f3a97c86 1838
2bdb9d38 1839 //ChiSquare of the track if needed
56982dd1 1840 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1841 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1842 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
4302e20f 1843
2bdb9d38 1844 // total number of muons inside a vertex cut
56982dd1 1845 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 1846 nMuons++;
56982dd1 1847 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 1848 nMu1GeV++;
56982dd1 1849 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 1850 nMu3GeV++;
56982dd1 1851 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 1852 nMu10GeV++;
1853 }
1854 }
1855 }
1856 }
1857 }//muon track loop
1858
1859 // Fill the event tags
1860 if(ntrack != 0)
1861 meanPt = meanPt/ntrack;
1862
1863 evTag->SetEventId(iEventNumber+1);
547d75a6 1864 if (vertexIn) {
1865 evTag->SetVertexX(vertexIn->GetXv());
1866 evTag->SetVertexY(vertexIn->GetYv());
1867 evTag->SetVertexZ(vertexIn->GetZv());
1868 evTag->SetVertexZError(vertexIn->GetZRes());
1869 }
d1a50cb5 1870 evTag->SetVertexFlag(fVertexflag);
1871
2bdb9d38 1872 evTag->SetT0VertexZ(esd->GetT0zVertex());
1873
8bd8ac26 1874 evTag->SetTriggerMask(esd->GetTriggerMask());
1875 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 1876
32a5cab4 1877 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1878 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1879 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1880 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2bdb9d38 1881 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
1882 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1883
1884
1885 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1886 evTag->SetNumOfPosTracks(nPos);
1887 evTag->SetNumOfNegTracks(nNeg);
1888 evTag->SetNumOfNeutrTracks(nNeutr);
1889
1890 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1891 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1892 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1893 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1894
1895 evTag->SetNumOfProtons(nProtons);
1896 evTag->SetNumOfKaons(nKaons);
1897 evTag->SetNumOfPions(nPions);
1898 evTag->SetNumOfMuons(nMuons);
1899 evTag->SetNumOfElectrons(nElectrons);
1900 evTag->SetNumOfPhotons(nGamas);
1901 evTag->SetNumOfPi0s(nPi0s);
1902 evTag->SetNumOfNeutrons(nNeutrons);
1903 evTag->SetNumOfKaon0s(nK0s);
1904
1905 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1906 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1907 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1908 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1909 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1910 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1911 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1912 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1913 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1914
85c60a8e 1915 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1916 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 1917
1918 evTag->SetTotalMomentum(totalP);
1919 evTag->SetMeanPt(meanPt);
1920 evTag->SetMaxPt(maxPt);
1921
1387d165 1922 tag->SetRunId(iInitRunNumber);
2bdb9d38 1923 tag->AddEventTag(*evTag);
1924 }
bb0901a4 1925 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
d71aa170 1926 else lastEvent = fLastEvent;
f3a97c86 1927
1928 ttag.Fill();
1929 tag->Clear();
1930
1931 char fileName[256];
1932 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
d71aa170 1933 tag->GetRunId(),fFirstEvent,lastEvent );
f3a97c86 1934 AliInfo(Form("writing tags to file %s", fileName));
1935 AliDebug(1, Form("writing tags to file %s", fileName));
1936
1937 TFile* ftag = TFile::Open(fileName, "recreate");
1938 ftag->cd();
1939 ttag.Write();
1940 ftag->Close();
1941 file->cd();
1942 delete tag;
f3a97c86 1943 delete evTag;
1944}
1945
98937d93 1946void AliReconstruction::WriteAlignmentData(AliESD* esd)
1947{
1948 // Write space-points which are then used in the alignment procedures
1949 // For the moment only ITS, TRD and TPC
1950
1951 // Load TOF clusters
d528ee75 1952 if (fTracker[3]){
1953 fLoader[3]->LoadRecPoints("read");
1954 TTree* tree = fLoader[3]->TreeR();
1955 if (!tree) {
1956 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
1957 return;
1958 }
1959 fTracker[3]->LoadClusters(tree);
98937d93 1960 }
98937d93 1961 Int_t ntracks = esd->GetNumberOfTracks();
1962 for (Int_t itrack = 0; itrack < ntracks; itrack++)
1963 {
1964 AliESDtrack *track = esd->GetTrack(itrack);
1965 Int_t nsp = 0;
ef7253ac 1966 Int_t idx[200];
98937d93 1967 for (Int_t iDet = 3; iDet >= 0; iDet--)
1968 nsp += track->GetNcls(iDet);
1969 if (nsp) {
1970 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
1971 track->SetTrackPointArray(sp);
1972 Int_t isptrack = 0;
1973 for (Int_t iDet = 3; iDet >= 0; iDet--) {
1974 AliTracker *tracker = fTracker[iDet];
1975 if (!tracker) continue;
1976 Int_t nspdet = track->GetNcls(iDet);
98937d93 1977 if (nspdet <= 0) continue;
1978 track->GetClusters(iDet,idx);
1979 AliTrackPoint p;
1980 Int_t isp = 0;
1981 Int_t isp2 = 0;
1982 while (isp < nspdet) {
1983 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
160db090 1984 const Int_t kNTPCmax = 159;
1985 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
98937d93 1986 if (!isvalid) continue;
1987 sp->AddPoint(isptrack,&p); isptrack++; isp++;
1988 }
98937d93 1989 }
1990 }
1991 }
d528ee75 1992 if (fTracker[3]){
1993 fTracker[3]->UnloadClusters();
1994 fLoader[3]->UnloadRecPoints();
1995 }
98937d93 1996}