1 ////////////////////////////////////////////////////////////////////////
6 // author: Jiri Chudoba based on code of Jourij Belikov
9 // reconstructs of tracks in TPC inthe following steps:
10 // TPC cluster finding
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)
20 // default file names correspond to pp production (2002-04)
24 // 20.11.2002 ... Changes due to a changed interface of AliTPCtracker.
25 // Use Riostream.h instead of iostream.h
27 // 22.08.2002 ... first version
29 ////////////////////////////////////////////////////////////////////////
31 #if !defined(__CINT__) || defined(__MAKECINT__)
32 #include "Riostream.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"
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);
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();
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");
74 ////////////////////////////////////////////////////////////////////////
75 Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
77 TString fileNameDigits,
78 TString fileNameClusters,
79 TString fileNameTracks) {
81 SetFieldFactor(fileNameHits,kFALSE);
83 // ********** Find TPC clusters *********** //
84 if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
85 cerr<<"Failed to get TPC clusters: !\n";
89 // ********** Find TPC tracks *********** //
90 if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
91 cerr<<"Failed to get TPC tracks !\n";
98 ////////////////////////////////////////////////////////////////////////
99 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
100 TString fileNameDigits, TString fileNameClusters) {
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";
112 if (!fileClusters->IsOpen()) {
113 cerr<<"Cannnot open "<<fileNameClusters<<" !\n";
117 rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
120 fileClusters->Close();
123 if (gDEBUG>1) gBenchmark->Stop(name);
124 if (gDEBUG>1) gBenchmark->Show(name);
128 ////////////////////////////////////////////////////////////////////////
129 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
130 TFile* fileDigits, TFile* fileClusters,
131 AliTPCParam* paramTPC) {
134 if (!paramTPC) paramTPC = LoadTPCParam(fileDigits);
135 if (!paramTPC) return 1;
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);
143 ////////////////////////////////////////////////////////////////////////
144 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
145 TString fileNameClusters, TString fileNameTracks) {
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);
154 rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
156 fileClusters->Close();
160 if (gDEBUG>1) gBenchmark->Stop(name);
161 if (gDEBUG>1) gBenchmark->Show(name);
165 ////////////////////////////////////////////////////////////////////////
166 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
167 TFile *fileClusters, TFile * fileTracks,
168 AliTPCParam* paramTPC) {
171 if (!paramTPC) paramTPC = LoadTPCParam(fileClusters);
172 if (!paramTPC) return 1;
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);
179 FindVertex(iEvent,vertex);
180 tracker->SetVertex(vertex);
182 rc = tracker->Clusters2Tracks(0,fileTracks);
188 ////////////////////////////////////////////////////////////////////////
189 Int_t SetFieldFactor() {
191 AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
192 if (gDEBUG > 2) cout<<"Magnetic field in kGauss: "<<gAlice->Field()->SolenoidField()<<endl;
195 ////////////////////////////////////////////////////////////////////////
196 Int_t SetFieldFactor(TFile *file, Bool_t deletegAlice) {
198 if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
199 cerr<<"gAlice has not been found in file "<<file->GetName();
202 Int_t rc = SetFieldFactor();
209 ////////////////////////////////////////////////////////////////////////
210 Int_t SetFieldFactor(TString fileName, Bool_t closeFile) {
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();
218 ////////////////////////////////////////////////////////////////////////
219 AliTPCParam* LoadTPCParam(TFile *file) {
222 sprintf(paramName,"75x40_100x60_150x60");
223 AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
225 if (gDEBUG > 1) cout<<"TPC parameters "<<paramName<<" found."<<endl;
227 cerr<<"TPC parameters not found. Create new (they may be incorrect)."
229 paramTPC = new AliTPCParamSR;
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.
239 // char paramName[50];
240 // sprintf(paramName,"75x40_100x60");
241 // AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
243 // cout<<"TPC parameters "<<paramName<<" found."<<endl;
245 // sprintf(paramName,"75x40_100x60_150x60");
246 // paramTPC=(AliTPCParam*)in->Get(paramName);
248 // cout<<"TPC parameters "<<paramName<<" found."<<endl;
250 // cerr<<"TPC parameters not found. Create new (they may be incorrect)."
252 // paramTPC = new AliTPCParamSR;
259 ////////////////////////////////////////////////////////////////////////
260 void FindVertex(Int_t eventNr, Double_t *vertex) {
262 vertex[0] = vertex[1] = vertex[2] = 0.;
264 cerr<<"gAlice was not found! Using vertex position (0,0,0).\n";
268 gAlice->GetEvent(eventNr);
269 AliHeader *header = gAlice->GetHeader();
271 cerr<<"header was not found!\n";
274 AliGenEventHeader* genEventHeader = header->GenEventHeader();
275 if (!genEventHeader) {
276 cerr<<"AliGenEventHeader was not found!\n";
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]);
287 delete genEventHeader;
291 ////////////////////////////////////////////////////////////////////////
292 void PrintVertex(TArrayF &primaryVertex)
295 <<primaryVertex[0]<<" "
296 <<primaryVertex[1]<<" "
297 <<primaryVertex[2]<<" "<<endl;
300 ////////////////////////////////////////////////////////////////////////