Adding simple example to load default debug streamer
[u/mrichter/AliRoot.git] / TPC / AliTPCTracking.C
index e5b0e2e..c10c0ef 100644 (file)
@@ -6,21 +6,28 @@
 // author: Jiri Chudoba based on code of Jourij Belikov
 // version: 1.0
 // description: 
-//      reconstructs of tracks in TPC inthe following steps:
+//      reconstructs of tracks in TPC in the following steps:
 //         TPC cluster finding
 //         TPC track finding
 // input parameters: 
-//        Int_t nEvents      ... nr of events to process
+//        Int_t nEvents      ... nr of events to process (<0 means all)
 //        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)
+//        const char* fileName ... name of galice file
+//        Bool_t makeClusters ... run the cluster finder or not?
+//        Bool_t makeTracks ... run the track finder or not?
+//        const char* fileNameRaw ... if not NULL, the cluster finder uses
+//                                    the given file as raw data input
 //
 // History:
 //
+//     21.07.2003 ... NewIO
+//
+//     18.03.2003 ... Char_t* replaced by const char*
+//
+//     03.03.2003 ... SetFieldFactor moved to AliTracker class and
+//                    LoadTPCParam moved to AliTPC class
+//                    TString replaced by Char_t*
+//
 //     20.11.2002 ... Changes due to a changed interface of AliTPCtracker. 
 //                    Use Riostream.h instead of iostream.h
 //
 ////////////////////////////////////////////////////////////////////////
 
 #if !defined(__CINT__) || defined(__MAKECINT__)
-#include "Riostream.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"
-
+#include <Riostream.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TBenchmark.h>
+#include "AliRunLoader.h"
+#include "AliTPC.h"
+#include "AliTPCParam.h"
+#include "AliTPCclustererMI.h"
+#include "AliTPCtrackerMI.h"
+#include "AliRawReaderRoot.h"
 #endif
 
-Int_t gDEBUG = 2;
+Bool_t AliTPCTracking(Int_t nEvents = -1, Int_t firstEvent = 0,
+                     const char* fileName = "galice.root",
+                     Bool_t makeClusters = kTRUE,
+                     Bool_t makeTracks = kTRUE,
+                     const char* fileNameRaw = NULL);
 
-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 TPCRefitInward(Int_t nEvents=1, Int_t firstEvent=0,
+                    const char* fileNameClusters="tpc.clusters.root",
+                    const char* fileNameTracks="tpc.tracks.root",
+                    const char* fileNameTRDTracks="trd.tracks.root",
+                    const char* fileNameRefittedTracks="tpc.refitted.tracks.root");
+*/
 
 ////////////////////////////////////////////////////////////////////////
-Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
-                     TString fileNameHits,
-                     TString fileNameDigits,
-                     TString fileNameClusters,
-                     TString fileNameTracks) {
-
-  SetFieldFactor(fileNameHits,kFALSE);
+Bool_t AliTPCTracking(Int_t nEvents, Int_t firstEvent,
+                     const char* fileName,
+                     Bool_t makeClusters,
+                     Bool_t makeTracks,
+                     const char* fileNameRaw) 
+{
+  // get the loaders
+  AliRunLoader* runLoader = AliRunLoader::Open(fileName);
+  if (!runLoader) {
+    cerr << "AliTPCTracking: no run loader found\n";
+    return kFALSE;
+  }
+  AliLoader* tpcLoader = runLoader->GetLoader("TPCLoader");
+  if (!tpcLoader) {
+    cerr << "AliTPCTracking: no TPC loader found\n";
+    return kFALSE;
+  }
 
-// ********** Find TPC clusters *********** //
-  if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
-    cerr<<"Failed to get TPC clusters: !\n";
-    return 1;
-  }      
+  // get the TPC parameters
+  runLoader->CdGAFile();
+  AliTPCParam* param = AliTPC::LoadTPCParam(gFile);
+  if (!param) {
+    cerr << "AliTPCTracking: no TPC parameters found\n";
+    return kFALSE;
+  }
 
-// ********** Find TPC tracks *********** //
-  if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
-    cerr<<"Failed to get TPC tracks !\n";
-    return 2;
+  // create the clusterer object
+  AliTPCclustererMI* clusterer = NULL;
+  if (makeClusters) {
+    clusterer = new AliTPCclustererMI(param);
+    if (!fileNameRaw) tpcLoader->LoadDigits();
+    tpcLoader->LoadRecPoints("recreate");
   }
 
-  return 0;
-}
+  // create the tracker object
+  AliTPCtrackerMI* tracker = NULL;
+  if (makeTracks) {
+//    tracker = new AliTPCtrackerMI(param);
+    if (!makeClusters) tpcLoader->LoadRecPoints();
+    tpcLoader->LoadTracks("recreate");
+  }
 
