[u/mrichter/AliRoot.git] / TPC / AliTPCTracking.C

////////////////////////////////////////////////////////////////////////
//
// 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;
} 
////////////////////////////////////////////////////////////////////////
Commit	Line	Data
171e422d	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	////////////////////////////////////////////////////////////////////////