TPC/AliTPCTracking.C

   1 ////////////////////////////////////////////////////////////////////////
   2 //
   3 // AliTPCTracking.C
   4 //
   5 // date: 22.08.2002
   6 // author: Jiri Chudoba based on code of Jourij Belikov
   7 // version: 1.0
   8 // description:
   9 //      reconstructs of tracks in TPC inthe following steps:
  10 //         TPC cluster finding
  11 //         TPC track finding
  12 // input parameters:
  13 //        Int_t nEvents      ... nr of events to process
  14 //        Int_t firstEventNr ... first event number (starts from 0)
  15 //        TString fileNameHits ... name of file with hits
  16 //        TString fileNameDigits .. name of file with TPC digits
  17 //        TString fileNameClusters .. name of file with TPC clusters (output)
  18 //        TString fileNameTracks .. name of file with TPC tracks (output)
  19 //
  20 //        default file names correspond to pp production (2002-04)
  21 //
  22 // History:
  23 //
  24 //     22.08.2002 ... first version
  25 //
  26 ////////////////////////////////////////////////////////////////////////
  27
  28 #if !defined(__CINT__) || defined(__MAKECINT__)
  29 #include "iostream.h"
  30 #include "TTree.h"
  31 #include "TSystem.h"
  32 #include "TArrayF.h"
  33 #include "TPC/alles.h"
  34 #include "TPC/AliTPCtracker.h"
  35 #include "TPC/AliTPCclusterer.h"
  36 #include "STEER/AliRun.h"
  37 #include "STEER/AliHeader.h"
  38 #include "STEER/AliGenEventHeader.h"
  39 #include "STEER/AliMagF.h"
  40
  41 #endif
  42
  43 Int_t gDEBUG = 2;
  44
  45 Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
  46                       TString fileNameDigits="rfio:galiceSDR.root",
  47                       TString fileNameClusters="tpc.clusters.root");
  48 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
  49                       TFile* fileDigits, TFile* fileClusters,
  50                       AliTPCParam* paramTPC=0);
  51 Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
  52                     TString fileNameClusters="tpc.clusters.root",
  53                     TString fileNameTracks="tpc.tracks.root");
  54 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
  55                     TFile* fileClusters, TFile* fileTracks,
  56                     AliTPCParam* paramTPC=0);
  57
  58 AliTPCParam* LoadTPCParam(TFile *file);
  59 void FindVertex(Int_t iEvent, Double_t *vertex);
  60 void PrintVertex(TArrayF &primaryVertex);
  61 Int_t SetFieldFactor(TString fileName, Bool_t closeFile = kTRUE);
  62 Int_t SetFieldFactor(TFile* file, Bool_t deletegAlice = kTRUE);
  63 Int_t SetFieldFactor();
  64
  65 Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
  66                      TString fileNameHits="rfio:galice.root",
  67                      TString fileNameDigits="rfio:galiceSDR.root",
  68                      TString fileNameClusters="tpc.clusters.root",
  69                      TString fileNameTracks="tpc.tracks.root");
  70
  71 ////////////////////////////////////////////////////////////////////////
  72 Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
  73                       TString fileNameHits,
  74                       TString fileNameDigits,
  75                       TString fileNameClusters,
  76                       TString fileNameTracks) {
  77
  78   SetFieldFactor(fileNameHits,kFALSE);
  79
  80 // ********** Find TPC clusters *********** //
  81   if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
  82     cerr<<"Failed to get TPC clusters: !\n";
  83     return 1;
  84   }
  85
  86 // ********** Find TPC tracks *********** //
  87   if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
  88     cerr<<"Failed to get TPC tracks !\n";
  89     return 2;
  90   }
  91
  92   return 0;
  93 }
  94
  95 ////////////////////////////////////////////////////////////////////////
  96 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
  97                       TString fileNameDigits, TString fileNameClusters) {
  98
  99   Int_t rc;
 100   const Char_t *name="TPCFindClusters";
 101   if (gDEBUG>1) cout<<name<<" starts...\n";
 102   if (gDEBUG>1) gBenchmark->Start(name);
 103   TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
 104   TFile *fileDigits = TFile::Open(fileNameDigits);
 105   if (!fileDigits->IsOpen()) {
 106     cerr<<"Cannnot open "<<fileNameDigits<<" !\n";
 107     return 1;
 108   }
 109   if (!fileClusters->IsOpen()) {
 110     cerr<<"Cannnot open "<<fileNameClusters<<" !\n";
 111     return 1;
 112   }
 113
 114   rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
 115
 116   fileDigits->Close();
 117   fileClusters->Close();
 118   delete fileDigits;
 119   delete fileClusters;
 120   if (gDEBUG>1) gBenchmark->Stop(name);
 121   if (gDEBUG>1) gBenchmark->Show(name);
 122
 123   return rc;
 124 }
 125 ////////////////////////////////////////////////////////////////////////
 126 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
 127                       TFile* fileDigits, TFile* fileClusters,
 128                       AliTPCParam* paramTPC) {
 129
 130   fileDigits->cd();
 131   if (!paramTPC) paramTPC = LoadTPCParam(fileDigits);
 132   if (!paramTPC) return 1;
 133
 134   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
 135     if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
 136     AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
 137   }
 138   return 0;
 139 }
 140 ////////////////////////////////////////////////////////////////////////
 141 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
 142                     TString fileNameClusters, TString fileNameTracks) {
 143
 144   Int_t rc = 0;
 145   const Char_t *name="TPCFindTracks";
 146   if (gDEBUG>1) cout<<name<<" starts"<<endl;
 147   if (gDEBUG>1) gBenchmark->Start(name);
 148   TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
 149   TFile *fileClusters =TFile::Open(fileNameClusters);
 150
 151   rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
 152
 153   fileClusters->Close();
 154   fileTracks->Close();
 155   delete fileClusters;
 156   delete fileTracks;
 157   if (gDEBUG>1) gBenchmark->Stop(name);
 158   if (gDEBUG>1) gBenchmark->Show(name);
 159   return rc;
 160
 161 }
 162 ////////////////////////////////////////////////////////////////////////
 163 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
 164                     TFile *fileClusters, TFile * fileTracks,
 165                     AliTPCParam* paramTPC) {
 166
 167   Int_t rc = 0;
 168   if (!paramTPC) paramTPC = LoadTPCParam(fileClusters);
 169   if (!paramTPC) return 1;
 170
 171   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
 172     if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
 173     AliTPCtracker *tracker = new AliTPCtracker(paramTPC,iEvent);
 174     Double_t vertex[3];
 175     FindVertex(iEvent,vertex);
 176     tracker->SetVertex(vertex);
 177     fileClusters->cd();
 178     rc = tracker->Clusters2Tracks(0,fileTracks);
 179     delete tracker;
 180   }
 181   return rc;
 182 }
 183
 184 ////////////////////////////////////////////////////////////////////////
 185 Int_t SetFieldFactor() {
 186
 187    AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
 188    if (gDEBUG > 2) cout<<"Magnetic field in kGauss: "<<gAlice->Field()->SolenoidField()<<endl;
 189    return 0;
 190 }
 191 ////////////////////////////////////////////////////////////////////////
 192 Int_t SetFieldFactor(TFile *file, Bool_t deletegAlice) {
 193
 194   if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
 195     cerr<<"gAlice has not been found in file "<<file->GetName();
 196     return 1;
 197   }
 198   Int_t rc = SetFieldFactor();
 199   if (deletegAlice) {
 200     delete gAlice;
 201     gAlice = 0;
 202   }
 203   return rc;
 204 }
 205 ////////////////////////////////////////////////////////////////////////
 206 Int_t SetFieldFactor(TString fileName, Bool_t closeFile) {
 207
 208    TFile *file=TFile::Open(fileName);
 209    if (!file->IsOpen()) {cerr<<"Cannnot open "<<fileName<<" !\n"; return 1;}
 210    Int_t rc = SetFieldFactor(file, closeFile) ;
 211    if (closeFile) file->Close();
 212    return rc;
 213 }
 214 ////////////////////////////////////////////////////////////////////////
 215 AliTPCParam* LoadTPCParam(TFile *file) {
 216
 217   char paramName[50];
 218   sprintf(paramName,"75x40_100x60_150x60");
 219   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
 220   if (paramTPC) {
 221     if (gDEBUG > 1) cout<<"TPC parameters "<<paramName<<" found."<<endl;
 222   } else {
 223     cerr<<"TPC parameters not found. Create new (they may be incorrect)."
 224         <<endl;
 225     paramTPC = new AliTPCParamSR;
 226   }
 227   return paramTPC;
 228
 229 // the older version of parameters can be accessed with this code.
 230 // In some cases, we have old parameters saved in the file but
 231 // digits were created with new parameters, it can be distinguish
 232 // by the name of TPC TreeD. The code here is just for the case
 233 // we would need to compare with old data, uncomment it if needed.
 234 //
 235 //  char paramName[50];
 236 //  sprintf(paramName,"75x40_100x60");
 237 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
 238 //  if (paramTPC) {
 239 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
 240 //  } else {
 241 //    sprintf(paramName,"75x40_100x60_150x60");
 242 //    paramTPC=(AliTPCParam*)in->Get(paramName);
 243 //    if (paramTPC) {
 244 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
 245 //    } else {
 246 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
 247 //          <<endl;
 248 //      paramTPC = new AliTPCParamSR;
 249 //    }
 250 //  }
 251 //  return paramTPC;
 252
 253 }
 254
 255 ////////////////////////////////////////////////////////////////////////
 256 void FindVertex(Int_t eventNr, Double_t *vertex) {
 257
 258   vertex[0] = vertex[1] = vertex[2] = 0.;
 259   if (!gAlice) {
 260     cerr<<"gAlice was not found! Using vertex position (0,0,0).\n";
 261     return;
 262   }
 263
 264   gAlice->GetEvent(eventNr);
 265   AliHeader *header = gAlice->GetHeader();
 266   if (!header) {
 267     cerr<<"header was not found!\n";
 268     return;
 269   }
 270   AliGenEventHeader* genEventHeader = header->GenEventHeader();
 271   if (!genEventHeader) {
 272     cerr<<"AliGenEventHeader was not found!\n";
 273     return;
 274   }
 275
 276   TArrayF primaryVertex(3);
 277   genEventHeader->PrimaryVertex(primaryVertex);
 278   PrintVertex(primaryVertex);
 279   vertex[0] = static_cast<Double_t>(primaryVertex[0]);
 280   vertex[1] = static_cast<Double_t>(primaryVertex[1]);
 281   vertex[2] = static_cast<Double_t>(primaryVertex[2]);
 282 //  delete header;
 283   delete genEventHeader;
 284   return;
 285
 286 }
 287 ////////////////////////////////////////////////////////////////////////
 288 void PrintVertex(TArrayF &primaryVertex)
 289 {
 290   cout <<"Vertex: "
 291        <<primaryVertex[0]<<" "
 292        <<primaryVertex[1]<<" "
 293        <<primaryVertex[2]<<" "<<endl;
 294   exit;
 295 }
 296 ////////////////////////////////////////////////////////////////////////