]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliReconstruction.cxx
Correct overloading of virtual functions in the derived classes (icc)
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
index fd5096f4a49251d05e0a997c6a8dcfe6e7e7a235..ff8731f5f873f95529b3ae15f6f97354278be48b 100644 (file)
 //                                                                           //
 //   rec.SetRunTracking("...");                                              //
 //                                                                           //
+// Uniform/nonuniform field tracking switches (default: uniform field)       //
+//                                                                           //
+//   rec.SetUniformFieldTracking();  ( rec.SetNonuniformFieldTracking(); )   //
+//                                                                           //
 // The filling of additional ESD information can be steered by               //
 //                                                                           //
 //   rec.SetFillESD("...");                                                  //
 #include "AliRawReaderFile.h"
 #include "AliRawReaderDate.h"
 #include "AliRawReaderRoot.h"
-#include "AliTracker.h"
 #include "AliESD.h"
 #include "AliESDVertex.h"
+#include "AliTracker.h"
 #include "AliVertexer.h"
 #include "AliHeader.h"
 #include "AliGenEventHeader.h"
 #include "AliPID.h"
 #include "AliESDpid.h"
-#include "AliMagF.h"
+//#include "AliMagF.h"
+
+
+
+#include "AliRunTag.h"
+//#include "AliLHCTag.h"
+#include "AliDetectorTag.h"
+#include "AliEventTag.h"
+
+
 
 ClassImp(AliReconstruction)
 
@@ -138,6 +151,7 @@ AliReconstruction::AliReconstruction(const char* gAliceFilename,
   TNamed(name, title),
 
   fRunLocalReconstruction("ALL"),
+  fUniformField(kTRUE),
   fRunVertexFinder(kTRUE),
   fRunHLTTracking(kFALSE),
   fRunTracking("ALL"),
@@ -170,6 +184,7 @@ AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
   TNamed(rec),
 
   fRunLocalReconstruction(rec.fRunLocalReconstruction),
+  fUniformField(rec.fUniformField),
   fRunVertexFinder(rec.fRunVertexFinder),
   fRunHLTTracking(rec.fRunHLTTracking),
   fRunTracking(rec.fRunTracking),
@@ -407,7 +422,8 @@ Bool_t AliReconstruction::Run(const char* input,
     // write HLT ESD
     hlttree->Fill();
 
-    if (fCheckPointLevel > 0) WriteESD(esd, "final");
+    if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
     delete esd; delete hltesd;
     esd = NULL; hltesd = NULL;
   }
@@ -415,6 +431,10 @@ Bool_t AliReconstruction::Run(const char* input,
   file->cd();
   tree->Write();
   hlttree->Write();
+
+  // Create tags for the events in the ESD tree (the ESD tree is always present)
+  // In case of empty events the tags will contain dummy values
+  CreateTag(file);
   CleanUp(file, fileOld);
 
   return kTRUE;
@@ -445,8 +465,9 @@ Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
     } else {
       reconstructor->Reconstruct(fRunLoader);
     }
-    AliInfo(Form("execution time for %s:", fgkDetectorName[iDet]));
-    ToAliInfo(stopwatchDet.Print());
+    AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
+                fgkDetectorName[iDet],
+                stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
   }
 
   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
@@ -455,8 +476,8 @@ Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
     if (fStopOnError) return kFALSE;
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -489,8 +510,9 @@ Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
       reconstructor->ConvertDigits(fRawReader, digitsTree);
       loader->WriteDigits("OVERWRITE");
       loader->UnloadDigits();
-      AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet]));
-      ToAliDebug(1, stopwatchDet.Print());
+      AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
+                  fgkDetectorName[iDet],
+                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
     }
 
     // local reconstruction
@@ -517,8 +539,9 @@ Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
     }
     loader->WriteRecPoints("OVERWRITE");
     loader->UnloadRecPoints();
-    AliDebug(1, Form("execution time for %s:", fgkDetectorName[iDet]));
-    ToAliDebug(1, stopwatchDet.Print());
+    AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
+                   fgkDetectorName[iDet],
+                   stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
   }
 
   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
@@ -526,9 +549,9 @@ Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
                   detStr.Data()));
     if (fStopOnError) return kFALSE;
   }
