General tracking macro by J. Chudoba.

diff --git a/TPC/AliTPCTracking.C b/TPC/AliTPCTracking.C
new file mode 100644 (file)
index 0000000..6395d28
--- /dev/null
+++ b/TPC/AliTPCTracking.C
@@ -0,0 +1,296 @@
+////////////////////////////////////////////////////////////////////////
+//
+// AliTPCTracking.C 
+//
+// date: 22.08.2002
+// author: Jiri Chudoba based on code of Jourij Belikov
+// version: 1.0
+// description: 
+//      reconstructs of tracks in TPC inthe following steps:
+//         TPC cluster finding
+//         TPC track finding
+// input parameters: 
+//        Int_t nEvents      ... nr of events to process
+//        Int_t firstEventNr ... first event number (starts from 0)
+//        TString fileNameHits ... name of file with hits
+//        TString fileNameDigits .. name of file with TPC digits
+//        TString fileNameClusters .. name of file with TPC clusters (output)
+//        TString fileNameTracks .. name of file with TPC tracks (output)
+//
+//        default file names correspond to pp production (2002-04)
+//
+// History:
+//
+//     22.08.2002 ... first version
+//
+////////////////////////////////////////////////////////////////////////
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "iostream.h"
+#include "TTree.h"
+#include "TSystem.h"
+#include "TArrayF.h"
+#include "TPC/alles.h"
+#include "TPC/AliTPCtracker.h"
+#include "TPC/AliTPCclusterer.h"
+#include "STEER/AliRun.h"
+#include "STEER/AliHeader.h"
+#include "STEER/AliGenEventHeader.h"
+#include "STEER/AliMagF.h"
+
+#endif
+
+Int_t gDEBUG = 2;
+
+Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
+                     TString fileNameDigits="rfio:galiceSDR.root", 
+                     TString fileNameClusters="tpc.clusters.root");
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TFile* fileDigits, TFile* fileClusters, 
+                     AliTPCParam* paramTPC=0);
+Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
+                   TString fileNameClusters="tpc.clusters.root",
+                   TString fileNameTracks="tpc.tracks.root");
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TFile* fileClusters, TFile* fileTracks,
+                   AliTPCParam* paramTPC=0);
+
+AliTPCParam* LoadTPCParam(TFile *file);
+void FindVertex(Int_t iEvent, Double_t *vertex);
+void PrintVertex(TArrayF &primaryVertex);
+Int_t SetFieldFactor(TString fileName, Bool_t closeFile = kTRUE);
+Int_t SetFieldFactor(TFile* file, Bool_t deletegAlice = kTRUE);
+Int_t SetFieldFactor();
+
+Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
+                    TString fileNameHits="rfio:galice.root",
+                    TString fileNameDigits="rfio:galiceSDR.root",
+                    TString fileNameClusters="tpc.clusters.root",
+                    TString fileNameTracks="tpc.tracks.root");
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
+                     TString fileNameHits,
+                     TString fileNameDigits,
+                     TString fileNameClusters,
+                     TString fileNameTracks) {
+
+  SetFieldFactor(fileNameHits,kFALSE);
+
+// ********** Find TPC clusters *********** //
+  if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
+    cerr<<"Failed to get TPC clusters: !\n";
+    return 1;
+  }      
+
+// ********** Find TPC tracks *********** //
+  if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
+    cerr<<"Failed to get TPC tracks !\n";
+    return 2;
+  }
+
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TString fileNameDigits, TString fileNameClusters) {
+  
+  Int_t rc;
+  const Char_t *name="TPCFindClusters";
+  if (gDEBUG>1) cout<<name<<" starts...\n";
+  if (gDEBUG>1) gBenchmark->Start(name);
+  TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
+  TFile *fileDigits = TFile::Open(fileNameDigits);
+  if (!fileDigits->IsOpen()) {
+    cerr<<"Cannnot open "<<fileNameDigits<<" !\n"; 
+    return 1;
+  }
+  if (!fileClusters->IsOpen()) {
+    cerr<<"Cannnot open "<<fileNameClusters<<" !\n"; 
+    return 1;
+  }
+
+  rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
+
+  fileDigits->Close();
+  fileClusters->Close();
+  delete fileDigits;
+  delete fileClusters;
+  if (gDEBUG>1) gBenchmark->Stop(name);
+  if (gDEBUG>1) gBenchmark->Show(name);
+
+  return rc;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TFile* fileDigits, TFile* fileClusters, 
+                     AliTPCParam* paramTPC) {
+
+  fileDigits->cd();
+  if (!paramTPC) paramTPC = LoadTPCParam(fileDigits);
+  if (!paramTPC) return 1;
+
+  for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
+    if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
+    AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
+  }
+  return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TString fileNameClusters, TString fileNameTracks) {
+
+  Int_t rc = 0;
+  const Char_t *name="TPCFindTracks";
+  if (gDEBUG>1) cout<<name<<" starts"<<endl;
+  if (gDEBUG>1) gBenchmark->Start(name);
+  TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
+  TFile *fileClusters =TFile::Open(fileNameClusters);
+
+  rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
+
+  fileClusters->Close();
+  fileTracks->Close();
+  delete fileClusters;
+  delete fileTracks;
+  if (gDEBUG>1) gBenchmark->Stop(name);
+  if (gDEBUG>1) gBenchmark->Show(name);
+  return rc;
+
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TFile *fileClusters, TFile * fileTracks,
+                   AliTPCParam* paramTPC) {
+
+  Int_t rc = 0;
+  if (!paramTPC) paramTPC = LoadTPCParam(fileClusters);
+  if (!paramTPC) return 1;
+
+  for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
+    if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
+    AliTPCtracker *tracker = new AliTPCtracker(paramTPC,iEvent);
+    Double_t vertex[3];
+    FindVertex(iEvent,vertex);
+    tracker->SetVertex(vertex);
+    fileClusters->cd();
+    rc = tracker->Clusters2Tracks(0,fileTracks);
+    delete tracker;
+  }
+  return rc;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t SetFieldFactor() {
+
+   AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
+   if (gDEBUG > 2) cout<<"Magnetic field in kGauss: "<<gAlice->Field()->SolenoidField()<<endl;
+   return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t SetFieldFactor(TFile *file, Bool_t deletegAlice) {
+
+  if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
+    cerr<<"gAlice has not been found in file "<<file->GetName();
+    return 1;
+  }   
+  Int_t rc = SetFieldFactor();
+  if (deletegAlice) {
+    delete gAlice;  
+    gAlice = 0;
+  }
+  return rc;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t SetFieldFactor(TString fileName, Bool_t closeFile) {
+
+   TFile *file=TFile::Open(fileName);
+   if (!file->IsOpen()) {cerr<<"Cannnot open "<<fileName<<" !\n"; return 1;}
+   Int_t rc = SetFieldFactor(file, closeFile) ;
+   if (closeFile) file->Close();
+   return rc;
+}
+////////////////////////////////////////////////////////////////////////
+AliTPCParam* LoadTPCParam(TFile *file) {
+
+  char paramName[50];
+  sprintf(paramName,"75x40_100x60_150x60");
+  AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
+  if (paramTPC) {
+    if (gDEBUG > 1) cout<<"TPC parameters "<<paramName<<" found."<<endl;
+  } else {
+    cerr<<"TPC parameters not found. Create new (they may be incorrect)."
+       <<endl;    
+    paramTPC = new AliTPCParamSR;
+  }
+  return paramTPC;
+
+// the older version of parameters can be accessed with this code.
+// In some cases, we have old parameters saved in the file but 
+// digits were created with new parameters, it can be distinguish 
+// by the name of TPC TreeD. The code here is just for the case 
+// we would need to compare with old data, uncomment it if needed.
+//
+//  char paramName[50];
+//  sprintf(paramName,"75x40_100x60");
+//  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
+//  if (paramTPC) {
+//    cout<<"TPC parameters "<<paramName<<" found."<<endl;
+//  } else {
+//    sprintf(paramName,"75x40_100x60_150x60");
+//    paramTPC=(AliTPCParam*)in->Get(paramName);
+//    if (paramTPC) {
+//     cout<<"TPC parameters "<<paramName<<" found."<<endl;
+//    } else {
+//     cerr<<"TPC parameters not found. Create new (they may be incorrect)."
+//         <<endl;    
+//     paramTPC = new AliTPCParamSR;
+//    }
+//  }
+//  return paramTPC;
+
+}
+
+////////////////////////////////////////////////////////////////////////
+void FindVertex(Int_t eventNr, Double_t *vertex) {
+
+  vertex[0] = vertex[1] = vertex[2] = 0.;
+  if (!gAlice) {
+    cerr<<"gAlice was not found! Using vertex position (0,0,0).\n";
+    return;
+  }
+  
+  gAlice->GetEvent(eventNr);
+  AliHeader *header = gAlice->GetHeader();
+  if (!header) {
+    cerr<<"header was not found!\n";
+    return;
+  } 
+  AliGenEventHeader* genEventHeader = header->GenEventHeader();
+  if (!genEventHeader) {
+    cerr<<"AliGenEventHeader was not found!\n";
+    return;
+  } 
+
+  TArrayF primaryVertex(3);
+  genEventHeader->PrimaryVertex(primaryVertex);
+  PrintVertex(primaryVertex);
+  vertex[0] = static_cast<Double_t>(primaryVertex[0]);
+  vertex[1] = static_cast<Double_t>(primaryVertex[1]);
+  vertex[2] = static_cast<Double_t>(primaryVertex[2]);
+//  delete header;
+  delete genEventHeader;
+  return;
+     
+}
+////////////////////////////////////////////////////////////////////////
+void PrintVertex(TArrayF &primaryVertex) 
+{
+  cout <<"Vertex: "
+       <<primaryVertex[0]<<" "
+       <<primaryVertex[1]<<" "
+       <<primaryVertex[2]<<" "<<endl;
+  exit;
+} 
+////////////////////////////////////////////////////////////////////////