[u/mrichter/AliRoot.git] / ITS / AliITSFindTracksANN.C

#if !defined(__CINT__) || defined(__MAKECINT__)
  #include "Riostream.h"
  #include "TFile.h"
  #include "TTree.h"
  #include "TStopwatch.h"

  #include "AliRun.h"
  #include "AliMagF.h"
  #include "AliRunLoader.h"
  #include "AliTPCLoader.h"
  #include "AliITSLoader.h"
  #include "AliITS.h"
  #include "AliITSgeom.h"
  #include "AliITStrackerANN.h"
  #include "AliESD.h"
#endif

extern AliRun *gAlice;

Int_t AliITSFindTracksANN 
(Int_t nsectors = 10,      // number of azimutal sectors
 Int_t nev = 5)            // number of events to process
{

	// ==================================
	// ==== EVENT READING ===============
	// ==================================

	// Remove an already existing Run Loader
	if (gAlice) {
		delete gAlice->GetRunLoader();
		delete gAlice; 
		gAlice=0;
	}
	
	// Instance the new Run Loader
	AliRunLoader* rl = AliRunLoader::Open("galice.root");
	if (rl == 0x0) {
		cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
		return 3;
	}
	
	// Instance the ITS Loader
	AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
	if (itsl == 0x0) {
		cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
		return 4;
	}
	
	// Load the gAlice
	if (rl->LoadgAlice()) {
		cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
		delete rl;
		return 3;
	}
	
	// Set NECESSARY conversion constant for magnetic field
	AliKalmanTrack::SetConvConst
	(
		1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
	);
	
	// Get the ITS data & geometry & recpoints
	AliITS *dITS = (AliITS*)rl->GetAliRun()->GetDetector("ITS");
	if (!dITS) {
		cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
		return 6;
	}
	AliITSgeom *geom = dITS->GetITSgeom();
	itsl->LoadRecPoints("read");
	
	// ==================================
	// ==== CURVATURE CUT DEFINITION ====
	// ==================================

	// These values define the curvature cuts for all steps
	// within a sector.
	// For a greater clarity, the cuts are indicated in units
	// of transverse momentum (GeV/c) but these value have no
	// exact physical meaning, but are useful to understand
	// well what means a choice in the value of a certain
	// curvature constraint
	// NOTE: be careful to make sure that the 'ncuts' variable
	//       have the same value of the dimension of the allocated arrays

	Int_t ncuts;
	Double_t *p, *cut;

	ncuts = 6;
	p = new Double_t[6];
	cut = new Double_t[6];
	cut[ 0] = 0.0006;
	cut[ 1] = 0.0010;
	cut[ 2] = 0.0013;
	cut[ 3] = 0.0018;
	cut[ 4] = 0.0020;
	cut[ 5] = 0.0030;


	// ==========================
	// ==== OTHER PARAMETERS ====
	// ==========================

	
	Double_t helix[5]   = { 0.03, 0.02, 0.01, 0.03, 0.2 };
	Double_t theta2D[5] = { 1.5, 1.0, 1.0, 3.0, 10.0 };
	Double_t theta3D[5] = { 200., 200., 200., 200.0, 200.0 };
	
	
	/*
	Double_t helix[5]   = { 0.3, 0.3, 0.3, 0.5, 0.5 };
	Double_t theta2D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
	Double_t theta3D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
	*/
	
	Double_t temp   = 1.0;     // temperature parameter
	Double_t var    = 0.00001; // stabilization threshold

	Double_t exp    = 20.0;    // straight-line excitator
	Double_t gtoc   = 6.0;     // gain/cost contribution ratio

	Double_t min    = 0.4;     // minimum in random activation initialization
	Double_t max    = 0.6;     // maximum in random activation initialization
	Double_t actmin = 0.55;    // activation threshold for binary map conversion
	
	// Instance the Tracker
	AliITStrackerANN tracker(geom, 2);
	
	// Set cuts
	tracker.SetVertex(0.0, 0.0, 0.0);

	tracker.SetCuts(ncuts, cut, theta2D, theta3D, helix);
	tracker.SetTemperature(temp);
	tracker.SetVariationLimit(var);
	tracker.SetGainToCostRatio(gtoc);
	tracker.SetWeightExponent(exp);
	tracker.SetInitInterval(min, max);
	tracker.SetActThreshold(actmin);
	
	tracker.SetPolarInterval(45.0); 
		
	// Set number of events
	if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
	Int_t rc = 0;
	
	// Get ESD files
	TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
	if ((!itsf)||(!itsf->IsOpen())) {
		cerr << "Can't AliESDits.root !\n"; 
		return 1;
	}
	TFile *tpcf=TFile::Open("AliESDtpc.root");
	if ((!tpcf)||(!tpcf->IsOpen())) {
		cerr<<"Can't AliESDtpc.root !\n";
		return 1;
	}
	AliESD* event = new AliESD;
	TTree* esdTree = (TTree*) tpcf->Get("esdTree");
	if (!esdTree) {
	  cerr<<"no ESD tree found !\n";
	  return 1;
	}
	esdTree->SetBranchAddress("ESD", &event);
	
	// Loop on events
	TStopwatch timer; 
	for (Int_t i = 0; i < nev; i++) {
		cerr << "Processing event number: " << i << endl;
		esdTree->GetEvent(i);
		
		rl->GetEvent(i);
		
		// Get clusters from file
		TTree *cTree=itsl->TreeR();
		if (!cTree) {
			cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
			return 4;
		}
		
		// Load clusters in tracker
		tracker.LoadClusters(cTree);
		
		// Create array structure and arrange points in it
		tracker.CreateArrayStructure(nsectors);
		tracker.ArrangePoints("its.recpoints.txt");
		tracker.StoreOverallMatches();
		//tracker.PrintMatches(kFALSE);
		
		// ***************************
		//  NEURAL TRACKING OPERATION
		// ***************************
		
		Bool_t isStable = kFALSE;
		Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
		Int_t sectTracks = 0, totTracks = 0;

		// tracking through sectors
		cout << endl;
		Int_t sector, curv;
//		nUnits = tracker.CreateNetwork(0, ncuts - 1);
//		cout << endl << nUnits << " NEURONS CREATED" << endl;
//		return;
		
		for (sector = 0; sector < nsectors; sector++) {
			sectTracks = 0;
			for(curv = 0; curv < ncuts; curv++) {
				// units creation
				nUnits = tracker.CreateNetwork(sector, curv);
				if (!nUnits) {
					cout << "no neurons --> skipped" << endl;
					continue;
				}
				// units connection
				nLinks = tracker.LinkNeurons();
				if (!nLinks) {
					cout << "no connections --> skipped" << endl;
					continue;
				}
				// neural updating
				for (i = 0;; i++) {
					isStable = tracker.Update();
					if (isStable) break;
				}
				// tracks saving
				tracker.CleanNetwork();
				nTracks = tracker.StoreTracks();
				cout << nUnits << " units, ";
				cout << nLinks << " connections, ";
				cout << i << " cycles --> ";
				cout << nTracks << " tracks" << endl;
				sectTracks += nTracks;
			}
			cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
			totTracks += sectTracks;
		}
		tracker.ExportTracks("ITS.Neural.root");
		cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
		
		// ***************************
		
		// End of operations: unload clusters
		tracker.UnloadClusters();
		
	}
	timer.Stop(); 
	timer.Print();
	
	// Close files & delete objects
	delete event;
	tpcf->Close();
	itsf->Close();
	delete rl;
	return rc;
}
Commit	Line	Data
cec46807	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
cec46807	2	#include "Riostream.h"
8b462fd8	3	#include "TFile.h"
8b462fd8	4	#include "TTree.h"
cec46807	5	#include "TStopwatch.h"
	6
	7	#include "AliRun.h"
	8	#include "AliMagF.h"
	9	#include "AliRunLoader.h"
	10	#include "AliTPCLoader.h"
	11	#include "AliITSLoader.h"
	12	#include "AliITS.h"
	13	#include "AliITSgeom.h"