-
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -581,8 +604,8 @@ Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
   }  
   delete vertex;
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -615,6 +638,7 @@ Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
     if (!tracker) {
       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
       if (fStopOnError) return kFALSE;
+      continue;
     }
     Double_t vtxPos[3];
     Double_t vtxErr[3]={0.005,0.005,0.010};
@@ -640,8 +664,8 @@ Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
     delete tracker;
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -694,13 +718,8 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     // load clusters
     if (iDet > 1) {     // all except ITS, TPC
       TTree* tree = NULL;
-      if (iDet == 3) {   // TOF
-       fLoader[iDet]->LoadDigits("read");
-       tree = fLoader[iDet]->TreeD();
-      } else {
-       fLoader[iDet]->LoadRecPoints("read");
-       tree = fLoader[iDet]->TreeR();
-      }
+      fLoader[iDet]->LoadRecPoints("read");
+      tree = fLoader[iDet]->TreeR();
       if (!tree) {
        AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
        return kFALSE;
@@ -720,11 +739,7 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     // unload clusters
     if (iDet > 2) {     // all except ITS, TPC, TRD
       fTracker[iDet]->UnloadClusters();
-      if (iDet == 3) {   // TOF
-       fLoader[iDet]->UnloadDigits();
-      } else {
-       fLoader[iDet]->UnloadRecPoints();
-      }
+      fLoader[iDet]->UnloadRecPoints();
     }
   }
 
@@ -747,8 +762,8 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     fLoader[iDet]->UnloadRecPoints();
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -815,8 +830,8 @@ Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
     if (fStopOnError) return kFALSE;
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -882,7 +897,12 @@ Bool_t AliReconstruction::InitRunLoader()
     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
       if (fRunLoader->LoadgAlice() == 0) {
        gAlice = fRunLoader->GetAliRun();
-       AliTracker::SetFieldMap(gAlice->Field());
+       AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
+       AliExternalTrackParam::SetFieldMap(gAlice->Field());
+       if(fUniformField)
+         AliExternalTrackParam::SetUniformFieldTracking();
+       else
+         AliExternalTrackParam::SetNonuniformFieldTracking();
       }
     }
     if (!gAlice && !fRawReader) {
@@ -1107,6 +1127,7 @@ Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
          esd->GetRunNumber(), esd->GetEventNumber(), recStep);
   if (gSystem->AccessPathName(fileName)) return kFALSE;
 
+  AliInfo(Form("reading ESD from file %s", fileName));
   AliDebug(1, Form("reading ESD from file %s", fileName));
   TFile* file = TFile::Open(fileName);
   if (!file || !file->IsOpen()) {
@@ -1143,3 +1164,251 @@ void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
   }
   delete file;
 }