-////////////////////////////////////////////////////////////////////////
-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;
+  // get the event number range
+  Int_t maxEvent = 0;
+  if (fileNameRaw) {
+    TFile* file = TFile::Open(fileNameRaw);
+    if (file && file->IsOpen()) {
+      TTree* tree = (TTree*) file->Get("T");
+      if (tree) maxEvent = (Int_t) tree->GetEntries();
+    }
+  } else {
+    maxEvent = runLoader->GetNumberOfEvents();
   }
-  if (!fileClusters->IsOpen()) {
-    cerr<<"Cannnot open "<<fileNameClusters<<" !\n"; 
-    return 1;
+  if (nEvents < 0) nEvents = maxEvent - firstEvent;
+  Int_t lastEvent = firstEvent + nEvents;
+  if (lastEvent > maxEvent) lastEvent = maxEvent;
+
+  // loop over the events
+  for (Int_t iEvent = firstEvent; iEvent < lastEvent; iEvent++) {
+
+    runLoader->GetEvent(iEvent);
+
+    // run the cluster finder
+    if (makeClusters) {
+      if (!tpcLoader->TreeR()) tpcLoader->MakeRecPointsContainer();
+      clusterer->SetOutput(tpcLoader->TreeR());
+      if (fileNameRaw) {
+       AliRawReaderRoot rawReader(fileNameRaw, iEvent);
+       clusterer->Digits2Clusters(&rawReader);
+      } else {
+       clusterer->SetInput(tpcLoader->TreeD());
+       clusterer->Digits2Clusters();
+      }
+      tpcLoader->WriteRecPoints("OVERWRITE");
+    }
+
+    // run the track finder
+    if (makeTracks) {
+      tracker = new AliTPCtrackerMI(param);
+      tracker->Clusters2Tracks();
+      delete tracker;
+    }
   }
 
-  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);
+  if (tracker) delete tracker;
+  if (clusterer) delete clusterer;
 
-  return rc;
+  return kTRUE;
 }
-////////////////////////////////////////////////////////////////////////
-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 TPCRefitInward(Int_t nEvents, Int_t firstEvent,
+                    const char* fileNameClusters,
+                    const char* fileNameTracks,
+                    const char* fileNameTRDTracks,
+                    const char* fileNameRefittedTracks)
+{
   Int_t rc = 0;
-  const Char_t *name="TPCFindTracks";
+  const Char_t *name="TPCRefitInward";
   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;
+  TFile *fileClusters = TFile::Open(fileNameClusters);
+  TFile *fileTracks = TFile::Open(fileNameTracks);
+  TFile *fileTRDTracks = TFile::Open(fileNameTRDTracks);
+  TFile *fileRefittedTracks = TFile::Open(fileNameRefittedTracks, "recreate");
 
-}
-////////////////////////////////////////////////////////////////////////
-Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
-                   TFile *fileClusters, TFile * fileTracks,
-                   AliTPCParam* paramTPC) {
-
-  Int_t rc = 0;
-  if (!paramTPC) paramTPC = LoadTPCParam(fileClusters);
+  AliTPCParam* paramTPC = AliTPC::LoadTPCParam(fileClusters);
   if (!paramTPC) return 1;
 
   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
-    if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
+    if (gDEBUG > 2) cout<<"TPCRefitInward: event "<<iEvent<<endl;
     AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
     tracker->SetEventNumber(iEvent);
-    Double_t vertex[3];
-    FindVertex(iEvent,vertex);
-    tracker->SetVertex(vertex);
     fileClusters->cd();
-    rc = tracker->Clusters2Tracks(0,fileTracks);
+    rc = tracker->RefitInward(fileTRDTracks, fileTracks, fileRefittedTracks);
     delete tracker;
+    if (rc) return rc;
   }
-  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;
-  }
+  fileClusters->Close();
+  fileTracks->Close();
+  fileTRDTracks->Close();
+  fileRefittedTracks->Close();
+  delete fileClusters;
+  delete fileTracks;
+  delete fileTRDTracks;
+  delete fileRefittedTracks;
+  if (gDEBUG>1) gBenchmark->Show(name);
   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;
-} 
-////////////////////////////////////////////////////////////////////////
+*/