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

/******************************************************************************

  "AliITSOccupancy.C"
  
  this macro calculates the mean occupancy of each ITS layer, counting the
  "fired" digits of each module, and making two overall calculi:
     1) the calculus of the mean overall layer occupancy (as the ratio
	     between the total number of active digits and the total number of 
		  channels in the module;
	  2) a distribution of the occupancies as a funcion of z, to obtain some
		  kind of most significand data for this value, along the z axis
	
  The macro creates also the GIF of the TCanvas where the histograms are plotted
			  
  The arguments that must be inserted are the following:
     1) the answer to the following question (as a Bool_t)
		  "Can I get the total channels distribution from the file totals.root?"
		  if the answer is NO (false) these histograms are re-created and re-stored 
		  in that file.
			 
	  2) the name of the file to analyze; 
	     NOTE: in this argument you MUST omit the ".root" extension, because the 
		        macro uses this argument to recall the file, but also to name the
				  GIF & PS of the output canvas, in order to save the plot of the 
				  distributions to the disk. If you put, as argument, the string
				  "file.root", the macro will look for an unexisting file named
				  "file.root.root"...
                                                                         
  It is also possible (and recommended) to compile this macro, 
  to execute it faster. To do this, it is necessary to:
  1) make sure that de "#define __COMPILED__" instruction below is UNcommented
     and the "#define __NOCOMPILED__" instruction is commented (remember that
	  if you comment or uncomment both these instructions, the macro can't work)
  2) type, at the root prompt, the instruction 
     gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
	  to adjoin the ALICE include libraries to the main include directory
  3) load the function via the ".L ITSOccupancy.C++" statement (only the first
     time, because the first time the macro must be compiled). Frome the 
	  second time you use the macro, you must only load it via the 
	  ".L AliITSOccupancy.C+" instruction (note the number of '+' signs in each case
  4) call the function with only its name, es: AliITSOccupancy("HP", "galice.root")

 Author: Alberto Pulvirenti                                             
******************************************************************************/

#define __COMPILED__
//#define __NOCOMPILED__

#ifdef __COMPILED__
	#include <iostream.h>
	#include <TROOT.h>
	#include <TArrayI.h>
	#include <TCanvas.h>
	#include <TClassTable.h>
	#include <TClonesArray.h>
	#include <TFile.h>
	#include <TObject.h>
	#include <TObjArray.h>
	#include <TTree.h>
	#include <TMath.h>
	#include <TString.h>
	#include <TH1.h>
	#include <AliRun.h>
	#include <AliITS.h>
	#include <AliITSgeom.h>
	#include <AliITSDetType.h>
	#include <AliITSRecPoint.h>
	#include <AliITSdigit.h>
	#include <AliITShit.h>
	#include <AliITSmodule.h> 
	#include <AliITSsegmentation.h>
	#include <AliITSsegmentationSPD.h> 
	#include <AliITSsegmentationSDD.h>
	#include <AliITSsegmentationSSD.h>
#endif