+
+
+
+
+//_____________________________________________________________________________
+void AliReconstruction::CreateTag(TFile* file)
+{
+  // Creates the tags for all the events in a given ESD file
+  Int_t ntrack;
+  Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
+  Int_t nPos, nNeg, nNeutr;
+  Int_t nK0s, nNeutrons, nPi0s, nGamas;
+  Int_t nCh1GeV, nCh3GeV, nCh10GeV;
+  Int_t nMu1GeV, nMu3GeV, nMu10GeV;
+  Int_t nEl1GeV, nEl3GeV, nEl10GeV;
+  Float_t maxPt = .0, meanPt = .0, totalP = .0;
+
+  AliRunTag *tag = new AliRunTag();
+  AliDetectorTag *detTag = new AliDetectorTag();
+  AliEventTag *evTag = new AliEventTag();
+  TTree ttag("T","A Tree with event tags");
+  TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
+  btag->SetCompressionLevel(9);
+
+  AliInfo(Form("Creating the tags......."));   
+  
+  if (!file || !file->IsOpen()) {
+    AliError(Form("opening failed"));
+    delete file;
+    return ;
+  }
+
+  TTree *t = (TTree*) file->Get("esdTree");
+  TBranch * b = t->GetBranch("ESD");
+  AliESD *esd = 0;
+  b->SetAddress(&esd);
+
+  tag->SetRunId(esd->GetRunNumber());
+
+  Int_t firstEvent = 0,lastEvent = 0;
+  Int_t iNumberOfEvents = b->GetEntries();
+  for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++)
+    {
+      ntrack = 0;
+      nPos = 0;
+      nNeg = 0;
+      nNeutr =0;
+      nK0s = 0;
+      nNeutrons = 0;
+      nPi0s = 0;
+      nGamas = 0;
+      nProtons = 0;
+      nKaons = 0;
+      nPions = 0;
+      nMuons = 0;
+      nElectrons = 0;    
+      nCh1GeV = 0;
+      nCh3GeV = 0;
+      nCh10GeV = 0;
+      nMu1GeV = 0;
+      nMu3GeV = 0;
+      nMu10GeV = 0;
+      nEl1GeV = 0;
+      nEl3GeV = 0;
+      nEl10GeV = 0;
+      maxPt = .0;
+      meanPt = .0;
+      totalP = .0;
+
+      b->GetEntry(iEventNumber);
+      const AliESDVertex * vertexIn = esd->GetVertex();
+
+      for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++)
+       {
+         AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
+         UInt_t status = esdTrack->GetStatus();
+         
+         //select only tracks with ITS refit
+         if ((status&AliESDtrack::kITSrefit)==0) continue;
+         
+         //select only tracks with TPC refit-->remove extremely high Pt tracks
+         if ((status&AliESDtrack::kTPCrefit)==0) continue;
+         
+         //select only tracks with the "combined PID"
+         if ((status&AliESDtrack::kESDpid)==0) continue;
+                 Double_t p[3];
+         esdTrack->GetPxPyPz(p);
+         Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
+         Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
+         totalP += momentum;
+         meanPt += fPt;
+         if(fPt > maxPt)
+           maxPt = fPt;
+         
+         if(esdTrack->GetSign() > 0)
+           {
+             nPos++;
+             if(fPt > 1.0)
+               nCh1GeV++;
+             if(fPt > 3.0)
+               nCh3GeV++;
+             if(fPt > 10.0)
+               nCh10GeV++;
+           }
+         if(esdTrack->GetSign() < 0)
+           {
+             nNeg++;
+             if(fPt > 1.0)
+               nCh1GeV++;
+             if(fPt > 3.0)
+               nCh3GeV++;
+             if(fPt > 10.0)
+               nCh10GeV++;
+           }
+         if(esdTrack->GetSign() == 0)
+           nNeutr++;
+         
+         //PID
+         Double_t prob[10];
+         esdTrack->GetESDpid(prob);
+                   
+         //K0s
+         if ((prob[8]>prob[7])&&(prob[8]>prob[6])&&(prob[8]>prob[5])&&(prob[8]>prob[4])&&(prob[8]>prob[3])&&(prob[8]>prob[2])&&(prob[8]>prob[1])&&(prob[8]>prob[0]))
+           nK0s++;
+         //neutrons
+         if ((prob[7]>prob[8])&&(prob[7]>prob[6])&&(prob[7]>prob[5])&&(prob[7]>prob[4])&&(prob[7]>prob[3])&&(prob[7]>prob[2])&&(prob[7]>prob[1])&&(prob[7]>prob[0]))
+           nNeutrons++; 
+         //pi0s
+         if ((prob[6]>prob[8])&&(prob[6]>prob[7])&&(prob[6]>prob[5])&&(prob[6]>prob[4])&&(prob[6]>prob[3])&&(prob[6]>prob[2])&&(prob[6]>prob[1])&&(prob[6]>prob[0]))
+           nPi0s++;
+         //gamas
+         if ((prob[5]>prob[8])&&(prob[5]>prob[7])&&(prob[5]>prob[6])&&(prob[5]>prob[4])&&(prob[5]>prob[3])&&(prob[5]>prob[2])&&(prob[5]>prob[1])&&(prob[5]>prob[0]))
+           nGamas++;
+         //protons
+         if ((prob[4]>prob[8])&&(prob[4]>prob[7])&&(prob[4]>prob[6])&&(prob[4]>prob[5])&&(prob[4]>prob[3])&&(prob[4]>prob[2])&&(prob[4]>prob[1])&&(prob[4]>prob[0]))
+           nProtons++;
+         //kaons
+         if ((prob[3]>prob[8])&&(prob[3]>prob[7])&&(prob[3]>prob[6])&&(prob[3]>prob[5])&&(prob[3]>prob[4])&&(prob[3]>prob[2])&&(prob[3]>prob[1])&&(prob[3]>prob[0]))
+           nKaons++;
+         //kaons
+         if ((prob[2]>prob[8])&&(prob[2]>prob[7])&&(prob[2]>prob[6])&&(prob[2]>prob[5])&&(prob[2]>prob[4])&&(prob[2]>prob[3])&&(prob[2]>prob[1])&&(prob[2]>prob[0]))
+           nPions++; 
+         //muons
+         if ((prob[1]>prob[8])&&(prob[1]>prob[7])&&(prob[1]>prob[6])&&(prob[1]>prob[5])&&(prob[1]>prob[4])&&(prob[1]>prob[3])&&(prob[1]>prob[2])&&(prob[1]>prob[0]))
+           {
+             nMuons++;
+             if(fPt > 1.0)
+               nMu1GeV++;
+             if(fPt > 3.0)
+               nMu3GeV++;
+             if(fPt > 10.0)
+               nMu10GeV++;
+           }
+         //electrons
+         if ((prob[0]>prob[8])&&(prob[0]>prob[7])&&(prob[0]>prob[6])&&(prob[0]>prob[5])&&(prob[0]>prob[4])&&(prob[0]>prob[3])&&(prob[0]>prob[2])&&(prob[0]>prob[1]))
+           {
+             nElectrons++;
+             if(fPt > 1.0)
+               nEl1GeV++;
+             if(fPt > 3.0)
+               nEl3GeV++;
+             if(fPt > 10.0)
+               nEl10GeV++;
+           }
+         
+         
+         
+         ntrack++;
+       }//track loop
+      // Fill the event tags 
+      if(ntrack != 0)
+       meanPt = meanPt/ntrack;
+      
+      evTag->SetEventId(iEventNumber+1);
+      evTag->SetVertexX(vertexIn->GetXv());
+      evTag->SetVertexY(vertexIn->GetYv());
+      evTag->SetVertexZ(vertexIn->GetZv());
+      
+      evTag->SetT0VertexZ(esd->GetT0zVertex());
+      
+      evTag->SetTrigger(esd->GetTrigger());
+      
+      evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
+      evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
+      evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
+      evTag->SetNumOfParticipants(esd->GetZDCParticipants());
+      
+      
+      evTag->SetNumOfTracks(esd->GetNumberOfTracks());
+      evTag->SetNumOfPosTracks(nPos);
+      evTag->SetNumOfNegTracks(nNeg);
+      evTag->SetNumOfNeutrTracks(nNeutr);
+      
+      evTag->SetNumOfV0s(esd->GetNumberOfV0s());
+      evTag->SetNumOfCascades(esd->GetNumberOfCascades());
+      evTag->SetNumOfKinks(esd->GetNumberOfKinks());
+      evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
+      
+      evTag->SetNumOfProtons(nProtons);
+      evTag->SetNumOfKaons(nKaons);
+      evTag->SetNumOfPions(nPions);
+      evTag->SetNumOfMuons(nMuons);
+      evTag->SetNumOfElectrons(nElectrons);
+      evTag->SetNumOfPhotons(nGamas);
+      evTag->SetNumOfPi0s(nPi0s);
+      evTag->SetNumOfNeutrons(nNeutrons);
+      evTag->SetNumOfKaon0s(nK0s);
+      
+      evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
+      evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
+      evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
+      evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
+      evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
+      evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
+      evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
+      evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
+      evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
+      
+      evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
+      evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
+      
+      evTag->SetTotalMomentum(totalP);
+      evTag->SetMeanPt(meanPt);
+      evTag->SetMaxPt(maxPt);
+  
+      tag->AddEventTag(*evTag);
+    }
+  lastEvent = iNumberOfEvents;
+       
+  ttag.Fill();
+  tag->Clear();
+
+  char fileName[256];
+  sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
+         tag->GetRunId(),firstEvent,lastEvent );
+  AliInfo(Form("writing tags to file %s", fileName));
+  AliDebug(1, Form("writing tags to file %s", fileName));
+  TFile* ftag = TFile::Open(fileName, "recreate");
+  ftag->cd();
+  ttag.Write();
+  ftag->Close();
+  file->cd();
+  delete tag;
+  delete detTag;
+  delete evTag;
+}
+