Moving the alignment-related static methods from AliAlignObj to the new geometry...
[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// 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// //
59// The name of the galice file can be changed from the default //
60// "galice.root" by passing it as argument to the AliReconstruction //
61// constructor or by //
62// //
63// rec.SetGAliceFile("..."); //
64// //
65// The local reconstruction can be switched on or off for individual //
66// detectors by //
67// //
68// rec.SetRunLocalReconstruction("..."); //
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// //
74// The reconstruction of the primary vertex position can be switched off by //
75// //
76// rec.SetRunVertexFinder(kFALSE); //
77// //
78// The tracking and the creation of ESD tracks can be switched on for //
79// selected detectors by //
80// //
81// rec.SetRunTracking("..."); //
82// //
83// Uniform/nonuniform field tracking switches (default: uniform field) //
84// //
85// rec.SetUniformFieldTracking(); ( rec.SetUniformFieldTracking(kFALSE); ) //
86// //
87// The filling of additional ESD information can be steered by //
88// //
89// rec.SetFillESD("..."); //
90// //
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. //
100// //
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. //
103// //
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// //
115///////////////////////////////////////////////////////////////////////////////
116
117#include <TArrayF.h>
118#include <TFile.h>
119#include <TSystem.h>
120#include <TROOT.h>
121#include <TPluginManager.h>
122#include <TStopwatch.h>
123#include <TGeoManager.h>
124#include <TLorentzVector.h>
125
126#include "AliReconstruction.h"
127#include "AliReconstructor.h"
128#include "AliLog.h"
129#include "AliRunLoader.h"
130#include "AliRun.h"
131#include "AliRawReaderFile.h"
132#include "AliRawReaderDate.h"
133#include "AliRawReaderRoot.h"
134#include "AliRawEventHeaderBase.h"
135#include "AliESD.h"
136#include "AliESDfriend.h"
137#include "AliESDVertex.h"
138#include "AliMultiplicity.h"
139#include "AliTracker.h"
140#include "AliVertexer.h"
141#include "AliVertexerTracks.h"
142#include "AliV0vertexer.h"
143#include "AliCascadeVertexer.h"
144#include "AliHeader.h"
145#include "AliGenEventHeader.h"
146#include "AliPID.h"
147#include "AliESDpid.h"
148#include "AliESDtrack.h"
149
150#include "AliRunTag.h"
151#include "AliDetectorTag.h"
152#include "AliEventTag.h"
153
154#include "AliTrackPointArray.h"
155#include "AliCDBManager.h"
156#include "AliCDBEntry.h"
157#include "AliAlignObj.h"
158
159#include "AliCentralTrigger.h"
160#include "AliCTPRawStream.h"
161
162#include "AliAODEvent.h"
163#include "AliAODHeader.h"
164#include "AliAODTrack.h"
165#include "AliAODVertex.h"
166#include "AliAODCluster.h"
167
168
169ClassImp(AliReconstruction)
170
171
172//_____________________________________________________________________________
173const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "HMPID", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "T0", "VZERO", "ACORDE", "HLT"};
174
175//_____________________________________________________________________________
176AliReconstruction::AliReconstruction(const char* gAliceFilename, const char* cdbUri,
177 const char* name, const char* title) :
178 TNamed(name, title),
179
180 fUniformField(kTRUE),
181 fRunVertexFinder(kTRUE),
182 fRunHLTTracking(kFALSE),
183 fRunMuonTracking(kFALSE),
184 fStopOnError(kFALSE),
185 fWriteAlignmentData(kFALSE),
186 fWriteESDfriend(kFALSE),
187 fWriteAOD(kFALSE),
188 fFillTriggerESD(kTRUE),
189
190 fRunLocalReconstruction("ALL"),
191 fRunTracking("ALL"),
192 fFillESD("ALL"),
193 fGAliceFileName(gAliceFilename),
194 fInput(""),
195 fEquipIdMap(""),
196 fFirstEvent(0),
197 fLastEvent(-1),
198 fNumberOfEventsPerFile(1),
199 fCheckPointLevel(0),
200 fOptions(),
201 fLoadAlignFromCDB(kTRUE),
202 fLoadAlignData("ALL"),
203
204 fRunLoader(NULL),
205 fRawReader(NULL),
206
207 fVertexer(NULL),
208 fDiamondProfile(NULL),
209
210 fAlignObjArray(NULL),
211 fCDBUri(cdbUri),
212 fSpecCDBUri()
213{
214// create reconstruction object with default parameters
215
216 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
217 fReconstructor[iDet] = NULL;
218 fLoader[iDet] = NULL;
219 fTracker[iDet] = NULL;
220 }
221 AliPID pid;
222}
223
224//_____________________________________________________________________________
225AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
226 TNamed(rec),
227
228 fUniformField(rec.fUniformField),
229 fRunVertexFinder(rec.fRunVertexFinder),
230 fRunHLTTracking(rec.fRunHLTTracking),
231 fRunMuonTracking(rec.fRunMuonTracking),
232 fStopOnError(rec.fStopOnError),
233 fWriteAlignmentData(rec.fWriteAlignmentData),
234 fWriteESDfriend(rec.fWriteESDfriend),
235 fWriteAOD(rec.fWriteAOD),
236 fFillTriggerESD(rec.fFillTriggerESD),
237
238 fRunLocalReconstruction(rec.fRunLocalReconstruction),
239 fRunTracking(rec.fRunTracking),
240 fFillESD(rec.fFillESD),
241 fGAliceFileName(rec.fGAliceFileName),
242 fInput(rec.fInput),
243 fEquipIdMap(rec.fEquipIdMap),
244 fFirstEvent(rec.fFirstEvent),
245 fLastEvent(rec.fLastEvent),
246 fNumberOfEventsPerFile(rec.fNumberOfEventsPerFile),
247 fCheckPointLevel(0),
248 fOptions(),
249 fLoadAlignFromCDB(rec.fLoadAlignFromCDB),
250 fLoadAlignData(rec.fLoadAlignData),
251
252 fRunLoader(NULL),
253 fRawReader(NULL),
254
255 fVertexer(NULL),
256 fDiamondProfile(NULL),
257
258 fAlignObjArray(rec.fAlignObjArray),
259 fCDBUri(rec.fCDBUri),
260 fSpecCDBUri()
261{
262// copy constructor
263
264 for (Int_t i = 0; i < rec.fOptions.GetEntriesFast(); i++) {
265 if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
266 }
267 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
268 fReconstructor[iDet] = NULL;
269 fLoader[iDet] = NULL;
270 fTracker[iDet] = NULL;
271 }
272 for (Int_t i = 0; i < rec.fSpecCDBUri.GetEntriesFast(); i++) {
273 if (rec.fSpecCDBUri[i]) fSpecCDBUri.Add(rec.fSpecCDBUri[i]->Clone());
274 }
275}
276
277//_____________________________________________________________________________
278AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
279{
280// assignment operator
281
282 this->~AliReconstruction();
283 new(this) AliReconstruction(rec);
284 return *this;
285}
286
287//_____________________________________________________________________________
288AliReconstruction::~AliReconstruction()
289{
290// clean up
291
292 CleanUp();
293 fOptions.Delete();
294 fSpecCDBUri.Delete();
295}
296
297//_____________________________________________________________________________
298void AliReconstruction::InitCDBStorage()
299{
300// activate a default CDB storage
301// First check if we have any CDB storage set, because it is used
302// to retrieve the calibration and alignment constants
303
304 AliCDBManager* man = AliCDBManager::Instance();
305 if (man->IsDefaultStorageSet())
306 {
307 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
308 AliWarning("Default CDB storage has been already set !");
309 AliWarning(Form("Ignoring the default storage declared in AliReconstruction: %s",fCDBUri.Data()));
310 AliWarning("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
311 fCDBUri = "";
312 }
313 else {
314 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
315 AliDebug(2, Form("Default CDB storage is set to: %s",fCDBUri.Data()));
316 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
317 man->SetDefaultStorage(fCDBUri);
318 }
319
320 // Now activate the detector specific CDB storage locations
321 for (Int_t i = 0; i < fSpecCDBUri.GetEntriesFast(); i++) {
322 TObject* obj = fSpecCDBUri[i];
323 if (!obj) continue;
324 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
325 AliDebug(2, Form("Specific CDB storage for %s is set to: %s",obj->GetName(),obj->GetTitle()));
326 AliDebug(2, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
327 man->SetSpecificStorage(obj->GetName(), obj->GetTitle());
328 }
329 man->Print();
330}
331
332//_____________________________________________________________________________
333void AliReconstruction::SetDefaultStorage(const char* uri) {
334// Store the desired default CDB storage location
335// Activate it later within the Run() method
336
337 fCDBUri = uri;
338
339}
340
341//_____________________________________________________________________________
342void AliReconstruction::SetSpecificStorage(const char* calibType, const char* uri) {
343// Store a detector-specific CDB storage location
344// Activate it later within the Run() method
345
346 AliCDBPath aPath(calibType);
347 if(!aPath.IsValid()){
348 // if calibType is not wildcard but it is a valid detector, add "/*" to make it a valid path
349 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
350 if(!strcmp(calibType, fgkDetectorName[iDet])) {
351 aPath.SetPath(Form("%s/*", calibType));
352 AliInfo(Form("Path for specific storage set to %s", aPath.GetPath().Data()));
353 break;
354 }
355 }
356 if(!aPath.IsValid()){
357 AliError(Form("Not a valid path or detector: %s", calibType));
358 return;
359 }
360 }
361
362 // check that calibType refers to a "valid" detector name
363 Bool_t isDetector = kFALSE;
364 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
365 TString detName = fgkDetectorName[iDet];
366 if(aPath.GetLevel0() == detName) {
367 isDetector = kTRUE;
368 break;
369 }
370 }
371
372 if(!isDetector) {
373 AliError(Form("Not a valid detector: %s", aPath.GetLevel0().Data()));
374 return;
375 }
376
377 TObject* obj = fSpecCDBUri.FindObject(aPath.GetPath().Data());
378 if (obj) fSpecCDBUri.Remove(obj);
379 fSpecCDBUri.Add(new TNamed(aPath.GetPath().Data(), uri));
380
381}
382
383//_____________________________________________________________________________
384Bool_t AliReconstruction::SetRunNumber()
385{
386 // The method is called in Run() in order
387 // to set a correct run number.
388 // In case of raw data reconstruction the
389 // run number is taken from the raw data header
390
391 if(AliCDBManager::Instance()->GetRun() < 0) {
392 if (!fRunLoader) {
393 AliError("No run loader is found !");
394 return kFALSE;
395 }
396 // read run number from gAlice
397 if(fRunLoader->GetAliRun())
398 AliCDBManager::Instance()->SetRun(fRunLoader->GetAliRun()->GetRunNumber());
399 else {
400 if(fRawReader) {
401 if(fRawReader->NextEvent()) {
402 AliCDBManager::Instance()->SetRun(fRawReader->GetRunNumber());
403 fRawReader->RewindEvents();
404 }
405 else {
406 AliError("No raw-data events found !");
407 return kFALSE;
408 }
409 }
410 else {
411 AliError("Neither gAlice nor RawReader objects are found !");
412 return kFALSE;
413 }
414 }
415 AliInfo(Form("CDB Run number: %d",AliCDBManager::Instance()->GetRun()));
416 }
417 return kTRUE;
418}
419
420//_____________________________________________________________________________
421Bool_t AliReconstruction::ApplyAlignObjsToGeom(TObjArray* alObjArray)
422{
423 // Read collection of alignment objects (AliAlignObj derived) saved
424 // in the TClonesArray ClArrayName and apply them to the geometry
425 // manager singleton.
426 //
427 alObjArray->Sort();
428 Int_t nvols = alObjArray->GetEntriesFast();
429
430 Bool_t flag = kTRUE;
431
432 for(Int_t j=0; j<nvols; j++)
433 {
434 AliAlignObj* alobj = (AliAlignObj*) alObjArray->UncheckedAt(j);
435 if (alobj->ApplyToGeometry() == kFALSE) flag = kFALSE;
436 }
437
438 if (AliDebugLevelClass() >= 1) {
439 gGeoManager->GetTopNode()->CheckOverlaps(20);
440 TObjArray* ovexlist = gGeoManager->GetListOfOverlaps();
441 if(ovexlist->GetEntriesFast()){
442 AliError("The application of alignment objects to the geometry caused huge overlaps/extrusions!");
443 }
444 }
445
446 return flag;
447
448}
449
450//_____________________________________________________________________________
451Bool_t AliReconstruction::SetAlignObjArraySingleDet(const char* detName)
452{
453 // Fills array of single detector's alignable objects from CDB
454
455 AliDebug(2, Form("Loading alignment data for detector: %s",detName));
456
457 AliCDBEntry *entry;
458
459 AliCDBPath path(detName,"Align","Data");
460
461 entry=AliCDBManager::Instance()->Get(path.GetPath());
462 if(!entry){
463 AliDebug(2,Form("Couldn't load alignment data for detector %s",detName));
464 return kFALSE;
465 }
466 entry->SetOwner(1);
467 TClonesArray *alignArray = (TClonesArray*) entry->GetObject();
468 alignArray->SetOwner(0);
469 AliDebug(2,Form("Found %d alignment objects for %s",
470 alignArray->GetEntries(),detName));
471
472 AliAlignObj *alignObj=0;
473 TIter iter(alignArray);
474
475 // loop over align objects in detector
476 while( ( alignObj=(AliAlignObj *) iter.Next() ) ){
477 fAlignObjArray->Add(alignObj);
478 }
479 // delete entry --- Don't delete, it is cached!
480
481 AliDebug(2, Form("fAlignObjArray entries: %d",fAlignObjArray->GetEntries() ));
482 return kTRUE;
483
484}
485
486//_____________________________________________________________________________
487Bool_t AliReconstruction::MisalignGeometry(const TString& detectors)
488{
489 // Read the alignment objects from CDB.
490 // Each detector is supposed to have the
491 // alignment objects in DET/Align/Data CDB path.
492 // All the detector objects are then collected,
493 // sorted by geometry level (starting from ALIC) and
494 // then applied to the TGeo geometry.
495 // Finally an overlaps check is performed.
496
497 // Load alignment data from CDB and fill fAlignObjArray
498 if(fLoadAlignFromCDB){
499 if(!fAlignObjArray) fAlignObjArray = new TObjArray();
500
501 //fAlignObjArray->RemoveAll();
502 fAlignObjArray->Clear();
503 fAlignObjArray->SetOwner(0);
504
505 TString detStr = detectors;
506 TString dataNotLoaded="";
507 TString dataLoaded="";
508
509 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
510 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
511 if(!SetAlignObjArraySingleDet(fgkDetectorName[iDet])){
512 dataNotLoaded += fgkDetectorName[iDet];
513 dataNotLoaded += " ";
514 } else {
515 dataLoaded += fgkDetectorName[iDet];
516 dataLoaded += " ";
517 }
518 } // end loop over detectors
519
520 if ((detStr.CompareTo("ALL") == 0)) detStr = "";
521 dataNotLoaded += detStr;
522 if(!dataLoaded.IsNull()) AliInfo(Form("Alignment data loaded for: %s",
523 dataLoaded.Data()));
524 if(!dataNotLoaded.IsNull()) AliInfo(Form("Didn't/couldn't load alignment data for: %s",
525 dataNotLoaded.Data()));
526 } // fLoadAlignFromCDB flag
527
528 // Check if the array with alignment objects was
529 // provided by the user. If yes, apply the objects
530 // to the present TGeo geometry
531 if (fAlignObjArray) {
532 if (gGeoManager && gGeoManager->IsClosed()) {
533 if (ApplyAlignObjsToGeom(fAlignObjArray) == kFALSE) {
534 AliError("The misalignment of one or more volumes failed!"
535 "Compare the list of simulated detectors and the list of detector alignment data!");
536 return kFALSE;
537 }
538 }
539 else {
540 AliError("Can't apply the misalignment! gGeoManager doesn't exist or it is still opened!");
541 return kFALSE;
542 }
543 }
544
545 delete fAlignObjArray; fAlignObjArray=0;
546
547 // Update the TGeoPhysicalNodes
548 gGeoManager->RefreshPhysicalNodes();
549
550 return kTRUE;
551}
552
553//_____________________________________________________________________________
554void AliReconstruction::SetGAliceFile(const char* fileName)
555{
556// set the name of the galice file
557
558 fGAliceFileName = fileName;
559}
560
561//_____________________________________________________________________________
562void AliReconstruction::SetOption(const char* detector, const char* option)
563{
564// set options for the reconstruction of a detector
565
566 TObject* obj = fOptions.FindObject(detector);
567 if (obj) fOptions.Remove(obj);
568 fOptions.Add(new TNamed(detector, option));
569}
570
571
572//_____________________________________________________________________________
573Bool_t AliReconstruction::Run(const char* input)
574{
575// run the reconstruction
576
577 // set the input
578 if (!input) input = fInput.Data();
579 TString fileName(input);
580 if (fileName.EndsWith("/")) {
581 fRawReader = new AliRawReaderFile(fileName);
582 } else if (fileName.EndsWith(".root")) {
583 fRawReader = new AliRawReaderRoot(fileName);
584 } else if (!fileName.IsNull()) {
585 fRawReader = new AliRawReaderDate(fileName);
586 fRawReader->SelectEvents(7);
587 }
588 if (!fEquipIdMap.IsNull() && fRawReader)
589 fRawReader->LoadEquipmentIdsMap(fEquipIdMap);
590
591
592 // get the run loader
593 if (!InitRunLoader()) return kFALSE;
594
595 // Initialize the CDB storage
596 InitCDBStorage();
597
598 // Set run number in CDBManager (if it is not already set by the user)
599 if (!SetRunNumber()) if (fStopOnError) return kFALSE;
600
601 // Import ideal TGeo geometry and apply misalignment
602 if (!gGeoManager) {
603 TString geom(gSystem->DirName(fGAliceFileName));
604 geom += "/geometry.root";
605 TGeoManager::Import(geom.Data());
606 if (!gGeoManager) if (fStopOnError) return kFALSE;
607 }
608
609 AliCDBManager* man = AliCDBManager::Instance();
610 if (!MisalignGeometry(fLoadAlignData)) if (fStopOnError) return kFALSE;
611
612 // local reconstruction
613 if (!fRunLocalReconstruction.IsNull()) {
614 if (!RunLocalReconstruction(fRunLocalReconstruction)) {
615 if (fStopOnError) {CleanUp(); return kFALSE;}
616 }
617 }
618// if (!fRunVertexFinder && fRunTracking.IsNull() &&
619// fFillESD.IsNull()) return kTRUE;
620
621 // get vertexer
622 if (fRunVertexFinder && !CreateVertexer()) {
623 if (fStopOnError) {
624 CleanUp();
625 return kFALSE;
626 }
627 }
628
629 // get trackers
630 if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
631 if (fStopOnError) {
632 CleanUp();
633 return kFALSE;
634 }
635 }
636
637
638 TStopwatch stopwatch;
639 stopwatch.Start();
640
641 // get the possibly already existing ESD file and tree
642 AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
643 TFile* fileOld = NULL;
644 TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
645 if (!gSystem->AccessPathName("AliESDs.root")){
646 gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
647 fileOld = TFile::Open("AliESDs.old.root");
648 if (fileOld && fileOld->IsOpen()) {
649 treeOld = (TTree*) fileOld->Get("esdTree");
650 if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
651 hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
652 if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
653 }
654 }
655
656 // create the ESD output file and tree
657 TFile* file = TFile::Open("AliESDs.root", "RECREATE");
658 if (!file->IsOpen()) {
659 AliError("opening AliESDs.root failed");
660 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
661 }
662 TTree* tree = new TTree("esdTree", "Tree with ESD objects");
663 tree->Branch("ESD", "AliESD", &esd);
664 TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
665 hlttree->Branch("ESD", "AliESD", &hltesd);
666 delete esd; delete hltesd;
667 esd = NULL; hltesd = NULL;
668
669 // create the branch with ESD additions
670 AliESDfriend *esdf=0;
671 if (fWriteESDfriend) {
672 TBranch *br=tree->Branch("ESDfriend.", "AliESDfriend", &esdf);
673 br->SetFile("AliESDfriends.root");
674 }
675
676 // Get the diamond profile from OCDB
677 AliCDBEntry* entry = AliCDBManager::Instance()
678 ->Get("GRP/Calib/MeanVertex");
679
680 if(entry) {
681 fDiamondProfile = dynamic_cast<AliESDVertex*> (entry->GetObject());
682 } else {
683 AliError("No diamond profile found in OCDB!");
684 }
685
686 AliVertexerTracks tVertexer(AliTracker::GetBz());
687 if(fDiamondProfile) tVertexer.SetVtxStart(fDiamondProfile);
688
689 // loop over events
690 if (fRawReader) fRawReader->RewindEvents();
691
692 for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
693 if (fRawReader) fRawReader->NextEvent();
694 if ((iEvent < fFirstEvent) || ((fLastEvent >= 0) && (iEvent > fLastEvent))) {
695 // copy old ESD to the new one
696 if (treeOld) {
697 treeOld->SetBranchAddress("ESD", &esd);
698 treeOld->GetEntry(iEvent);
699 }
700 tree->Fill();
701 if (hlttreeOld) {
702 hlttreeOld->SetBranchAddress("ESD", &hltesd);
703 hlttreeOld->GetEntry(iEvent);
704 }
705 hlttree->Fill();
706 continue;
707 }
708
709 AliInfo(Form("processing event %d", iEvent));
710 fRunLoader->GetEvent(iEvent);
711
712 char aFileName[256];
713 sprintf(aFileName, "ESD_%d.%d_final.root",
714 fRunLoader->GetHeader()->GetRun(),
715 fRunLoader->GetHeader()->GetEventNrInRun());
716 if (!gSystem->AccessPathName(aFileName)) continue;
717
718 // local reconstruction
719 if (!fRunLocalReconstruction.IsNull()) {
720 if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
721 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
722 }
723 }
724
725 esd = new AliESD; hltesd = new AliESD;
726 esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
727 hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
728 esd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
729 hltesd->SetEventNumberInFile(fRunLoader->GetHeader()->GetEventNrInRun());
730
731 // Set magnetic field from the tracker
732 esd->SetMagneticField(AliTracker::GetBz());
733 hltesd->SetMagneticField(AliTracker::GetBz());
734
735 // Fill raw-data error log into the ESD
736 if (fRawReader) FillRawDataErrorLog(iEvent,esd);
737
738 // vertex finder
739 if (fRunVertexFinder) {
740 if (!ReadESD(esd, "vertex")) {
741 if (!RunVertexFinder(esd)) {
742 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
743 }
744 if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
745 }
746 }
747
748 // HLT tracking
749 if (!fRunTracking.IsNull()) {
750 if (fRunHLTTracking) {
751 hltesd->SetVertex(esd->GetVertex());
752 if (!RunHLTTracking(hltesd)) {
753 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
754 }
755 }
756 }
757
758 // Muon tracking
759 if (!fRunTracking.IsNull()) {
760 if (fRunMuonTracking) {
761 if (!RunMuonTracking()) {
762 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
763 }
764 }
765 }
766
767 // barrel tracking
768 if (!fRunTracking.IsNull()) {
769 if (!ReadESD(esd, "tracking")) {
770 if (!RunTracking(esd)) {
771 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
772 }
773 if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
774 }
775 }
776
777 // fill ESD
778 if (!fFillESD.IsNull()) {
779 if (!FillESD(esd, fFillESD)) {
780 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
781 }
782 }
783 // fill Event header information from the RawEventHeader
784 if (fRawReader){FillRawEventHeaderESD(esd);}
785
786 // combined PID
787 AliESDpid::MakePID(esd);
788 if (fCheckPointLevel > 1) WriteESD(esd, "PID");
789
790 if (fFillTriggerESD) {
791 if (!ReadESD(esd, "trigger")) {
792 if (!FillTriggerESD(esd)) {
793 if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
794 }
795 if (fCheckPointLevel > 1) WriteESD(esd, "trigger");
796 }
797 }
798
799 //Try to improve the reconstructed primary vertex position using the tracks
800 AliESDVertex *pvtx=tVertexer.FindPrimaryVertex(esd);
801 if(fDiamondProfile) esd->SetDiamond(fDiamondProfile);
802
803 if (pvtx)
804 if (pvtx->GetStatus()) {
805 // Store the improved primary vertex
806 esd->SetPrimaryVertex(pvtx);
807 // Propagate the tracks to the DCA to the improved primary vertex
808 Double_t somethingbig = 777.;
809 Double_t bz = esd->GetMagneticField();
810 Int_t nt=esd->GetNumberOfTracks();
811 while (nt--) {
812 AliESDtrack *t = esd->GetTrack(nt);
813 t->RelateToVertex(pvtx, bz, somethingbig);
814 }
815 }
816
817 {
818 // V0 finding
819 AliV0vertexer vtxer;
820 vtxer.Tracks2V0vertices(esd);
821
822 // Cascade finding
823 AliCascadeVertexer cvtxer;
824 cvtxer.V0sTracks2CascadeVertices(esd);
825 }
826
827 // write ESD
828 if (fWriteESDfriend) {
829 esdf=new AliESDfriend();
830 esd->GetESDfriend(esdf);
831 }
832 tree->Fill();
833
834 // write HLT ESD
835 hlttree->Fill();
836
837 if (fCheckPointLevel > 0) WriteESD(esd, "final");
838
839 delete esd; delete esdf; delete hltesd;
840 esd = NULL; esdf=NULL; hltesd = NULL;
841 }
842
843 AliInfo(Form("Execution time for filling ESD : R:%.2fs C:%.2fs",
844 stopwatch.RealTime(),stopwatch.CpuTime()));
845
846 file->cd();
847 if (fWriteESDfriend)
848 tree->SetBranchStatus("ESDfriend*",0);
849 tree->Write();
850 hlttree->Write();
851
852 if (fWriteAOD) {
853 TFile *aodFile = TFile::Open("AliAOD.root", "RECREATE");
854 ESDFile2AODFile(file, aodFile);
855 aodFile->Close();
856 }
857
858 // Create tags for the events in the ESD tree (the ESD tree is always present)
859 // In case of empty events the tags will contain dummy values
860 CreateTag(file);
861 CleanUp(file, fileOld);
862
863 return kTRUE;
864}
865
866
867//_____________________________________________________________________________
868Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
869{
870// run the local reconstruction
871
872 TStopwatch stopwatch;
873 stopwatch.Start();
874
875 AliCDBManager* man = AliCDBManager::Instance();
876 Bool_t origCache = man->GetCacheFlag();
877
878 TString detStr = detectors;
879 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
880 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
881 AliReconstructor* reconstructor = GetReconstructor(iDet);
882 if (!reconstructor) continue;
883 if (reconstructor->HasLocalReconstruction()) continue;
884
885 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
886 TStopwatch stopwatchDet;
887 stopwatchDet.Start();
888
889 AliInfo(Form("Loading calibration data from OCDB for %s", fgkDetectorName[iDet]));
890
891 man->SetCacheFlag(kTRUE);
892 TString calibPath = Form("%s/Calib/*", fgkDetectorName[iDet]);
893 man->GetAll(calibPath); // entries are cached!
894
895 if (fRawReader) {
896 fRawReader->RewindEvents();
897 reconstructor->Reconstruct(fRunLoader, fRawReader);
898 } else {
899 reconstructor->Reconstruct(fRunLoader);
900 }
901 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
902 fgkDetectorName[iDet],
903 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
904
905 // unload calibration data
906 man->ClearCache();
907 }
908
909 man->SetCacheFlag(origCache);
910
911 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
912 AliError(Form("the following detectors were not found: %s",
913 detStr.Data()));
914 if (fStopOnError) return kFALSE;
915 }
916
917 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
918 stopwatch.RealTime(),stopwatch.CpuTime()));
919
920 return kTRUE;
921}
922
923//_____________________________________________________________________________
924Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
925{
926// run the local reconstruction
927
928 TStopwatch stopwatch;
929 stopwatch.Start();
930
931 TString detStr = detectors;
932 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
933 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
934 AliReconstructor* reconstructor = GetReconstructor(iDet);
935 if (!reconstructor) continue;
936 AliLoader* loader = fLoader[iDet];
937
938 // conversion of digits
939 if (fRawReader && reconstructor->HasDigitConversion()) {
940 AliInfo(Form("converting raw data digits into root objects for %s",
941 fgkDetectorName[iDet]));
942 TStopwatch stopwatchDet;
943 stopwatchDet.Start();
944 loader->LoadDigits("update");
945 loader->CleanDigits();
946 loader->MakeDigitsContainer();
947 TTree* digitsTree = loader->TreeD();
948 reconstructor->ConvertDigits(fRawReader, digitsTree);
949 loader->WriteDigits("OVERWRITE");
950 loader->UnloadDigits();
951 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
952 fgkDetectorName[iDet],
953 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
954 }
955
956 // local reconstruction
957 if (!reconstructor->HasLocalReconstruction()) continue;
958 AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
959 TStopwatch stopwatchDet;
960 stopwatchDet.Start();
961 loader->LoadRecPoints("update");
962 loader->CleanRecPoints();
963 loader->MakeRecPointsContainer();
964 TTree* clustersTree = loader->TreeR();
965 if (fRawReader && !reconstructor->HasDigitConversion()) {
966 reconstructor->Reconstruct(fRawReader, clustersTree);
967 } else {
968 loader->LoadDigits("read");
969 TTree* digitsTree = loader->TreeD();
970 if (!digitsTree) {
971 AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
972 if (fStopOnError) return kFALSE;
973 } else {
974 reconstructor->Reconstruct(digitsTree, clustersTree);
975 }
976 loader->UnloadDigits();
977 }
978 loader->WriteRecPoints("OVERWRITE");
979 loader->UnloadRecPoints();
980 AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
981 fgkDetectorName[iDet],
982 stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
983 }
984
985 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
986 AliError(Form("the following detectors were not found: %s",
987 detStr.Data()));
988 if (fStopOnError) return kFALSE;
989 }
990
991 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
992 stopwatch.RealTime(),stopwatch.CpuTime()));
993
994 return kTRUE;
995}
996
997//_____________________________________________________________________________
998Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
999{
1000// run the barrel tracking
1001
1002 TStopwatch stopwatch;
1003 stopwatch.Start();
1004
1005 AliESDVertex* vertex = NULL;
1006 Double_t vtxPos[3] = {0, 0, 0};
1007 Double_t vtxErr[3] = {0.07, 0.07, 0.1};
1008 TArrayF mcVertex(3);
1009 if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
1010 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
1011 for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
1012 }
1013
1014 if (fVertexer) {
1015 if(fDiamondProfile) fVertexer->SetVtxStart(fDiamondProfile);
1016 AliInfo("running the ITS vertex finder");
1017 if (fLoader[0]) fLoader[0]->LoadRecPoints();
1018 vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
1019 if (fLoader[0]) fLoader[0]->UnloadRecPoints();
1020 if(!vertex){
1021 AliWarning("Vertex not found");
1022 vertex = new AliESDVertex();
1023 vertex->SetName("default");
1024 }
1025 else {
1026 vertex->SetTruePos(vtxPos); // store also the vertex from MC
1027 vertex->SetName("reconstructed");
1028 }
1029
1030 } else {
1031 AliInfo("getting the primary vertex from MC");
1032 vertex = new AliESDVertex(vtxPos, vtxErr);
1033 }
1034
1035 if (vertex) {
1036 vertex->GetXYZ(vtxPos);
1037 vertex->GetSigmaXYZ(vtxErr);
1038 } else {
1039 AliWarning("no vertex reconstructed");
1040 vertex = new AliESDVertex(vtxPos, vtxErr);
1041 }
1042 esd->SetVertex(vertex);
1043 // if SPD multiplicity has been determined, it is stored in the ESD
1044 if (fVertexer) {
1045 AliMultiplicity *mult= fVertexer->GetMultiplicity();
1046 if(mult)esd->SetMultiplicity(mult);
1047 }
1048
1049 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1050 if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
1051 }
1052 delete vertex;
1053
1054 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1055 stopwatch.RealTime(),stopwatch.CpuTime()));
1056
1057 return kTRUE;
1058}
1059
1060//_____________________________________________________________________________
1061Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
1062{
1063// run the HLT barrel tracking
1064
1065 TStopwatch stopwatch;
1066 stopwatch.Start();
1067
1068 if (!fRunLoader) {
1069 AliError("Missing runLoader!");
1070 return kFALSE;
1071 }
1072
1073 AliInfo("running HLT tracking");
1074
1075 // Get a pointer to the HLT reconstructor
1076 AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
1077 if (!reconstructor) return kFALSE;
1078
1079 // TPC + ITS
1080 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1081 TString detName = fgkDetectorName[iDet];
1082 AliDebug(1, Form("%s HLT tracking", detName.Data()));
1083 reconstructor->SetOption(detName.Data());
1084 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1085 if (!tracker) {
1086 AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
1087 if (fStopOnError) return kFALSE;
1088 continue;
1089 }
1090 Double_t vtxPos[3];
1091 Double_t vtxErr[3]={0.005,0.005,0.010};
1092 const AliESDVertex *vertex = esd->GetVertex();
1093 vertex->GetXYZ(vtxPos);
1094 tracker->SetVertex(vtxPos,vtxErr);
1095 if(iDet != 1) {
1096 fLoader[iDet]->LoadRecPoints("read");
1097 TTree* tree = fLoader[iDet]->TreeR();
1098 if (!tree) {
1099 AliError(Form("Can't get the %s cluster tree", detName.Data()));
1100 return kFALSE;
1101 }
1102 tracker->LoadClusters(tree);
1103 }
1104 if (tracker->Clusters2Tracks(esd) != 0) {
1105 AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
1106 return kFALSE;
1107 }
1108 if(iDet != 1) {
1109 tracker->UnloadClusters();
1110 }
1111 delete tracker;
1112 }
1113
1114 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1115 stopwatch.RealTime(),stopwatch.CpuTime()));
1116
1117 return kTRUE;
1118}
1119
1120//_____________________________________________________________________________
1121Bool_t AliReconstruction::RunMuonTracking()
1122{
1123// run the muon spectrometer tracking
1124
1125 TStopwatch stopwatch;
1126 stopwatch.Start();
1127
1128 if (!fRunLoader) {
1129 AliError("Missing runLoader!");
1130 return kFALSE;
1131 }
1132 Int_t iDet = 7; // for MUON
1133
1134 AliInfo("is running...");
1135
1136 // Get a pointer to the MUON reconstructor
1137 AliReconstructor *reconstructor = GetReconstructor(iDet);
1138 if (!reconstructor) return kFALSE;
1139
1140
1141 TString detName = fgkDetectorName[iDet];
1142 AliDebug(1, Form("%s tracking", detName.Data()));
1143 AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
1144 if (!tracker) {
1145 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1146 return kFALSE;
1147 }
1148
1149 // create Tracks
1150 fLoader[iDet]->LoadTracks("update");
1151 fLoader[iDet]->CleanTracks();
1152 fLoader[iDet]->MakeTracksContainer();
1153
1154 // read RecPoints
1155 fLoader[iDet]->LoadRecPoints("read");
1156
1157 if (!tracker->Clusters2Tracks(0x0)) {
1158 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1159 return kFALSE;
1160 }
1161 fLoader[iDet]->UnloadRecPoints();
1162
1163 fLoader[iDet]->WriteTracks("OVERWRITE");
1164 fLoader[iDet]->UnloadTracks();
1165
1166 delete tracker;
1167
1168
1169 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1170 stopwatch.RealTime(),stopwatch.CpuTime()));
1171
1172 return kTRUE;
1173}
1174
1175
1176//_____________________________________________________________________________
1177Bool_t AliReconstruction::RunTracking(AliESD*& esd)
1178{
1179// run the barrel tracking
1180
1181 TStopwatch stopwatch;
1182 stopwatch.Start();
1183
1184 AliInfo("running tracking");
1185
1186 //Fill the ESD with the T0 info (will be used by the TOF)
1187 if (fReconstructor[11])
1188 GetReconstructor(11)->FillESD(fRunLoader, esd);
1189
1190 // pass 1: TPC + ITS inwards
1191 for (Int_t iDet = 1; iDet >= 0; iDet--) {
1192 if (!fTracker[iDet]) continue;
1193 AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
1194
1195 // load clusters
1196 fLoader[iDet]->LoadRecPoints("read");
1197 TTree* tree = fLoader[iDet]->TreeR();
1198 if (!tree) {
1199 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1200 return kFALSE;
1201 }
1202 fTracker[iDet]->LoadClusters(tree);
1203
1204 // run tracking
1205 if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
1206 AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
1207 return kFALSE;
1208 }
1209 if (fCheckPointLevel > 1) {
1210 WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
1211 }
1212 // preliminary PID in TPC needed by the ITS tracker
1213 if (iDet == 1) {
1214 GetReconstructor(1)->FillESD(fRunLoader, esd);
1215 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1216 AliESDpid::MakePID(esd);
1217 }
1218 }
1219
1220 // pass 2: ALL backwards
1221 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1222 if (!fTracker[iDet]) continue;
1223 AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
1224
1225 // load clusters
1226 if (iDet > 1) { // all except ITS, TPC
1227 TTree* tree = NULL;
1228 fLoader[iDet]->LoadRecPoints("read");
1229 tree = fLoader[iDet]->TreeR();
1230 if (!tree) {
1231 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
1232 return kFALSE;
1233 }
1234 fTracker[iDet]->LoadClusters(tree);
1235 }
1236
1237 // run tracking
1238 if (fTracker[iDet]->PropagateBack(esd) != 0) {
1239 AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
1240 return kFALSE;
1241 }
1242 if (fCheckPointLevel > 1) {
1243 WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
1244 }
1245
1246 // unload clusters
1247 if (iDet > 2) { // all except ITS, TPC, TRD
1248 fTracker[iDet]->UnloadClusters();
1249 fLoader[iDet]->UnloadRecPoints();
1250 }
1251 // updated PID in TPC needed by the ITS tracker -MI
1252 if (iDet == 1) {
1253 GetReconstructor(1)->FillESD(fRunLoader, esd);
1254 GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
1255 AliESDpid::MakePID(esd);
1256 }
1257 }
1258
1259 // write space-points to the ESD in case alignment data output
1260 // is switched on
1261 if (fWriteAlignmentData)
1262 WriteAlignmentData(esd);
1263
1264 // pass 3: TRD + TPC + ITS refit inwards
1265 for (Int_t iDet = 2; iDet >= 0; iDet--) {
1266 if (!fTracker[iDet]) continue;
1267 AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
1268
1269 // run tracking
1270 if (fTracker[iDet]->RefitInward(esd) != 0) {
1271 AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
1272 return kFALSE;
1273 }
1274 if (fCheckPointLevel > 1) {
1275 WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
1276 }
1277
1278 // unload clusters
1279 fTracker[iDet]->UnloadClusters();
1280 fLoader[iDet]->UnloadRecPoints();
1281 }
1282 //
1283 // Propagate track to the vertex - if not done by ITS
1284 //
1285 Int_t ntracks = esd->GetNumberOfTracks();
1286 for (Int_t itrack=0; itrack<ntracks; itrack++){
1287 const Double_t kRadius = 3; // beam pipe radius
1288 const Double_t kMaxStep = 5; // max step
1289 const Double_t kMaxD = 123456; // max distance to prim vertex
1290 Double_t fieldZ = AliTracker::GetBz(); //
1291 AliESDtrack * track = esd->GetTrack(itrack);
1292 if (!track) continue;
1293 if (track->IsOn(AliESDtrack::kITSrefit)) continue;
1294 track->PropagateTo(kRadius, fieldZ, track->GetMass(),kMaxStep,kTRUE);
1295 track->RelateToVertex(esd->GetVertex(),fieldZ, kMaxD);
1296 }
1297
1298 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1299 stopwatch.RealTime(),stopwatch.CpuTime()));
1300
1301 return kTRUE;
1302}
1303
1304//_____________________________________________________________________________
1305Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
1306{
1307// fill the event summary data
1308
1309 TStopwatch stopwatch;
1310 stopwatch.Start();
1311 AliInfo("filling ESD");
1312
1313 TString detStr = detectors;
1314 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1315 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1316 AliReconstructor* reconstructor = GetReconstructor(iDet);
1317 if (!reconstructor) continue;
1318
1319 if (!ReadESD(esd, fgkDetectorName[iDet])) {
1320 AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
1321 TTree* clustersTree = NULL;
1322 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1323 fLoader[iDet]->LoadRecPoints("read");
1324 clustersTree = fLoader[iDet]->TreeR();
1325 if (!clustersTree) {
1326 AliError(Form("Can't get the %s clusters tree",
1327 fgkDetectorName[iDet]));
1328 if (fStopOnError) return kFALSE;
1329 }
1330 }
1331 if (fRawReader && !reconstructor->HasDigitConversion()) {
1332 reconstructor->FillESD(fRawReader, clustersTree, esd);
1333 } else {
1334 TTree* digitsTree = NULL;
1335 if (fLoader[iDet]) {
1336 fLoader[iDet]->LoadDigits("read");
1337 digitsTree = fLoader[iDet]->TreeD();
1338 if (!digitsTree) {
1339 AliError(Form("Can't get the %s digits tree",
1340 fgkDetectorName[iDet]));
1341 if (fStopOnError) return kFALSE;
1342 }
1343 }
1344 reconstructor->FillESD(digitsTree, clustersTree, esd);
1345 if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
1346 }
1347 if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
1348 fLoader[iDet]->UnloadRecPoints();
1349 }
1350
1351 if (fRawReader) {
1352 reconstructor->FillESD(fRunLoader, fRawReader, esd);
1353 } else {
1354 reconstructor->FillESD(fRunLoader, esd);
1355 }
1356 if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
1357 }
1358 }
1359
1360 if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
1361 AliError(Form("the following detectors were not found: %s",
1362 detStr.Data()));
1363 if (fStopOnError) return kFALSE;
1364 }
1365
1366 AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
1367 stopwatch.RealTime(),stopwatch.CpuTime()));
1368
1369 return kTRUE;
1370}
1371
1372//_____________________________________________________________________________
1373Bool_t AliReconstruction::FillTriggerESD(AliESD*& esd)
1374{
1375 // Reads the trigger decision which is
1376 // stored in Trigger.root file and fills
1377 // the corresponding esd entries
1378
1379 AliInfo("Filling trigger information into the ESD");
1380
1381 if (fRawReader) {
1382 AliCTPRawStream input(fRawReader);
1383 if (!input.Next()) {
1384 AliError("No valid CTP (trigger) DDL raw data is found ! The trigger information is not stored in the ESD !");
1385 return kFALSE;
1386 }
1387 esd->SetTriggerMask(input.GetClassMask());
1388 esd->SetTriggerCluster(input.GetClusterMask());
1389 }
1390 else {
1391 AliRunLoader *runloader = AliRunLoader::GetRunLoader();
1392 if (runloader) {
1393 if (!runloader->LoadTrigger()) {
1394 AliCentralTrigger *aCTP = runloader->GetTrigger();
1395 esd->SetTriggerMask(aCTP->GetClassMask());
1396 esd->SetTriggerCluster(aCTP->GetClusterMask());
1397 }
1398 else {
1399 AliWarning("No trigger can be loaded! The trigger information is not stored in the ESD !");
1400 return kFALSE;
1401 }
1402 }
1403 else {
1404 AliError("No run loader is available! The trigger information is not stored in the ESD !");
1405 return kFALSE;
1406 }
1407 }
1408
1409 return kTRUE;
1410}
1411
1412
1413
1414
1415
1416//_____________________________________________________________________________
1417Bool_t AliReconstruction::FillRawEventHeaderESD(AliESD*& esd)
1418{
1419 //
1420 // Filling information from RawReader Header
1421 //
1422
1423 AliInfo("Filling information from RawReader Header");
1424 esd->SetBunchCrossNumber(0);
1425 esd->SetOrbitNumber(0);
1426 esd->SetPeriodNumber(0);
1427 esd->SetTimeStamp(0);
1428 esd->SetEventType(0);
1429 const AliRawEventHeaderBase * eventHeader = fRawReader->GetEventHeader();
1430 if (eventHeader){
1431
1432 const UInt_t *id = eventHeader->GetP("Id");
1433 esd->SetBunchCrossNumber((id)[1]&0x00000fff);
1434 esd->SetOrbitNumber((((id)[0]<<20)&0xf00000)|(((id)[1]>>12)&0xfffff));
1435 esd->SetPeriodNumber(((id)[0]>>4)&0x0fffffff);
1436
1437 esd->SetTimeStamp((eventHeader->Get("Timestamp")));
1438 esd->SetEventType((eventHeader->Get("Type")));
1439 }
1440
1441 return kTRUE;
1442}
1443
1444
1445//_____________________________________________________________________________
1446Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
1447{
1448// check whether detName is contained in detectors
1449// if yes, it is removed from detectors
1450
1451 // check if all detectors are selected
1452 if ((detectors.CompareTo("ALL") == 0) ||
1453 detectors.BeginsWith("ALL ") ||
1454 detectors.EndsWith(" ALL") ||
1455 detectors.Contains(" ALL ")) {
1456 detectors = "ALL";
1457 return kTRUE;
1458 }
1459
1460 // search for the given detector
1461 Bool_t result = kFALSE;
1462 if ((detectors.CompareTo(detName) == 0) ||
1463 detectors.BeginsWith(detName+" ") ||
1464 detectors.EndsWith(" "+detName) ||
1465 detectors.Contains(" "+detName+" ")) {
1466 detectors.ReplaceAll(detName, "");
1467 result = kTRUE;
1468 }
1469
1470 // clean up the detectors string
1471 while (detectors.Contains(" ")) detectors.ReplaceAll(" ", " ");
1472 while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
1473 while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
1474
1475 return result;
1476}
1477
1478//_____________________________________________________________________________
1479Bool_t AliReconstruction::InitRunLoader()
1480{
1481// get or create the run loader
1482
1483 if (gAlice) delete gAlice;
1484 gAlice = NULL;
1485
1486 if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
1487 // load all base libraries to get the loader classes
1488 TString libs = gSystem->GetLibraries();
1489 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1490 TString detName = fgkDetectorName[iDet];
1491 if (detName == "HLT") continue;
1492 if (libs.Contains("lib" + detName + "base.so")) continue;
1493 gSystem->Load("lib" + detName + "base.so");
1494 }
1495 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
1496 if (!fRunLoader) {
1497 AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
1498 CleanUp();
1499 return kFALSE;
1500 }
1501 fRunLoader->CdGAFile();
1502 if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
1503 if (fRunLoader->LoadgAlice() == 0) {
1504 gAlice = fRunLoader->GetAliRun();
1505 AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
1506 }
1507 }
1508 if (!gAlice && !fRawReader) {
1509 AliError(Form("no gAlice object found in file %s",
1510 fGAliceFileName.Data()));
1511 CleanUp();
1512 return kFALSE;
1513 }
1514
1515 } else { // galice.root does not exist
1516 if (!fRawReader) {
1517 AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
1518 CleanUp();
1519 return kFALSE;
1520 }
1521 fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
1522 AliConfig::GetDefaultEventFolderName(),
1523 "recreate");
1524 if (!fRunLoader) {
1525 AliError(Form("could not create run loader in file %s",
1526 fGAliceFileName.Data()));
1527 CleanUp();
1528 return kFALSE;
1529 }
1530 fRunLoader->MakeTree("E");
1531 Int_t iEvent = 0;
1532 while (fRawReader->NextEvent()) {
1533 fRunLoader->SetEventNumber(iEvent);
1534 fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
1535 iEvent, iEvent);
1536 fRunLoader->MakeTree("H");
1537 fRunLoader->TreeE()->Fill();
1538 iEvent++;
1539 }
1540 fRawReader->RewindEvents();
1541 if (fNumberOfEventsPerFile > 0)
1542 fRunLoader->SetNumberOfEventsPerFile(fNumberOfEventsPerFile);
1543 else
1544 fRunLoader->SetNumberOfEventsPerFile(iEvent);
1545 fRunLoader->WriteHeader("OVERWRITE");
1546 fRunLoader->CdGAFile();
1547 fRunLoader->Write(0, TObject::kOverwrite);
1548// AliTracker::SetFieldMap(???);
1549 }
1550
1551 return kTRUE;
1552}
1553
1554//_____________________________________________________________________________
1555AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
1556{
1557// get the reconstructor object and the loader for a detector
1558
1559 if (fReconstructor[iDet]) return fReconstructor[iDet];
1560
1561 // load the reconstructor object
1562 TPluginManager* pluginManager = gROOT->GetPluginManager();
1563 TString detName = fgkDetectorName[iDet];
1564 TString recName = "Ali" + detName + "Reconstructor";
1565 if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
1566
1567 if (detName == "HLT") {
1568 if (!gROOT->GetClass("AliLevel3")) {
1569 gSystem->Load("libAliHLTSrc.so");
1570 gSystem->Load("libAliHLTMisc.so");
1571 gSystem->Load("libAliHLTHough.so");
1572 gSystem->Load("libAliHLTComp.so");
1573 }
1574 }
1575
1576 AliReconstructor* reconstructor = NULL;
1577 // first check if a plugin is defined for the reconstructor
1578 TPluginHandler* pluginHandler =
1579 pluginManager->FindHandler("AliReconstructor", detName);
1580 // if not, add a plugin for it
1581 if (!pluginHandler) {
1582 AliDebug(1, Form("defining plugin for %s", recName.Data()));
1583 TString libs = gSystem->GetLibraries();
1584 if (libs.Contains("lib" + detName + "base.so") ||
1585 (gSystem->Load("lib" + detName + "base.so") >= 0)) {
1586 pluginManager->AddHandler("AliReconstructor", detName,
1587 recName, detName + "rec", recName + "()");
1588 } else {
1589 pluginManager->AddHandler("AliReconstructor", detName,
1590 recName, detName, recName + "()");
1591 }
1592 pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
1593 }
1594 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1595 reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
1596 }
1597 if (reconstructor) {
1598 TObject* obj = fOptions.FindObject(detName.Data());
1599 if (obj) reconstructor->SetOption(obj->GetTitle());
1600 reconstructor->Init(fRunLoader);
1601 fReconstructor[iDet] = reconstructor;
1602 }
1603
1604 // get or create the loader
1605 if (detName != "HLT") {
1606 fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
1607 if (!fLoader[iDet]) {
1608 AliConfig::Instance()
1609 ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
1610 detName, detName);
1611 // first check if a plugin is defined for the loader
1612 pluginHandler =
1613 pluginManager->FindHandler("AliLoader", detName);
1614 // if not, add a plugin for it
1615 if (!pluginHandler) {
1616 TString loaderName = "Ali" + detName + "Loader";
1617 AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1618 pluginManager->AddHandler("AliLoader", detName,
1619 loaderName, detName + "base",
1620 loaderName + "(const char*, TFolder*)");
1621 pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1622 }
1623 if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1624 fLoader[iDet] =
1625 (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1626 fRunLoader->GetEventFolder());
1627 }
1628 if (!fLoader[iDet]) { // use default loader
1629 fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1630 }
1631 if (!fLoader[iDet]) {
1632 AliWarning(Form("couldn't get loader for %s", detName.Data()));
1633 if (fStopOnError) return NULL;
1634 } else {
1635 fRunLoader->AddLoader(fLoader[iDet]);
1636 fRunLoader->CdGAFile();
1637 if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1638 fRunLoader->Write(0, TObject::kOverwrite);
1639 }
1640 }
1641 }
1642
1643 return reconstructor;
1644}
1645
1646//_____________________________________________________________________________
1647Bool_t AliReconstruction::CreateVertexer()
1648{
1649// create the vertexer
1650
1651 fVertexer = NULL;
1652 AliReconstructor* itsReconstructor = GetReconstructor(0);
1653 if (itsReconstructor) {
1654 fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1655 }
1656 if (!fVertexer) {
1657 AliWarning("couldn't create a vertexer for ITS");
1658 if (fStopOnError) return kFALSE;
1659 }
1660
1661 return kTRUE;
1662}
1663
1664//_____________________________________________________________________________
1665Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1666{
1667// create the trackers
1668
1669 TString detStr = detectors;
1670 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1671 if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1672 AliReconstructor* reconstructor = GetReconstructor(iDet);
1673 if (!reconstructor) continue;
1674 TString detName = fgkDetectorName[iDet];
1675 if (detName == "HLT") {
1676 fRunHLTTracking = kTRUE;
1677 continue;
1678 }
1679 if (detName == "MUON") {
1680 fRunMuonTracking = kTRUE;
1681 continue;
1682 }
1683
1684
1685 fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1686 if (!fTracker[iDet] && (iDet < 7)) {
1687 AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1688 if (fStopOnError) return kFALSE;
1689 }
1690 }
1691
1692 return kTRUE;
1693}
1694
1695//_____________________________________________________________________________
1696void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1697{
1698// delete trackers and the run loader and close and delete the file
1699
1700 for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1701 delete fReconstructor[iDet];
1702 fReconstructor[iDet] = NULL;
1703 fLoader[iDet] = NULL;
1704 delete fTracker[iDet];
1705 fTracker[iDet] = NULL;
1706 }
1707 delete fVertexer;
1708 fVertexer = NULL;
1709 delete fDiamondProfile;
1710 fDiamondProfile = NULL;
1711
1712 delete fRunLoader;
1713 fRunLoader = NULL;
1714 delete fRawReader;
1715 fRawReader = NULL;
1716
1717 if (file) {
1718 file->Close();
1719 delete file;
1720 }
1721
1722 if (fileOld) {
1723 fileOld->Close();
1724 delete fileOld;
1725 gSystem->Unlink("AliESDs.old.root");
1726 }
1727}
1728
1729
1730//_____________________________________________________________________________
1731Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1732{
1733// read the ESD event from a file
1734
1735 if (!esd) return kFALSE;
1736 char fileName[256];
1737 sprintf(fileName, "ESD_%d.%d_%s.root",
1738 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1739 if (gSystem->AccessPathName(fileName)) return kFALSE;
1740
1741 AliInfo(Form("reading ESD from file %s", fileName));
1742 AliDebug(1, Form("reading ESD from file %s", fileName));
1743 TFile* file = TFile::Open(fileName);
1744 if (!file || !file->IsOpen()) {
1745 AliError(Form("opening %s failed", fileName));
1746 delete file;
1747 return kFALSE;
1748 }
1749
1750 gROOT->cd();
1751 delete esd;
1752 esd = (AliESD*) file->Get("ESD");
1753 file->Close();
1754 delete file;
1755 return kTRUE;
1756}
1757
1758//_____________________________________________________________________________
1759void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1760{
1761// write the ESD event to a file
1762
1763 if (!esd) return;
1764 char fileName[256];
1765 sprintf(fileName, "ESD_%d.%d_%s.root",
1766 esd->GetRunNumber(), esd->GetEventNumberInFile(), recStep);
1767
1768 AliDebug(1, Form("writing ESD to file %s", fileName));
1769 TFile* file = TFile::Open(fileName, "recreate");
1770 if (!file || !file->IsOpen()) {
1771 AliError(Form("opening %s failed", fileName));
1772 } else {
1773 esd->Write("ESD");
1774 file->Close();
1775 }
1776 delete file;
1777}
1778
1779
1780
1781
1782//_____________________________________________________________________________
1783void AliReconstruction::CreateTag(TFile* file)
1784{
1785 //GRP
1786 Float_t lhcLuminosity = 0.0;
1787 TString lhcState = "test";
1788 UInt_t detectorMask = 0;
1789
1790 /////////////
1791 //muon code//
1792 ////////////
1793 Double_t fMUONMASS = 0.105658369;
1794 //Variables
1795 Double_t fX,fY,fZ ;
1796 Double_t fThetaX, fThetaY, fPyz, fChisquare;
1797 Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1798 Int_t fCharge;
1799 TLorentzVector fEPvector;
1800
1801 Float_t fZVertexCut = 10.0;
1802 Float_t fRhoVertexCut = 2.0;
1803
1804 Float_t fLowPtCut = 1.0;
1805 Float_t fHighPtCut = 3.0;
1806 Float_t fVeryHighPtCut = 10.0;
1807 ////////////
1808
1809 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1810
1811 // Creates the tags for all the events in a given ESD file
1812 Int_t ntrack;
1813 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1814 Int_t nPos, nNeg, nNeutr;
1815 Int_t nK0s, nNeutrons, nPi0s, nGamas;
1816 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1817 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1818 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1819 Float_t maxPt = .0, meanPt = .0, totalP = .0;
1820 Int_t fVertexflag;
1821 Int_t iRunNumber = 0;
1822 TString fVertexName("default");
1823
1824 AliRunTag *tag = new AliRunTag();
1825 AliEventTag *evTag = new AliEventTag();
1826 TTree ttag("T","A Tree with event tags");
1827 TBranch * btag = ttag.Branch("AliTAG", &tag);
1828 btag->SetCompressionLevel(9);
1829
1830 AliInfo(Form("Creating the tags......."));
1831
1832 if (!file || !file->IsOpen()) {
1833 AliError(Form("opening failed"));
1834 delete file;
1835 return ;
1836 }
1837 Int_t lastEvent = 0;
1838 TTree *t = (TTree*) file->Get("esdTree");
1839 TBranch * b = t->GetBranch("ESD");
1840 AliESD *esd = 0;
1841 b->SetAddress(&esd);
1842
1843 b->GetEntry(fFirstEvent);
1844 Int_t iInitRunNumber = esd->GetRunNumber();
1845
1846 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1847 if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1848 for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1849 ntrack = 0;
1850 nPos = 0;
1851 nNeg = 0;
1852 nNeutr =0;
1853 nK0s = 0;
1854 nNeutrons = 0;
1855 nPi0s = 0;
1856 nGamas = 0;
1857 nProtons = 0;
1858 nKaons = 0;
1859 nPions = 0;
1860 nMuons = 0;
1861 nElectrons = 0;
1862 nCh1GeV = 0;
1863 nCh3GeV = 0;
1864 nCh10GeV = 0;
1865 nMu1GeV = 0;
1866 nMu3GeV = 0;
1867 nMu10GeV = 0;
1868 nEl1GeV = 0;
1869 nEl3GeV = 0;
1870 nEl10GeV = 0;
1871 maxPt = .0;
1872 meanPt = .0;
1873 totalP = .0;
1874 fVertexflag = 0;
1875
1876 b->GetEntry(iEventNumber);
1877 iRunNumber = esd->GetRunNumber();
1878 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1879 const AliESDVertex * vertexIn = esd->GetVertex();
1880 if (!vertexIn) AliError("ESD has not defined vertex.");
1881 if (vertexIn) fVertexName = vertexIn->GetName();
1882 if(fVertexName != "default") fVertexflag = 1;
1883 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1884 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1885 UInt_t status = esdTrack->GetStatus();
1886
1887 //select only tracks with ITS refit
1888 if ((status&AliESDtrack::kITSrefit)==0) continue;
1889 //select only tracks with TPC refit
1890 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1891
1892 //select only tracks with the "combined PID"
1893 if ((status&AliESDtrack::kESDpid)==0) continue;
1894 Double_t p[3];
1895 esdTrack->GetPxPyPz(p);
1896 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1897 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1898 totalP += momentum;
1899 meanPt += fPt;
1900 if(fPt > maxPt) maxPt = fPt;
1901
1902 if(esdTrack->GetSign() > 0) {
1903 nPos++;
1904 if(fPt > fLowPtCut) nCh1GeV++;
1905 if(fPt > fHighPtCut) nCh3GeV++;
1906 if(fPt > fVeryHighPtCut) nCh10GeV++;
1907 }
1908 if(esdTrack->GetSign() < 0) {
1909 nNeg++;
1910 if(fPt > fLowPtCut) nCh1GeV++;
1911 if(fPt > fHighPtCut) nCh3GeV++;
1912 if(fPt > fVeryHighPtCut) nCh10GeV++;
1913 }
1914 if(esdTrack->GetSign() == 0) nNeutr++;
1915
1916 //PID
1917 Double_t prob[5];
1918 esdTrack->GetESDpid(prob);
1919
1920 Double_t rcc = 0.0;
1921 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1922 if(rcc == 0.0) continue;
1923 //Bayes' formula
1924 Double_t w[5];
1925 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1926
1927 //protons
1928 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1929 //kaons
1930 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1931 //pions
1932 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1933 //electrons
1934 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1935 nElectrons++;
1936 if(fPt > fLowPtCut) nEl1GeV++;
1937 if(fPt > fHighPtCut) nEl3GeV++;
1938 if(fPt > fVeryHighPtCut) nEl10GeV++;
1939 }
1940 ntrack++;
1941 }//track loop
1942
1943 /////////////
1944 //muon code//
1945 ////////////
1946 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1947 // loop over all reconstructed tracks (also first track of combination)
1948 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
1949 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1950 if (muonTrack == 0x0) continue;
1951
1952 // Coordinates at vertex
1953 fZ = muonTrack->GetZ();
1954 fY = muonTrack->GetBendingCoor();
1955 fX = muonTrack->GetNonBendingCoor();
1956
1957 fThetaX = muonTrack->GetThetaX();
1958 fThetaY = muonTrack->GetThetaY();
1959
1960 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1961 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1962 fPxRec = fPzRec * TMath::Tan(fThetaX);
1963 fPyRec = fPzRec * TMath::Tan(fThetaY);
1964 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1965
1966 //ChiSquare of the track if needed
1967 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1968 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1969 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1970
1971 // total number of muons inside a vertex cut
1972 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1973 nMuons++;
1974 if(fEPvector.Pt() > fLowPtCut) {
1975 nMu1GeV++;
1976 if(fEPvector.Pt() > fHighPtCut) {
1977 nMu3GeV++;
1978 if (fEPvector.Pt() > fVeryHighPtCut) {
1979 nMu10GeV++;
1980 }
1981 }
1982 }
1983 }
1984 }//muon track loop
1985
1986 // Fill the event tags
1987 if(ntrack != 0)
1988 meanPt = meanPt/ntrack;
1989
1990 evTag->SetEventId(iEventNumber+1);
1991 if (vertexIn) {
1992 evTag->SetVertexX(vertexIn->GetXv());
1993 evTag->SetVertexY(vertexIn->GetYv());
1994 evTag->SetVertexZ(vertexIn->GetZv());
1995 evTag->SetVertexZError(vertexIn->GetZRes());
1996 }
1997 evTag->SetVertexFlag(fVertexflag);
1998
1999 evTag->SetT0VertexZ(esd->GetT0zVertex());
2000
2001 evTag->SetTriggerMask(esd->GetTriggerMask());
2002 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2003
2004 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
2005 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
2006 evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
2007 evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
2008 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
2009 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
2010
2011
2012 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
2013 evTag->SetNumOfPosTracks(nPos);
2014 evTag->SetNumOfNegTracks(nNeg);
2015 evTag->SetNumOfNeutrTracks(nNeutr);
2016
2017 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
2018 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
2019 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
2020 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
2021
2022 evTag->SetNumOfProtons(nProtons);
2023 evTag->SetNumOfKaons(nKaons);
2024 evTag->SetNumOfPions(nPions);
2025 evTag->SetNumOfMuons(nMuons);
2026 evTag->SetNumOfElectrons(nElectrons);
2027 evTag->SetNumOfPhotons(nGamas);
2028 evTag->SetNumOfPi0s(nPi0s);
2029 evTag->SetNumOfNeutrons(nNeutrons);
2030 evTag->SetNumOfKaon0s(nK0s);
2031
2032 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
2033 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
2034 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
2035 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
2036 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
2037 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
2038 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
2039 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
2040 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
2041
2042 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
2043 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2044
2045 evTag->SetTotalMomentum(totalP);
2046 evTag->SetMeanPt(meanPt);
2047 evTag->SetMaxPt(maxPt);
2048
2049 tag->SetLHCTag(lhcLuminosity,lhcState);
2050 tag->SetDetectorTag(detectorMask);
2051
2052 tag->SetRunId(iInitRunNumber);
2053 tag->AddEventTag(*evTag);
2054 }
2055 if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
2056 else lastEvent = fLastEvent;
2057
2058 ttag.Fill();
2059 tag->Clear();
2060
2061 char fileName[256];
2062 sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
2063 tag->GetRunId(),fFirstEvent,lastEvent );
2064 AliInfo(Form("writing tags to file %s", fileName));
2065 AliDebug(1, Form("writing tags to file %s", fileName));
2066
2067 TFile* ftag = TFile::Open(fileName, "recreate");
2068 ftag->cd();
2069 ttag.Write();
2070 ftag->Close();
2071 file->cd();
2072 delete tag;
2073 delete evTag;
2074}
2075
2076//_____________________________________________________________________________
2077void AliReconstruction::ESDFile2AODFile(TFile* esdFile, TFile* aodFile)
2078{
2079 // write all files from the given esd file to an aod file
2080
2081 // create an AliAOD object
2082 AliAODEvent *aod = new AliAODEvent();
2083 aod->CreateStdContent();
2084
2085 // go to the file
2086 aodFile->cd();
2087
2088 // create the tree
2089 TTree *aodTree = new TTree("AOD", "AliAOD tree");
2090 aodTree->Branch(aod->GetList());
2091
2092 // connect to ESD
2093 TTree *t = (TTree*) esdFile->Get("esdTree");
2094 TBranch *b = t->GetBranch("ESD");
2095 AliESD *esd = 0;
2096 b->SetAddress(&esd);
2097
2098 Int_t nEvents = b->GetEntries();
2099
2100 // set arrays and pointers
2101 Float_t posF[3];
2102 Double_t pos[3];
2103 Double_t p[3];
2104 Double_t covVtx[6];
2105 Double_t covTr[21];
2106 Double_t pid[10];
2107
2108 // loop over events and fill them
2109 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
2110 b->GetEntry(iEvent);
2111
2112 // Multiplicity information needed by the header (to be revised!)
2113 Int_t nTracks = esd->GetNumberOfTracks();
2114 Int_t nPosTracks = 0;
2115 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
2116 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
2117
2118 // create the header
2119 aod->AddHeader(new AliAODHeader(esd->GetRunNumber(),
2120 esd->GetBunchCrossNumber(),
2121 esd->GetOrbitNumber(),
2122 esd->GetPeriodNumber(),
2123 nTracks,
2124 nPosTracks,
2125 nTracks-nPosTracks,
2126 esd->GetMagneticField(),
2127 -999., // fill muon magnetic field
2128 -999., // centrality; to be filled, still
2129 esd->GetZDCN1Energy(),
2130 esd->GetZDCP1Energy(),
2131 esd->GetZDCN2Energy(),
2132 esd->GetZDCP2Energy(),
2133 esd->GetZDCEMEnergy(),
2134 esd->GetTriggerMask(),
2135 esd->GetTriggerCluster(),
2136 esd->GetEventType()));
2137
2138 Int_t nV0s = esd->GetNumberOfV0s();
2139 Int_t nCascades = esd->GetNumberOfCascades();
2140 Int_t nKinks = esd->GetNumberOfKinks();
2141 Int_t nVertices = nV0s + nCascades + nKinks;
2142
2143 aod->ResetStd(nTracks, nVertices);
2144 AliAODTrack *aodTrack;
2145
2146
2147 // Array to take into account the tracks already added to the AOD
2148 Bool_t * usedTrack = NULL;
2149 if (nTracks>0) {
2150 usedTrack = new Bool_t[nTracks];
2151 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
2152 }
2153 // Array to take into account the V0s already added to the AOD
2154 Bool_t * usedV0 = NULL;
2155 if (nV0s>0) {
2156 usedV0 = new Bool_t[nV0s];
2157 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
2158 }
2159 // Array to take into account the kinks already added to the AOD
2160 Bool_t * usedKink = NULL;
2161 if (nKinks>0) {
2162 usedKink = new Bool_t[nKinks];
2163 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
2164 }
2165
2166 // Access to the AOD container of vertices
2167 TClonesArray &vertices = *(aod->GetVertices());
2168 Int_t jVertices=0;
2169
2170 // Access to the AOD container of tracks
2171 TClonesArray &tracks = *(aod->GetTracks());
2172 Int_t jTracks=0;
2173
2174 // Add primary vertex. The primary tracks will be defined
2175 // after the loops on the composite objects (V0, cascades, kinks)
2176 const AliESDVertex *vtx = esd->GetPrimaryVertex();
2177
2178 vtx->GetXYZ(pos); // position
2179 vtx->GetCovMatrix(covVtx); //covariance matrix
2180
2181 AliAODVertex * primary = new(vertices[jVertices++])
2182 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
2183
2184 // Create vertices starting from the most complex objects
2185
2186 // Cascades
2187 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
2188 AliESDcascade *cascade = esd->GetCascade(nCascade);
2189
2190 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2191 cascade->GetPosCovXi(covVtx);
2192
2193 // Add the cascade vertex
2194 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
2195 covVtx,
2196 cascade->GetChi2Xi(), // = chi2/NDF since NDF = 2*2-3
2197 primary,
2198 AliAODVertex::kCascade);
2199
2200 primary->AddDaughter(vcascade);
2201
2202 // Add the V0 from the cascade. The ESD class have to be optimized...
2203 // Now we have to search for the corresponding Vo in the list of V0s
2204 // using the indeces of the positive and negative tracks
2205
2206 Int_t posFromV0 = cascade->GetPindex();
2207 Int_t negFromV0 = cascade->GetNindex();
2208
2209
2210 AliESDv0 * v0 = 0x0;
2211 Int_t indV0 = -1;
2212
2213 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
2214
2215 v0 = esd->GetV0(iV0);
2216 Int_t posV0 = v0->GetPindex();
2217 Int_t negV0 = v0->GetNindex();
2218
2219 if (posV0==posFromV0 && negV0==negFromV0) {
2220 indV0 = iV0;
2221 break;
2222 }
2223 }
2224
2225 AliAODVertex * vV0FromCascade = 0x0;
2226
2227 if (indV0>-1 && !usedV0[indV0] ) {
2228
2229 // the V0 exists in the array of V0s and is not used
2230
2231 usedV0[indV0] = kTRUE;
2232
2233 v0->GetXYZ(pos[0], pos[1], pos[2]);
2234 v0->GetPosCov(covVtx);
2235
2236 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2237 covVtx,
2238 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2239 vcascade,
2240 AliAODVertex::kV0);
2241 } else {
2242
2243 // the V0 doesn't exist in the array of V0s or was used
2244 cerr << "Error: event " << iEvent << " cascade " << nCascade
2245 << " The V0 " << indV0
2246 << " doesn't exist in the array of V0s or was used!" << endl;
2247
2248 cascade->GetXYZ(pos[0], pos[1], pos[2]);
2249 cascade->GetPosCov(covVtx);
2250
2251 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
2252 covVtx,
2253 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2254 vcascade,
2255 AliAODVertex::kV0);
2256 vcascade->AddDaughter(vV0FromCascade);
2257 }
2258
2259 // Add the positive tracks from the V0
2260
2261 if (! usedTrack[posFromV0]) {
2262
2263 usedTrack[posFromV0] = kTRUE;
2264
2265 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2266 esdTrack->GetPxPyPz(p);
2267 esdTrack->GetXYZ(pos);
2268 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2269 esdTrack->GetESDpid(pid);
2270
2271 vV0FromCascade->AddDaughter(aodTrack =
2272 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2273 esdTrack->GetLabel(),
2274 p,
2275 kTRUE,
2276 pos,
2277 kFALSE,
2278 covTr,
2279 (Short_t)esdTrack->GetSign(),
2280 esdTrack->GetITSClusterMap(),
2281 pid,
2282 vV0FromCascade,
2283 kTRUE, // check if this is right
2284 kFALSE, // check if this is right
2285 AliAODTrack::kSecondary)
2286 );
2287 aodTrack->ConvertAliPIDtoAODPID();
2288 }
2289 else {
2290 cerr << "Error: event " << iEvent << " cascade " << nCascade
2291 << " track " << posFromV0 << " has already been used!" << endl;
2292 }
2293
2294 // Add the negative tracks from the V0
2295
2296 if (!usedTrack[negFromV0]) {
2297
2298 usedTrack[negFromV0] = kTRUE;
2299
2300 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2301 esdTrack->GetPxPyPz(p);
2302 esdTrack->GetXYZ(pos);
2303 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2304 esdTrack->GetESDpid(pid);
2305
2306 vV0FromCascade->AddDaughter(aodTrack =
2307 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2308 esdTrack->GetLabel(),
2309 p,
2310 kTRUE,
2311 pos,
2312 kFALSE,
2313 covTr,
2314 (Short_t)esdTrack->GetSign(),
2315 esdTrack->GetITSClusterMap(),
2316 pid,
2317 vV0FromCascade,
2318 kTRUE, // check if this is right
2319 kFALSE, // check if this is right
2320 AliAODTrack::kSecondary)
2321 );
2322 aodTrack->ConvertAliPIDtoAODPID();
2323 }
2324 else {
2325 cerr << "Error: event " << iEvent << " cascade " << nCascade
2326 << " track " << negFromV0 << " has already been used!" << endl;
2327 }
2328
2329 // Add the bachelor track from the cascade
2330
2331 Int_t bachelor = cascade->GetBindex();
2332
2333 if(!usedTrack[bachelor]) {
2334
2335 usedTrack[bachelor] = kTRUE;
2336
2337 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
2338 esdTrack->GetPxPyPz(p);
2339 esdTrack->GetXYZ(pos);
2340 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2341 esdTrack->GetESDpid(pid);
2342
2343 vcascade->AddDaughter(aodTrack =
2344 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2345 esdTrack->GetLabel(),
2346 p,
2347 kTRUE,
2348 pos,
2349 kFALSE,
2350 covTr,
2351 (Short_t)esdTrack->GetSign(),
2352 esdTrack->GetITSClusterMap(),
2353 pid,
2354 vcascade,
2355 kTRUE, // check if this is right
2356 kFALSE, // check if this is right
2357 AliAODTrack::kSecondary)
2358 );
2359 aodTrack->ConvertAliPIDtoAODPID();
2360 }
2361 else {
2362 cerr << "Error: event " << iEvent << " cascade " << nCascade
2363 << " track " << bachelor << " has already been used!" << endl;
2364 }
2365
2366 // Add the primary track of the cascade (if any)
2367
2368 } // end of the loop on cascades
2369
2370 // V0s
2371
2372 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
2373
2374 if (usedV0[nV0]) continue; // skip if aready added to the AOD
2375
2376 AliESDv0 *v0 = esd->GetV0(nV0);
2377
2378 v0->GetXYZ(pos[0], pos[1], pos[2]);
2379 v0->GetPosCov(covVtx);
2380
2381 AliAODVertex * vV0 =
2382 new(vertices[jVertices++]) AliAODVertex(pos,
2383 covVtx,
2384 v0->GetChi2V0(), // = chi2/NDF since NDF = 2*2-3
2385 primary,
2386 AliAODVertex::kV0);
2387 primary->AddDaughter(vV0);
2388
2389 Int_t posFromV0 = v0->GetPindex();
2390 Int_t negFromV0 = v0->GetNindex();
2391
2392 // Add the positive tracks from the V0
2393
2394 if (!usedTrack[posFromV0]) {
2395
2396 usedTrack[posFromV0] = kTRUE;
2397
2398 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
2399 esdTrack->GetPxPyPz(p);
2400 esdTrack->GetXYZ(pos);
2401 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2402 esdTrack->GetESDpid(pid);
2403
2404 vV0->AddDaughter(aodTrack =
2405 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2406 esdTrack->GetLabel(),
2407 p,
2408 kTRUE,
2409 pos,
2410 kFALSE,
2411 covTr,
2412 (Short_t)esdTrack->GetSign(),
2413 esdTrack->GetITSClusterMap(),
2414 pid,
2415 vV0,
2416 kTRUE, // check if this is right
2417 kFALSE, // check if this is right
2418 AliAODTrack::kSecondary)
2419 );
2420 aodTrack->ConvertAliPIDtoAODPID();
2421 }
2422 else {
2423 cerr << "Error: event " << iEvent << " V0 " << nV0
2424 << " track " << posFromV0 << " has already been used!" << endl;
2425 }
2426
2427 // Add the negative tracks from the V0
2428
2429 if (!usedTrack[negFromV0]) {
2430
2431 usedTrack[negFromV0] = kTRUE;
2432
2433 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
2434 esdTrack->GetPxPyPz(p);
2435 esdTrack->GetXYZ(pos);
2436 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2437 esdTrack->GetESDpid(pid);
2438
2439 vV0->AddDaughter(aodTrack =
2440 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2441 esdTrack->GetLabel(),
2442 p,
2443 kTRUE,
2444 pos,
2445 kFALSE,
2446 covTr,
2447 (Short_t)esdTrack->GetSign(),
2448 esdTrack->GetITSClusterMap(),
2449 pid,
2450 vV0,
2451 kTRUE, // check if this is right
2452 kFALSE, // check if this is right
2453 AliAODTrack::kSecondary)
2454 );
2455 aodTrack->ConvertAliPIDtoAODPID();
2456 }
2457 else {
2458 cerr << "Error: event " << iEvent << " V0 " << nV0
2459 << " track " << negFromV0 << " has already been used!" << endl;
2460 }
2461
2462 } // end of the loop on V0s
2463
2464 // Kinks: it is a big mess the access to the information in the kinks
2465 // The loop is on the tracks in order to find the mother and daugther of each kink
2466
2467
2468 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
2469
2470
2471 AliESDtrack * esdTrack = esd->GetTrack(iTrack);
2472
2473 Int_t ikink = esdTrack->GetKinkIndex(0);
2474
2475 if (ikink) {
2476 // Negative kink index: mother, positive: daughter
2477
2478 // Search for the second track of the kink
2479
2480 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
2481
2482 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
2483
2484 Int_t jkink = esdTrack1->GetKinkIndex(0);
2485
2486 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
2487
2488 // The two tracks are from the same kink
2489
2490 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
2491
2492 Int_t imother = -1;
2493 Int_t idaughter = -1;
2494
2495 if (ikink<0 && jkink>0) {
2496
2497 imother = iTrack;
2498 idaughter = jTrack;
2499 }
2500 else if (ikink>0 && jkink<0) {
2501
2502 imother = jTrack;
2503 idaughter = iTrack;
2504 }
2505 else {
2506 cerr << "Error: Wrong combination of kink indexes: "
2507 << ikink << " " << jkink << endl;
2508 continue;
2509 }
2510
2511 // Add the mother track
2512
2513 AliAODTrack * mother = NULL;
2514
2515 if (!usedTrack[imother]) {
2516
2517 usedTrack[imother] = kTRUE;
2518
2519 AliESDtrack *esdTrack = esd->GetTrack(imother);
2520 esdTrack->GetPxPyPz(p);
2521 esdTrack->GetXYZ(pos);
2522 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2523 esdTrack->GetESDpid(pid);
2524
2525 mother =
2526 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2527 esdTrack->GetLabel(),
2528 p,
2529 kTRUE,
2530 pos,
2531 kFALSE,
2532 covTr,
2533 (Short_t)esdTrack->GetSign(),
2534 esdTrack->GetITSClusterMap(),
2535 pid,
2536 primary,
2537 kTRUE, // check if this is right
2538 kTRUE, // check if this is right
2539 AliAODTrack::kPrimary);
2540 primary->AddDaughter(mother);
2541 mother->ConvertAliPIDtoAODPID();
2542 }
2543 else {
2544 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2545 << " track " << imother << " has already been used!" << endl;
2546 }
2547
2548 // Add the kink vertex
2549 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
2550
2551 AliAODVertex * vkink =
2552 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
2553 NULL,
2554 0.,
2555 mother,
2556 AliAODVertex::kKink);
2557 // Add the daughter track
2558
2559 AliAODTrack * daughter = NULL;
2560
2561 if (!usedTrack[idaughter]) {
2562
2563 usedTrack[idaughter] = kTRUE;
2564
2565 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
2566 esdTrack->GetPxPyPz(p);
2567 esdTrack->GetXYZ(pos);
2568 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2569 esdTrack->GetESDpid(pid);
2570
2571 daughter =
2572 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2573 esdTrack->GetLabel(),
2574 p,
2575 kTRUE,
2576 pos,
2577 kFALSE,
2578 covTr,
2579 (Short_t)esdTrack->GetSign(),
2580 esdTrack->GetITSClusterMap(),
2581 pid,
2582 vkink,
2583 kTRUE, // check if this is right
2584 kTRUE, // check if this is right
2585 AliAODTrack::kPrimary);
2586 vkink->AddDaughter(daughter);
2587 daughter->ConvertAliPIDtoAODPID();
2588 }
2589 else {
2590 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
2591 << " track " << idaughter << " has already been used!" << endl;
2592 }
2593
2594
2595 }
2596 }
2597
2598 }
2599
2600 }
2601
2602
2603 // Tracks (primary and orphan)
2604
2605 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
2606
2607
2608 if (usedTrack[nTrack]) continue;
2609
2610 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
2611 esdTrack->GetPxPyPz(p);
2612 esdTrack->GetXYZ(pos);
2613 esdTrack->GetCovarianceXYZPxPyPz(covTr);
2614 esdTrack->GetESDpid(pid);
2615
2616 Float_t impactXY, impactZ;
2617
2618 esdTrack->GetImpactParameters(impactXY,impactZ);
2619
2620 if (impactXY<3) {
2621 // track inside the beam pipe
2622
2623 primary->AddDaughter(aodTrack =
2624 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2625 esdTrack->GetLabel(),
2626 p,
2627 kTRUE,
2628 pos,
2629 kFALSE,
2630 covTr,
2631 (Short_t)esdTrack->GetSign(),
2632 esdTrack->GetITSClusterMap(),
2633 pid,
2634 primary,
2635 kTRUE, // check if this is right
2636 kTRUE, // check if this is right
2637 AliAODTrack::kPrimary)
2638 );
2639 aodTrack->ConvertAliPIDtoAODPID();
2640 }
2641 else {
2642 // outside the beam pipe: orphan track
2643 aodTrack =
2644 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
2645 esdTrack->GetLabel(),
2646 p,
2647 kTRUE,
2648 pos,
2649 kFALSE,
2650 covTr,
2651 (Short_t)esdTrack->GetSign(),
2652 esdTrack->GetITSClusterMap(),
2653 pid,
2654 NULL,
2655 kFALSE, // check if this is right
2656 kFALSE, // check if this is right
2657 AliAODTrack::kOrphan);
2658 aodTrack->ConvertAliPIDtoAODPID();
2659 }
2660 } // end of loop on tracks
2661
2662 // muon tracks
2663 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
2664 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
2665
2666 AliESDMuonTrack *esdMuTrack = esd->GetMuonTrack(nMuTrack);
2667 p[0] = esdMuTrack->Px();
2668 p[1] = esdMuTrack->Py();
2669 p[2] = esdMuTrack->Pz();
2670 pos[0] = primary->GetX();
2671 pos[1] = primary->GetY();
2672 pos[2] = primary->GetZ();
2673
2674 // has to be changed once the muon pid is provided by the ESD
2675 for (Int_t i = 0; i < 10; pid[i++] = 0.); pid[AliAODTrack::kMuon]=1.;
2676
2677 primary->AddDaughter(
2678 new(tracks[jTracks++]) AliAODTrack(0, // no ID provided
2679 0, // no label provided
2680 p,
2681 kTRUE,
2682 pos,
2683 kFALSE,
2684 NULL, // no covariance matrix provided
2685 (Short_t)-99, // no charge provided
2686 0, // no ITSClusterMap
2687 pid,
2688 primary,
2689 kTRUE, // check if this is right
2690 kTRUE, // not used for vertex fit
2691 AliAODTrack::kPrimary)
2692 );
2693 }
2694
2695 // Access to the AOD container of clusters
2696 TClonesArray &clusters = *(aod->GetClusters());
2697 Int_t jClusters=0;
2698
2699 // Calo Clusters
2700 Int_t nClusters = esd->GetNumberOfCaloClusters();
2701
2702 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
2703
2704 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
2705
2706 Int_t id = cluster->GetID();
2707 Int_t label = -1;
2708 Float_t energy = cluster->GetClusterEnergy();
2709 cluster->GetGlobalPosition(posF);
2710 AliAODVertex *prodVertex = primary;
2711 AliAODTrack *primTrack = NULL;
2712 Char_t ttype=AliAODCluster::kUndef;
2713
2714 if (cluster->IsPHOS()) ttype=AliAODCluster::kPHOSNeutral;
2715 else if (cluster->IsEMCAL()) {
2716
2717 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
2718 ttype = AliAODCluster::kEMCALPseudoCluster;
2719 else
2720 ttype = AliAODCluster::kEMCALClusterv1;
2721
2722 }
2723
2724 new(clusters[jClusters++]) AliAODCluster(id,
2725 label,
2726 energy,
2727 pos,
2728 NULL, // no covariance matrix provided
2729 NULL, // no pid for clusters provided
2730 prodVertex,
2731 primTrack,
2732 ttype);
2733
2734 } // end of loop on calo clusters
2735
2736 delete [] usedTrack;
2737 delete [] usedV0;
2738 delete [] usedKink;
2739
2740 // fill the tree for this event
2741 aodTree->Fill();
2742 } // end of event loop
2743
2744 aodTree->GetUserInfo()->Add(aod);
2745
2746 // close ESD file
2747 esdFile->Close();
2748
2749 // write the tree to the specified file
2750 aodFile = aodTree->GetCurrentFile();
2751 aodFile->cd();
2752 aodTree->Write();
2753
2754 return;
2755}
2756
2757
2758void AliReconstruction::WriteAlignmentData(AliESD* esd)
2759{
2760 // Write space-points which are then used in the alignment procedures
2761 // For the moment only ITS, TRD and TPC
2762
2763 // Load TOF clusters
2764 if (fTracker[3]){
2765 fLoader[3]->LoadRecPoints("read");
2766 TTree* tree = fLoader[3]->TreeR();
2767 if (!tree) {
2768 AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
2769 return;
2770 }
2771 fTracker[3]->LoadClusters(tree);
2772 }
2773 Int_t ntracks = esd->GetNumberOfTracks();
2774 for (Int_t itrack = 0; itrack < ntracks; itrack++)
2775 {
2776 AliESDtrack *track = esd->GetTrack(itrack);
2777 Int_t nsp = 0;
2778 Int_t idx[200];
2779 for (Int_t iDet = 3; iDet >= 0; iDet--)
2780 nsp += track->GetNcls(iDet);
2781 if (nsp) {
2782 AliTrackPointArray *sp = new AliTrackPointArray(nsp);
2783 track->SetTrackPointArray(sp);
2784 Int_t isptrack = 0;
2785 for (Int_t iDet = 3; iDet >= 0; iDet--) {
2786 AliTracker *tracker = fTracker[iDet];
2787 if (!tracker) continue;
2788 Int_t nspdet = track->GetNcls(iDet);
2789 if (nspdet <= 0) continue;
2790 track->GetClusters(iDet,idx);
2791 AliTrackPoint p;
2792 Int_t isp = 0;
2793 Int_t isp2 = 0;
2794 while (isp < nspdet) {
2795 Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
2796 const Int_t kNTPCmax = 159;
2797 if (iDet==1 && isp2>kNTPCmax) break; // to be fixed
2798 if (!isvalid) continue;
2799 sp->AddPoint(isptrack,&p); isptrack++; isp++;
2800 }
2801 }
2802 }
2803 }
2804 if (fTracker[3]){
2805 fTracker[3]->UnloadClusters();
2806 fLoader[3]->UnloadRecPoints();
2807 }
2808}
2809
2810//_____________________________________________________________________________
2811void AliReconstruction::FillRawDataErrorLog(Int_t iEvent, AliESD* esd)
2812{
2813 // The method reads the raw-data error log
2814 // accumulated within the rawReader.
2815 // It extracts the raw-data errors related to
2816 // the current event and stores them into
2817 // a TClonesArray inside the esd object.
2818
2819 if (!fRawReader) return;
2820
2821 for(Int_t i = 0; i < fRawReader->GetNumberOfErrorLogs(); i++) {
2822
2823 AliRawDataErrorLog *log = fRawReader->GetErrorLog(i);
2824 if (!log) continue;
2825 if (iEvent != log->GetEventNumber()) continue;
2826
2827 esd->AddRawDataErrorLog(log);
2828 }
2829
2830}