Int_t AliITSOccupancy(char *name = "galice", Int_t evNum = 0) {       
		 
	// Evaluate the name of the file to open...
	Int_t filelen = strlen(name);
	Text_t *fileroot = new Text_t[filelen + 6];
	// ...and copy this name to the data ROOT file
	strcpy (fileroot, name);
	strcat (fileroot, ".root");
	
	// Open the Root input file containing the AliRun data and 
	// the Root output data file
	TFile *froot = (TFile*)gROOT->GetListOfFiles()->FindObject(fileroot);
	if (!froot) froot = new TFile(fileroot, "READ");
	// Load the gAlice shared libs if not already in memory and
	// clears an eventually existing AliRun object
#ifdef __NOCOMPILED__
	if (gClassTable->GetID("AliRun") < 0) {
    	gROOT->LoadMacro("loadlibs.C");
    	loadlibs();
  	}
#endif
   if (gAlice) {
      delete gAlice;
	   gAlice = 0;
   }
	// Get the AliRun object from file (if there's none, the macro is stopped)
	Int_t nparticles = 0;
	gAlice = (AliRun*)froot->Get("gAlice");
	if (gAlice) {
		cout << "Found an AliRun object in `" << fileroot << "'" << endl;
		// In this case, select the event number specified (default is 0)
		nparticles = gAlice->GetEvent(evNum);
		if (nparticles)
			cout << "\nNumber of particles   = " << nparticles << endl;
		else {
			cout << "NO PARTICLES!!!!" << endl;
			return 0;
		}
	}
	else {
		cout<<"Not found any AliRun object in `" << fileroot << "'" << endl;
		return 0;
	}
	// Initialize the ITS and its geometry
	AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
	AliITSgeom *gm = ITS->GetITSgeom();

	// Variables and objects definition and setting
	Float_t zmax[] = {20.0,20.0,25.0,35.0,49.5,54.0}; // z edges for TH1F (cm)
	Int_t nbins[]  = {40,40,50,70,99,108};             // bins number for TH1Fs
	
	// Histos for plotting occupancy distributions
	TH1F *above[6], *below[6];
	
	for (Int_t lay = 0; lay < 6; lay++) {
		Int_t nlads = gm->GetNladders(lay+1);
		Int_t ndets = gm->GetNdetectors(lay+1);
		Int_t dtype = lay / 2;
		Int_t minID = gm->GetModuleIndex(lay+1, 1, 1);
		Int_t maxID = gm->GetModuleIndex(lay+1, nlads, ndets);
		Text_t ffname[20];
		sprintf(ffname, "h_%d", lay+1);
		below[lay] = new TH1F(ffname, "Total z distribution of digits", nbins[lay], -zmax[lay], zmax[lay]);
		cout << "Evaluating total channels number of layer " << lay+1 << endl;
		for (Int_t I = minID; I <= maxID; I++) {		
			AliITSDetType *det = ITS->DetType(dtype);
			AliITSsegmentation *seg = det->GetSegmentationModel();
			Int_t NX = seg->Npx();
			if(lay == 2 || lay == 3) NX = 192;
			Int_t NZ = seg->Npz();
			//			cout << "ID: " << I << ", Layer: " << lay+1 << ", NX = " << NX << ", NZ = " << NZ << endl;
			for (Int_t ix = 0; ix <= NX; ix++) {
				for (Int_t iz = 0; iz <= NZ; iz++) {
					Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
					seg->DetToLocal(ix, iz, lx[0], lx[2]);
					gm->LtoG(I, lx, gx);
					below[lay]->Fill(gx[2]);
				}
			}
		}
		sprintf(ffname, "H_%d", lay+1);
		above[lay] = new TH1F(ffname, "histo", nbins[lay], -zmax[lay], zmax[lay]);
	}
		
	// Counting the hits, digits and recpoints contained in the ITS
	TTree *TD = gAlice->TreeD();
	ITS->ResetDigits();
	Float_t mean[6];
	Float_t average[6];

	for (Int_t L = 0; L < 6; L++) {
		
		cout << "Layer " << L + 1 << ": " << flush;				

		// To avoid two nested loops, are calculated 
		// the ID of the first and last module of the L
		// (remember the L goest from 0 to 5 (not from 1 to 6)
		Int_t first, last;
		first = gm->GetModuleIndex(L + 1, 1, 1);
		last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
				
		// Start loop on modules
		for (Int_t ID = first; ID <= last; ID++) {
			Int_t dtype = L / 2;
			AliITSDetType *det = ITS->DetType(dtype);
			AliITSsegmentation *seg = det->GetSegmentationModel();
			if (dtype == 2) seg->SetLayer(L+1);
			
			TD->GetEvent(ID);
			TClonesArray *digits_array = ITS->DigitsAddress(dtype);
			Int_t digits_num = digits_array->GetEntries();
			// Get the coordinates of the module
		  	for (Int_t j = 0; j < digits_num; j++) {
				Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
				AliITSdigit *digit = (AliITSdigit*)digits_array->UncheckedAt(j);
				Int_t iz=digit->fCoord1;  // cell number z
				Int_t ix=digit->fCoord2;  // cell number x
    			// Get local coordinates of the element (microns)
 		   	seg->DetToLocal(ix, iz, lx[0], lx[2]);
				gm->LtoG(ID, lx, gx);
				above[L]->Fill(gx[2]);
			}
		}
		
		Float_t occupied = above[L]->GetEntries();
		Float_t total = below[L]->GetEntries();
		cout << "Entries: " << occupied << ", " << total << endl;
		average[L] = 100.*occupied/total;
		above[L]->Divide(above[L], below[L], 100.0, 1.0);
		mean[L] = above[L]->GetSumOfWeights() / above[L]->GetNbinsX(); 
		cout.setf(ios::fixed);
		cout.precision(2);
		cout << " digits occupancy = " << mean[L] << "%" << endl;
		cout << " average digits occupancy = " << average[L] << "%" << endl;
	}
	
	TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
	view->Divide(2, 3);
	
	for (Int_t I = 1; I < 7; I++) {
		view->cd(I);
		Text_t title[50];
		sprintf(title, "Layer %d: %4.2f%c", I, mean[I-1], '%');
		title[6] = (Text_t)I + '0';
		above[I-1]->SetTitle(title);
		above[I-1]->SetStats(kFALSE);
		above[I-1]->SetXTitle("z (cm)");
		above[I-1]->SetYTitle("%");
		above[I-1]->Draw();
		view->Update();
	}
	Text_t *filegif  = new Text_t[filelen + 23];
	Text_t *fileps  = new Text_t[filelen + 23];
	strcpy (filegif, name);
	strcpy (fileps, name);
	strcat (filegif, "_digit_occupancy.gif");
	strcat (fileps, "_digit_occupancy.eps");
	view->SaveAs(filegif);
	view->SaveAs(fileps);
	TFile *fOc = new TFile("ITS_occupancy.root","NEW");
        fOc->cd();
	for (Int_t I = 0; I < 6; I++) above[I]->Write();
        fOc->Close();

	return 1;
}
Commit	Line	Data
af8e1c2d	1	/******************************************************************************
	2
	3	"AliITSOccupancy.C"
	4
	5	this macro calculates the mean occupancy of each ITS layer, counting the
	6	"fired" digits of each module, and making two overall calculi:
	7	1) the calculus of the mean overall layer occupancy (as the ratio
	8	between the total number of active digits and the total number of
	9	channels in the module;
	10	2) a distribution of the occupancies as a funcion of z, to obtain some
	11	kind of most significand data for this value, along the z axis
	12
	13	The macro creates also the GIF of the TCanvas where the histograms are plotted
	14
	15	The arguments that must be inserted are the following:
	16	1) the answer to the following question (as a Bool_t)
	17	"Can I get the total channels distribution from the file totals.root?"
	18	if the answer is NO (false) these histograms are re-created and re-stored
	19	in that file.
	20
	21	2) the name of the file to analyze;
	22	NOTE: in this argument you MUST omit the ".root" extension, because the
	23	macro uses this argument to recall the file, but also to name the
	24	GIF & PS of the output canvas, in order to save the plot of the
	25	distributions to the disk. If you put, as argument, the string
	26	"file.root", the macro will look for an unexisting file named
	27	"file.root.root"...
	28
	29	It is also possible (and recommended) to compile this macro,
	30	to execute it faster. To do this, it is necessary to:
	31	1) make sure that de "#define __COMPILED__" instruction below is UNcommented
	32	and the "#define __NOCOMPILED__" instruction is commented (remember that
	33	if you comment or uncomment both these instructions, the macro can't work)
	34	2) type, at the root prompt, the instruction
	35	gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
	36	to adjoin the ALICE include libraries to the main include directory
	37	3) load the function via the ".L ITSOccupancy.C++" statement (only the first
	38	time, because the first time the macro must be compiled). Frome the
	39	second time you use the macro, you must only load it via the
	40	".L AliITSOccupancy.C+" instruction (note the number of '+' signs in each case
	41	4) call the function with only its name, es: AliITSOccupancy("HP", "galice.root")
	42
	43	Author: Alberto Pulvirenti
	44	******************************************************************************/
	45
	46	#define __COMPILED__
	47	//#define __NOCOMPILED__
	48
	49	#ifdef __COMPILED__
	50	#include <iostream.h>
	51	#include <TROOT.h>
	52	#include <TArrayI.h>
	53	#include <TCanvas.h>
	54	#include <TClassTable.h>
	55	#include <TClonesArray.h>
	56	#include <TFile.h>
	57	#include <TObject.h>
	58	#include <TObjArray.h>
	59	#include <TTree.h>
	60	#include <TMath.h>
	61	#include <TString.h>
	62	#include <TH1.h>
	63	#include <AliRun.h>
	64	#include <AliITS.h>
