[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:
//
//     20.11.2002 ... Changes due to a changed interface of AliTPCtracker. 
//                    Use Riostream.h instead of iostream.h
//
//     22.08.2002 ... first version
//
////////////////////////////////////////////////////////////////////////

#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"

#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);
    tracker->SetEventNumber(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	//
73f6cbc2	24	// 20.11.2002 ... Changes due to a changed interface of AliTPCtracker.
	25	// Use Riostream.h instead of iostream.h
	26	//
171e422d	27	// 22.08.2002 ... first version
	28	//
	29	////////////////////////////////////////////////////////////////////////
	30
	31	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
73f6cbc2	32	#include "Riostream.h"
171e422d	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;
73f6cbc2	176	AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
73f6cbc2	177	tracker->SetEventNumber(iEvent);
171e422d	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	////////////////////////////////////////////////////////////////////////