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