65	#include <AliITSgeom.h>
66	#include <AliITSDetType.h>
67	#include <AliITSRecPoint.h>
68	#include <AliITSdigit.h>
69	#include <AliITShit.h>
70	#include <AliITSmodule.h>
71	#include <AliITSsegmentation.h>
72	#include <AliITSsegmentationSPD.h>
73	#include <AliITSsegmentationSDD.h>
74	#include <AliITSsegmentationSSD.h>
75	#endif
76
77	Int_t AliITSOccupancy(char *name = "galice", Int_t evNum = 0) {
78
79	// Evaluate the name of the file to open...
80	Int_t filelen = strlen(name);
81	Text_t *fileroot = new Text_t[filelen + 6];
82	// ...and copy this name to the data ROOT file
83	strcpy (fileroot, name);
84	strcat (fileroot, ".root");
85
86	// Open the Root input file containing the AliRun data and
87	// the Root output data file
88	TFile froot = (TFile)gROOT->GetListOfFiles()->FindObject(fileroot);
89	if (!froot) froot = new TFile(fileroot, "READ");
90	// Load the gAlice shared libs if not already in memory and
91	// clears an eventually existing AliRun object
92	#ifdef __NOCOMPILED__
93	if (gClassTable->GetID("AliRun") < 0) {
94	gROOT->LoadMacro("loadlibs.C");
95	loadlibs();
96	}
97	#endif
98	if (gAlice) {
99	delete gAlice;
100	gAlice = 0;
101	}
102	// Get the AliRun object from file (if there's none, the macro is stopped)
103	Int_t nparticles = 0;
104	gAlice = (AliRun*)froot->Get("gAlice");
105	if (gAlice) {
106	cout << "Found an AliRun object in `" << fileroot << "'" << endl;
107	// In this case, select the event number specified (default is 0)
108	nparticles = gAlice->GetEvent(evNum);
109	if (nparticles)
110	cout << "\nNumber of particles = " << nparticles << endl;
111	else {
112	cout << "NO PARTICLES!!!!" << endl;
113	return 0;
114	}
115	}
116	else {
117	cout<<"Not found any AliRun object in `" << fileroot << "'" << endl;
118	return 0;
119	}
120	// Initialize the ITS and its geometry
121	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
122	AliITSgeom *gm = ITS->GetITSgeom();
123
124	// Variables and objects definition and setting
125	Float_t zmax[] = {20.0,20.0,25.0,35.0,49.5,54.0}; // z edges for TH1F (cm)
126	Int_t nbins[] = {40,40,50,70,99,108}; // bins number for TH1Fs
127
128	// Histos for plotting occupancy distributions
129	TH1F above[6], below[6];
130
131	for (Int_t lay = 0; lay < 6; lay++) {
132	Int_t nlads = gm->GetNladders(lay+1);
133	Int_t ndets = gm->GetNdetectors(lay+1);
134	Int_t dtype = lay / 2;
135	Int_t minID = gm->GetModuleIndex(lay+1, 1, 1);
136	Int_t maxID = gm->GetModuleIndex(lay+1, nlads, ndets);
137	Text_t ffname[20];
138	sprintf(ffname, "h_%d", lay+1);
139	below[lay] = new TH1F(ffname, "Total z distribution of digits", nbins[lay], -zmax[lay], zmax[lay]);
140	cout << "Evaluating total channels number of layer " << lay+1 << endl;
141	for (Int_t I = minID; I <= maxID; I++) {
142	AliITSDetType *det = ITS->DetType(dtype);
143	AliITSsegmentation *seg = det->GetSegmentationModel();
144	Int_t NX = seg->Npx();
145	if(lay == 2 \|\| lay == 3) NX = 192;
146	Int_t NZ = seg->Npz();
147	// cout << "ID: " << I << ", Layer: " << lay+1 << ", NX = " << NX << ", NZ = " << NZ << endl;
148	for (Int_t ix = 0; ix <= NX; ix++) {
149	for (Int_t iz = 0; iz <= NZ; iz++) {
150	Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
151	seg->DetToLocal(ix, iz, lx[0], lx[2]);
152	gm->LtoG(I, lx, gx);
153	below[lay]->Fill(gx[2]);
154	}
155	}
156	}
157	sprintf(ffname, "H_%d", lay+1);
158	above[lay] = new TH1F(ffname, "histo", nbins[lay], -zmax[lay], zmax[lay]);
159	}
160
161	// Counting the hits, digits and recpoints contained in the ITS
162	TTree *TD = gAlice->TreeD();
163	ITS->ResetDigits();
164	Float_t mean[6];
165	Float_t average[6];
166
167	for (Int_t L = 0; L < 6; L++) {
168
169	cout << "Layer " << L + 1 << ": " << flush;
170
171	// To avoid two nested loops, are calculated
172	// the ID of the first and last module of the L
173	// (remember the L goest from 0 to 5 (not from 1 to 6)
174	Int_t first, last;
175	first = gm->GetModuleIndex(L + 1, 1, 1);
176	last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
177
178	// Start loop on modules
179	for (Int_t ID = first; ID <= last; ID++) {
180	Int_t dtype = L / 2;
181	AliITSDetType *det = ITS->DetType(dtype);
182	AliITSsegmentation *seg = det->GetSegmentationModel();
183	if (dtype == 2) seg->SetLayer(L+1);
184
185	TD->GetEvent(ID);
186	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
187	Int_t digits_num = digits_array->GetEntries();
188	// Get the coordinates of the module
189	for (Int_t j = 0; j < digits_num; j++) {
190	Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
191	AliITSdigit digit = (AliITSdigit)digits_array->UncheckedAt(j);
192	Int_t iz=digit->fCoord1; // cell number z
193	Int_t ix=digit->fCoord2; // cell number x
194	// Get local coordinates of the element (microns)
195	seg->DetToLocal(ix, iz, lx[0], lx[2]);
196	gm->LtoG(ID, lx, gx);
197	above[L]->Fill(gx[2]);
198	}
199	}
200
201	Float_t occupied = above[L]->GetEntries();
202	Float_t total = below[L]->GetEntries();
203	cout << "Entries: " << occupied << ", " << total << endl;
204	average[L] = 100.*occupied/total;
205	above[L]->Divide(above[L], below[L], 100.0, 1.0);
206	mean[L] = above[L]->GetSumOfWeights() / above[L]->GetNbinsX();
207	cout.setf(ios::fixed);
208	cout.precision(2);
209	cout << " digits occupancy = " << mean[L] << "%" << endl;
210	cout << " average digits occupancy = " << average[L] << "%" << endl;
211	}
212
213	TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
214	view->Divide(2, 3);
215
216	for (Int_t I = 1; I < 7; I++) {
217	view->cd(I);
218	Text_t title[50];
219	sprintf(title, "Layer %d: %4.2f%c", I, mean[I-1], '%');
220	title[6] = (Text_t)I + '0';
221	above[I-1]->SetTitle(title);
222	above[I-1]->SetStats(kFALSE);
223	above[I-1]->SetXTitle("z (cm)");
224	above[I-1]->SetYTitle("%");
225	above[I-1]->Draw();
226	view->Update();
227	}
228	Text_t *filegif = new Text_t[filelen + 23];
229	Text_t *fileps = new Text_t[filelen + 23];
230	strcpy (filegif, name);
231	strcpy (fileps, name);
232	strcat (filegif, "_digit_occupancy.gif");
233	strcat (fileps, "_digit_occupancy.eps");
234	view->SaveAs(filegif);
235	view->SaveAs(fileps);
236	TFile *fOc = new TFile("ITS_occupancy.root","NEW");
237	fOc->cd();
238	for (Int_t I = 0; I < 6; I++) above[I]->Write();
239	fOc->Close();
240
241	return 1;
242	}