8b462fd8	14	#include "AliITStrackerANN.h"
cec46807	15	#include "AliESD.h"
	16	#endif
	17
	18	extern AliRun *gAlice;
	19
	20	Int_t AliITSFindTracksANN
	21	(Int_t nsectors = 10, // number of azimutal sectors
	22	Int_t nev = 5) // number of events to process
	23	{
	24
	25	// ==================================
	26	// ==== EVENT READING ===============
	27	// ==================================
	28
	29	// Remove an already existing Run Loader
	30	if (gAlice) {
	31	delete gAlice->GetRunLoader();
	32	delete gAlice;
	33	gAlice=0;
	34	}
	35
	36	// Instance the new Run Loader
	37	AliRunLoader* rl = AliRunLoader::Open("galice.root");
	38	if (rl == 0x0) {
	39	cerr<<"AliITSFindTracks.C : Can not open session RL=NULL"<< endl;
	40	return 3;
	41	}
	42
	43	// Instance the ITS Loader
	44	AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
	45	if (itsl == 0x0) {
	46	cerr<<"AliITSFindTracksV2.C : Can not get ITS loader"<<endl;
	47	return 4;
	48	}
	49
	50	// Load the gAlice
	51	if (rl->LoadgAlice()) {
	52	cerr<<"AliITSFindTracksV2.C : LoadgAlice returned error"<<endl;
	53	delete rl;
	54	return 3;
	55	}
	56
	57	// Set NECESSARY conversion constant for magnetic field
	58	AliKalmanTrack::SetConvConst
	59	(
	60	1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
	61	);
	62
	63	// Get the ITS data & geometry & recpoints
	64	AliITS dITS = (AliITS)rl->GetAliRun()->GetDetector("ITS");
	65	if (!dITS) {
	66	cerr<<"AliITSFindClusters.C : Can not find the ITS detector !"<<endl;
	67	return 6;
	68	}
	69	AliITSgeom *geom = dITS->GetITSgeom();
	70	itsl->LoadRecPoints("read");
	71
	72	// ==================================
	73	// ==== CURVATURE CUT DEFINITION ====
	74	// ==================================
	75
	76	// These values define the curvature cuts for all steps
	77	// within a sector.
	78	// For a greater clarity, the cuts are indicated in units
79	// of transverse momentum (GeV/c) but these value have no
80	// exact physical meaning, but are useful to understand
81	// well what means a choice in the value of a certain
82	// curvature constraint
83	// NOTE: be careful to make sure that the 'ncuts' variable
84	// have the same value of the dimension of the allocated arrays
85
86	Int_t ncuts;
87	Double_t p, cut;
88
89	ncuts = 6;
90	p = new Double_t[6];
91	cut = new Double_t[6];
92	cut[ 0] = 0.0006;
93	cut[ 1] = 0.0010;
94	cut[ 2] = 0.0013;
95	cut[ 3] = 0.0018;
96	cut[ 4] = 0.0020;
97	cut[ 5] = 0.0030;
98
99
100	// ==========================
101	// ==== OTHER PARAMETERS ====
102	// ==========================
103
104
105	Double_t helix[5] = { 0.03, 0.02, 0.01, 0.03, 0.2 };
106	Double_t theta2D[5] = { 1.5, 1.0, 1.0, 3.0, 10.0 };
107	Double_t theta3D[5] = { 200., 200., 200., 200.0, 200.0 };
108
109
110	/*
111	Double_t helix[5] = { 0.3, 0.3, 0.3, 0.5, 0.5 };
112	Double_t theta2D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
113	Double_t theta3D[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
114	*/
115
116	Double_t temp = 1.0; // temperature parameter
117	Double_t var = 0.00001; // stabilization threshold
118
119	Double_t exp = 20.0; // straight-line excitator
120	Double_t gtoc = 6.0; // gain/cost contribution ratio
121
122	Double_t min = 0.4; // minimum in random activation initialization
123	Double_t max = 0.6; // maximum in random activation initialization
124	Double_t actmin = 0.55; // activation threshold for binary map conversion
125
126	// Instance the Tracker
127	AliITStrackerANN tracker(geom, 2);
128
129	// Set cuts
130	tracker.SetVertex(0.0, 0.0, 0.0);
131
132	tracker.SetCuts(ncuts, cut, theta2D, theta3D, helix);
133	tracker.SetTemperature(temp);
134	tracker.SetVariationLimit(var);
135	tracker.SetGainToCostRatio(gtoc);
136	tracker.SetWeightExponent(exp);
137	tracker.SetInitInterval(min, max);
138	tracker.SetActThreshold(actmin);
139
140	tracker.SetPolarInterval(45.0);
141
142	// Set number of events
143	if (nev > rl->GetNumberOfEvents()) nev = rl->GetNumberOfEvents();
144	Int_t rc = 0;
145
146	// Get ESD files
147	TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
148	if ((!itsf)\|\|(!itsf->IsOpen())) {
149	cerr << "Can't AliESDits.root !\n";
150	return 1;
151	}
152	TFile *tpcf=TFile::Open("AliESDtpc.root");
153	if ((!tpcf)\|\|(!tpcf->IsOpen())) {
154	cerr<<"Can't AliESDtpc.root !\n";
155	return 1;
156	}
8b462fd8	157	AliESD* event = new AliESD;
	158	TTree* esdTree = (TTree*) tpcf->Get("esdTree");
	159	if (!esdTree) {
	160	cerr<<"no ESD tree found !\n";
	161	return 1;
	162	}
	163	esdTree->SetBranchAddress("ESD", &event);
cec46807	164
cec46807	165	// Loop on events
cec46807	166	TStopwatch timer;
cec46807	167	for (Int_t i = 0; i < nev; i++) {
cec46807	168	cerr << "Processing event number: " << i << endl;
8b462fd8	169	esdTree->GetEvent(i);
cec46807	170
	171	rl->GetEvent(i);
	172
	173	// Get clusters from file
	174	TTree *cTree=itsl->TreeR();
	175	if (!cTree) {
	176	cerr<<"AliITSFindTracksANN.C : NO clusters tree!" << endl;
	177	return 4;
	178	}
	179
	180	// Load clusters in tracker
	181	tracker.LoadClusters(cTree);
	182
	183	// Create array structure and arrange points in it
	184	tracker.CreateArrayStructure(nsectors);
	185	tracker.ArrangePoints("its.recpoints.txt");
	186	tracker.StoreOverallMatches();
	187	//tracker.PrintMatches(kFALSE);
	188
	189	// ***************************
	190	// NEURAL TRACKING OPERATION
	191	// ***************************
	192
	193	Bool_t isStable = kFALSE;
	194	Int_t i, nUnits = 0, nLinks = 0, nTracks = 0;
	195	Int_t sectTracks = 0, totTracks = 0;
	196
	197	// tracking through sectors
	198	cout << endl;
	199	Int_t sector, curv;
	200	// nUnits = tracker.CreateNetwork(0, ncuts - 1);
	201	// cout << endl << nUnits << " NEURONS CREATED" << endl;
	202	// return;
	203
	204	for (sector = 0; sector < nsectors; sector++) {
	205	sectTracks = 0;
	206	for(curv = 0; curv < ncuts; curv++) {
	207	// units creation
	208	nUnits = tracker.CreateNetwork(sector, curv);
	209	if (!nUnits) {
	210	cout << "no neurons --> skipped" << endl;
	211	continue;
	212	}
	213	// units connection
	214	nLinks = tracker.LinkNeurons();
	215	if (!nLinks) {
	216	cout << "no connections --> skipped" << endl;
	217	continue;
	218	}
	219	// neural updating
	220	for (i = 0;; i++) {
	221	isStable = tracker.Update();
	222	if (isStable) break;
	223	}
	224	// tracks saving
	225	tracker.CleanNetwork();
	226	nTracks = tracker.StoreTracks();
	227	cout << nUnits << " units, ";
	228	cout << nLinks << " connections, ";
	229	cout << i << " cycles --> ";
	230	cout << nTracks << " tracks" << endl;
	231	sectTracks += nTracks;
	232	}
	233	cout << "\n >>> Total tracks found in sector: " << sectTracks << endl;
234	totTracks += sectTracks;
235	}
236	tracker.ExportTracks("ITS.Neural.root");
237	cout << "\n\n--- Totally found " << totTracks << " tracks ---\n\n" << endl;
238
239	// ***************************
240
241	// End of operations: unload clusters
242	tracker.UnloadClusters();
243
cec46807	244	}
	245	timer.Stop();
	246	timer.Print();
	247
	248	// Close files & delete objects
8b462fd8	249	delete event;
cec46807	250	tpcf->Close();
	251	itsf->Close();
	252	delete rl;
	253	return rc;
	254	}
	255
	256