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

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

  "ITSOccupancy.C"
  
  this macro calculates the mean occupancy of each ITS layer,       
  making also a distribution in function of the z-value of the           
  "fired" digits for each layer                                              
                                                                         
  Because of evaluation problems, it is necessary to permit the          
  analysis for both hits, digits and recpoints, so it's provided an     
  option data which must contain                                         
  - "H" for hits                                                         
  - "D" for digits                                                       
  - "R" for recpoints                                                    
  but it's possible to use several option in the same time               
  (like "HD", "HRD", and so on...)                                       
                                                                         
  to avoid the problem of drawing unrequested windows, the graph are     
  plotted only if it is explictly requested with the option "P", that    
  can be inserted together to the ones described above.                  
                                                                         
  It is also possible 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 ITSOccupancy.C+" instruction (note the number of '+' signs in each case
  4) call the function with only its name, es: ITSOccupancy("HP", "galice.root")

 NOTE 1: option interpretation is case-insensitive ("h" eqv "H", etc.)    
 NOTE 2: if you insert a filename as argument, it must be inserted only the name
         (the macro attaches the extension ".root" to it, every time)

 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 ITSOccupancy(char* opt, char *name = "galice", Int_t evNum = 0) {       

	extern void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter);
	extern void GetModuleHits (TObject* its, Int_t ID, TH1D *counter);
	extern void GetModuleRecPoints (TObject *its, Int_t ID, TH1D *counter);
	extern Int_t DType(Int_t layer);
		 
	// Evaluating options
	TString option(opt);
	Bool_t HITS = (option.Contains("h") || option.Contains("H"));
	Bool_t DIGS = (option.Contains("d") || option.Contains("D"));
	Bool_t RECS = (option.Contains("r") || option.Contains("R"));
	Bool_t PLOT = (option.Contains("p") || option.Contains("P"));
	
	Text_t filename[100];
	char  gifH[100], gifD[100], gifR[100];
	strcpy (filename, name);
	strcpy (gifH, name);
	strcpy (gifD, name);
	strcpy (gifR, name);
	strcat (filename, ".root");
	strcat (gifH, "H.gif");
	strcat (gifD, "D.gif");
	strcat (gifR, "R.gif");
	
	if (!HITS && !DIGS && !RECS) {
		cout << "\n\n--- NO DATA-TYPE SPECIFIED!!! ---\n\n";
		return 0;
	}
	
	// --------------------------------
	//  1. ALICE objects setting phase
	// --------------------------------
	
	// Load the gAlice shared libs if not already in memory
#ifdef __NOCOMPILED__
	if (gClassTable->GetID("AliRun") < 0) {
    	gROOT->LoadMacro("loadlibs.C");
    	loadlibs();
  	}
#endif
   if(gAlice){
      delete gAlice;
	   gAlice=0;
   }
	// Open the Root input file containing the AliRun data 
	// (default name is "galice.root")
	TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
	if (!file) file = new TFile(filename);
	file->ls();
	// Get the AliRun object from file (if there's none, the macro is stopped)
	cout << "\nSearching in '" << filename << "'" << endl;
	gAlice = (AliRun*)file->Get("gAlice");
	if (gAlice)
		cout << "Ok, I found an AliRun object..." << endl;
	else {
		cout<<"Sorry, there isn't any AliRun object here..." << endl;
		return 0;
	}
	// Select the event number specified. Default is 0.
	Int_t nparticles = gAlice->GetEvent(evNum);
	cout << "\nNumber of particles   = " << nparticles << endl;
	if (!nparticles) {
		cout << "With no particles I can't do much..." << endl;
		return 0;
	}
	// Initialize the ITS and its geometry
	AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
	AliITSgeom *gm = ITS->GetITSgeom();
	// Fill the AliITSmodule objects (only for hits)
	Int_t nmodules;
	if (HITS) {
	  	ITS->InitModules(-1, nmodules);
  		cout << "Number of ITS modules = " << nmodules << endl;
		cout << "\nFilling modules (it takes a while, now)..." << flush;
		ITS->FillModules(0, 0, nmodules, " ", " ");
	  	cout << "DONE!" << endl;	
	}

	// -------------------------------------------------
	//  B. Variables and objects definition and setting
	// -------------------------------------------------
	
	Double_t zmax[6] = {16.5, 16.5, 22.2, 29.7, 45.1, 50.8};  // max z measured (cm)
	Double_t binw[] = {2., 2., 2., 2., 5., 5.}; // bin width for TH1Ds (cm)
	
	Double_t tot[6];  // total number of channels 
	
	// Histos for plotting hits [1], digits [2] and points [3], and divisor [0]
	TH1D *hist[4][6];
	
	// Histograms initialization:
	// Using the AliITSsegmentation object, here is deduced
	// the maximum z that can be measurable for each layer, simply
	// obtaining the globar coordinates of the first digit of the first module 
	// of each layer (that is placed at the maximum z).
	for (Int_t i = 0; i < 6; i++) {
		AliITSDetType *det = ITS->DetType(DType(i));	
		AliITSsegmentation *seg = det->GetSegmentationModel();
		tot[i] = gm->GetNladders(i+1) * gm->GetNdetectors(i+1) * seg->GetNPads();
		cout.setf(ios::fixed);
		cout << tot[i] << endl;
		Text_t name[7];
		for (Int_t j = 0; j < 4; j++) {
			name[0] = 'h';
			name[1] = 'i';
			name[2] = 's';
			name[3] = 't';
			name[4] = '0' + j;
			name[5] = '0' + i;
			name[6] = '\0';
			Int_t nbins = (Int_t)(2. * zmax[i] / binw[i]);
			hist[j][i] = new TH1D(name, "histo", nbins, -zmax[i], zmax[i]);
		}
	}
	// Here is put the value of the density for each bin of each layer's TH1D
	for (Int_t i = 0; i < 6; i++) {
		for (Int_t j = 0; j < hist[0][i]->GetNbinsX(); j++) {
			Double_t z = hist[0][i]->GetBinLowEdge(j);
			Double_t bincontent = tot[i] / (zmax[i] * 2.0) * hist[0][i]->GetBinWidth(j);
			hist[0][i]->Fill(z + 0.1, bincontent / 100.0);
			z += hist[0][i]->GetBinWidth(j);
		}
	}
	
	// --------------------------------------------------------------
	//  Counting the hits, digits and recpoints contained in the ITS
	// --------------------------------------------------------------

	for (Int_t L = 0; L < 6; L++) {
		
		// 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
		cout << "Examinating layer " << L + 1 << " ... " << flush;
		for (Int_t ID = first; ID <= last; ID++) {
			
			// Hits analysis (if requested)
			if (HITS)
				GetModuleHits(ITS, ID, hist[2][L]);
			
			// Digits analysis (if requested)
			if (DIGS)
				GetModuleDigits(ITS, ID, DType(L), hist[1][L]);
			
			// Recpoints analysis (if requested)
			if (RECS) 
				GetModuleRecPoints(ITS, ID, hist[3][L]);
		}
		
		hist[1][L]->Divide(hist[0][L]);
		hist[2][L]->Divide(hist[0][L]);
		hist[3][L]->Divide(hist[0][L]);
		cout << "DONE" << endl;
	}

	// Write on the screen the total means
	cout << endl;
	cout << "********* MEAN OCCUPANCIES *********" << endl;
	for (Int_t L = 0; L < 6; L++) {
		Double_t mean[3];
		cout << " LAYER " << L << ": " << endl;
		for (Int_t i = 0; i < 3; i++) {
			mean[i] = hist[i+1][L]->GetSumOfWeights() / hist[i+1][L]->GetNbinsX(); 
			Text_t title[50];
			sprintf(title, " - Layer %d, mean %4.2f - ", L + 1, mean[i]);
			hist[i+1][L]->SetTitle(title);
			cout.setf(ios::fixed);
			cout.precision(3);
			if (HITS && i == 1) {
				cout << "     hits occupancy = ";
				if (mean[i] < 10.0) cout << ' ';
				cout << mean[i] << "%" << endl;
			}
			else if (DIGS && i == 0) {
				cout << "   digits occupancy = ";
				if (mean[i] < 10.0) cout << ' ';
				cout << mean[i] << "%" << endl;
			}
			if (RECS && i == 2) {
				cout << "   recpts occupancy = ";
				if (mean[i] < 10.0) cout << ' ';
				cout << mean[i] << "%" << endl;
			}
		}
		cout << "------------------------------------" << endl;
	}
	
	// ----------------------
	//  Plots (if requested)
	// ----------------------
	if (!PLOT) return 1;
	
	TCanvas *view[4];
	if (HITS) {
		view[2] = new TCanvas("viewH", "Occupancy view (HITS)", 0, 0, 1050, 700);
		view[2]->Divide(3,2, 0.001, 0.001);
	}	
	if (DIGS) {
		view[1] = new TCanvas("viewD", "Occupancy view (DIGITS)", 20, 40, 1050, 700);
		view[1]->Divide(3,2, 0.001, 0.001);
	}
	if (RECS) {
		view[3] = new TCanvas("viewR", "Occupancy view (RECPOINTS)", 40, 60, 1050, 700);
		view[3]->Divide(3,2, 0.001, 0.001);
	}
	
	for (Int_t L = 0; L < 6; L++) {
		for (Int_t i = 1; i <= 3; i++) {
			if (DIGS && i == 1) 
				hist[i][L]->SetFillColor(kBlue);
			else if (HITS && i == 2) 
				hist[i][L]->SetFillColor(kRed);
			else if (RECS && i == 3) 
				hist[i][L]->SetFillColor(kGreen);
			
			if ((HITS && i == 2) || (DIGS && i == 1) || (RECS && i == 3)) {			
				view[i]->cd(L+1);			
				hist[i][L]->SetStats(kFALSE);
				hist[i][L]->Draw();
				hist[i][L]->GetXaxis()->SetTitle("zeta");
				hist[i][L]->GetYaxis()->SetTitle("%");
				view[i]->Update();
			}
		}
	}
	
	if (HITS) view[2]->SaveAs(gifH);
	if (DIGS) view[1]->SaveAs(gifD);
	if (RECS) view[3]->SaveAs(gifR);
	
	return 1;
}

void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter) {
	// First of all, the macro selects the specified module,
	// then gets the array of recpoints in it and their number.
	AliITS *ITS = (AliITS*)its;
	TTree *TD = gAlice->TreeD();
	ITS->ResetDigits();
	TD->GetEvent(ID);
	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
	AliITSgeom *gm = ITS->GetITSgeom();
	AliITSDetType *det = ITS->DetType(dtype);	
	AliITSsegmentation *seg = det->GetSegmentationModel();	
	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->GetPadCxz(ix, iz, lx[0], lx[2]);
		if (dtype != 1) {
			// !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
			if (!dtype) {
				lx[0] -= seg->Dx() / 2.0;
				lx[2] -= seg->Dz() / 2.0;
			}
			lx[0] /= 10000.0; // convert from microns to cm (if not SDD)
			lx[2] /= 10000.0; // convert from microns to cm (if not SDD)
		}
		gm->LtoG(ID, lx, gx);
		counter->Fill(gx[2]);
	}
}

/*void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter) {
	// First of all, the macro selects the specified module,
	// then gets the array of recpoints in it and their number.
	AliITS *ITS = (AliITS*)its;
	TTree *TD = gAlice->TreeD();
	ITS->ResetDigits();
	TD->GetEvent(ID);
	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
	AliITSgeom *gm = ITS->GetITSgeom();
	AliITSDetType *det = ITS->DetType(dtype);	
	AliITSsegmentation *seg = det->GetSegmentationModel();	
	TArrayI ssdone(5000);  // used to match p and n side digits of strips
	TArrayI pair(5000);    // as above 	
	Int_t digits_num = digits_array->GetEntries();
	
	// Get the coordinates of the module
	if (dtype == 2) {
		for (Int_t j = 0; j < digits_num; j++) {
			ssdone[j] = 0;			
			pair[j] = 0;
		}
	}
  	for (Int_t j = 0; j < digits_num; j++) {
		Double_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)
		if(dtype < 2) {
    		seg->GetPadCxz(ix, iz, lx[0], lx[2]);
			if (dtype == 0) {
				// !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
				lx[0] -= seg->Dx() / 2.0;
				lx[2] -= seg->Dz() / 2.0;
				lx[0] /= 10000.0; // convert from microns to cm (SPD)
				lx[2] /= 10000.0; // convert from microns to cm (SPD)
			}
			gm->LtoG(ID, lx, gx);
			counter->Fill(gx[2]);
		}
	   else {
			// SSD: if iz==0 ---> N side; if iz==1 P side
      	if (ssdone[j] == 0) {
				ssdone[j]=1;
				pair[j]=-1;
				Bool_t pside = (iz == 1);
				Bool_t impaired = kTRUE;
				Int_t pstrip = 0;
				Int_t nstrip = 0;
				if (pside) pstrip = ix; else nstrip = ix;
				for (Int_t k = 0; k < digits_num; k++) {
					if (ssdone[k] == 0 && impaired) {
						AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
						if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
							ssdone[k]=2;
							pair[j]=k;
							if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
							impaired=kFALSE;
						}
					}
				}
        		if (!impaired) {
					seg->GetPadCxz(pstrip, nstrip, lx[0], lx[2]); 
					lx[0] /= 10000.0;
					lx[2] /= 10000.0;
					gm->LtoG(ID, lx, gx);
					counter->Fill(gx[2]);
				}
			}
		}
	}
}*/

void GetModuleHits (TObject* its, Int_t ID, TH1D *counter) {	
	// First of all, the macro selects the specified module,
	// then gets the array of hits in it and their number.
	AliITS *ITS = (AliITS*) its;
	AliITSmodule *module = ITS->GetModule(ID);
	TObjArray *hits_array = module->GetHits();
	Int_t hits_num = hits_array->GetEntriesFast();	
	// next, fills the counters with the hits coordinates' calcula
	for (Int_t j = 0; j < hits_num; j++) {
		AliITShit *hit = (AliITShit*) hits_array->At(j);
		Double_t Z = hit->GetZG();
		if (!hit->StatusEntering()) counter->Fill(Z);
	}
}

void GetModuleRecPoints (TObject *its, Int_t ID, TH1D *counter) {
	// First of all, the macro selects the specified module,
	// then gets the array of recpoints in it and their number.
	AliITS *ITS = (AliITS*) its;
	AliITSgeom *gm = ITS->GetITSgeom();	
	TTree *TR = gAlice->TreeR();
	if (!TR || TR->GetEntries() == 0) {
		cout << "The ITS wasn't clustered..." << endl;
		return;
	}
	ITS->ResetRecPoints();
	TR->GetEvent(ID);
	TClonesArray *recs_array = ITS->RecPoints();
	Int_t recs_num = recs_array->GetEntries();

	for(Int_t j = 0; j < recs_num; j++) {
		Double_t lx[3] = {0.,0.,0.}, gx[3] = {0.,0.,0.};
		AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
	   lx[0] = recp->GetX();
		lx[1] = 0.0;
		lx[2] = recp->GetZ();
		gm->LtoG(ID, lx, gx);
		counter->Fill(gx[2]);
	}
}

Int_t DType(Int_t layer) {
	if (layer == 0 || layer == 1)
		return 0;
	else if (layer == 2 || layer == 3)
		return 1;
	else 
		return 2;
}
Commit	Line	Data
2bba24f4	1	/******************************************************************************
	2
	3	"ITSOccupancy.C"
	4
	5	this macro calculates the mean occupancy of each ITS layer,
	6	making also a distribution in function of the z-value of the
	7	"fired" digits for each layer
	8
	9	Because of evaluation problems, it is necessary to permit the
	10	analysis for both hits, digits and recpoints, so it's provided an
	11	option data which must contain
	12	- "H" for hits
	13	- "D" for digits
	14	- "R" for recpoints
	15	but it's possible to use several option in the same time
	16	(like "HD", "HRD", and so on...)
	17
	18	to avoid the problem of drawing unrequested windows, the graph are
	19	plotted only if it is explictly requested with the option "P", that
	20	can be inserted together to the ones described above.
	21
	22	It is also possible to compile this macro, to execute it faster.
	23	To do this, it is necessary to:
	24	1) make sure that de "#define __COMPILED__" instruction below is UNcommented
	25	and the "#define __NOCOMPILED__" instruction is commented (remember that
	26	if you comment or uncomment both these instructions, the macro can't work)
	27	2) type, at the root prompt, the instruction
	28	gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
	29	to adjoin the ALICE include libraries to the main include directory
	30	3) load the function via the ".L ITSOccupancy.C++" statement (only the first
	31	time, because the first time the macro must be compiled). Frome the
	32	second time you use the macro, you must only load it via the
	33	".L ITSOccupancy.C+" instruction (note the number of '+' signs in each case
	34	4) call the function with only its name, es: ITSOccupancy("HP", "galice.root")
	35
	36	NOTE 1: option interpretation is case-insensitive ("h" eqv "H", etc.)
	37	NOTE 2: if you insert a filename as argument, it must be inserted only the name
	38	(the macro attaches the extension ".root" to it, every time)
	39
	40	Author: Alberto Pulvirenti
	41	******************************************************************************/
	42
	43	#define __COMPILED__
	44	//#define __NOCOMPILED__
	45
	46	#ifdef __COMPILED__
	47	#include <iostream.h>
	48	#include <TROOT.h>
	49	#include <TArrayI.h>
	50	#include <TCanvas.h>
	51	#include <TClassTable.h>
	52	#include <TClonesArray.h>
	53	#include <TFile.h>
	54	#include <TObject.h>
	55	#include <TObjArray.h>
	56	#include <TTree.h>
	57	#include <TMath.h>
	58	#include <TString.h>
	59	#include <TH1.h>
	60	#include <AliRun.h>
	61	#include <AliITS.h>
	62	#include <AliITSgeom.h>
	63	#include <AliITSDetType.h>
	64	#include <AliITSRecPoint.h>
65	#include <AliITSdigit.h>
66	#include <AliITShit.h>
67	#include <AliITSmodule.h>
68	#include <AliITSsegmentation.h>
69	#include <AliITSsegmentationSPD.h>
70	#include <AliITSsegmentationSDD.h>
71	#include <AliITSsegmentationSSD.h>
72	#endif
73
74	Int_t ITSOccupancy(char* opt, char *name = "galice", Int_t evNum = 0) {
75
76	extern void GetModuleDigits(TObject its, Int_t ID, Int_t dtype, TH1D counter);
77	extern void GetModuleHits (TObject* its, Int_t ID, TH1D *counter);
78	extern void GetModuleRecPoints (TObject its, Int_t ID, TH1D counter);
79	extern Int_t DType(Int_t layer);
80
81	// Evaluating options
82	TString option(opt);
83	Bool_t HITS = (option.Contains("h") \|\| option.Contains("H"));
84	Bool_t DIGS = (option.Contains("d") \|\| option.Contains("D"));
85	Bool_t RECS = (option.Contains("r") \|\| option.Contains("R"));
86	Bool_t PLOT = (option.Contains("p") \|\| option.Contains("P"));
87
88	Text_t filename[100];
89	char gifH[100], gifD[100], gifR[100];
90	strcpy (filename, name);
91	strcpy (gifH, name);
92	strcpy (gifD, name);
93	strcpy (gifR, name);
94	strcat (filename, ".root");
95	strcat (gifH, "H.gif");
96	strcat (gifD, "D.gif");
97	strcat (gifR, "R.gif");
98
99	if (!HITS && !DIGS && !RECS) {
100	cout << "\n\n--- NO DATA-TYPE SPECIFIED!!! ---\n\n";
101	return 0;
102	}
103
104	// --------------------------------
105	// 1. ALICE objects setting phase
106	// --------------------------------
107
108	// Load the gAlice shared libs if not already in memory
109	#ifdef __NOCOMPILED__
110	if (gClassTable->GetID("AliRun") < 0) {
111	gROOT->LoadMacro("loadlibs.C");
112	loadlibs();
113	}
114	#endif
115	if(gAlice){
116	delete gAlice;
117	gAlice=0;
118	}
119	// Open the Root input file containing the AliRun data
120	// (default name is "galice.root")
121	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(filename);
122	if (!file) file = new TFile(filename);
123	file->ls();
124	// Get the AliRun object from file (if there's none, the macro is stopped)
125	cout << "\nSearching in '" << filename << "'" << endl;
126	gAlice = (AliRun*)file->Get("gAlice");
127	if (gAlice)
128	cout << "Ok, I found an AliRun object..." << endl;
129	else {
130	cout<<"Sorry, there isn't any AliRun object here..." << endl;
131	return 0;
132	}
133	// Select the event number specified. Default is 0.
134	Int_t nparticles = gAlice->GetEvent(evNum);
135	cout << "\nNumber of particles = " << nparticles << endl;
136	if (!nparticles) {
137	cout << "With no particles I can't do much..." << endl;
138	return 0;
139	}
140	// Initialize the ITS and its geometry
141	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
142	AliITSgeom *gm = ITS->GetITSgeom();
143	// Fill the AliITSmodule objects (only for hits)
144	Int_t nmodules;
145	if (HITS) {
146	ITS->InitModules(-1, nmodules);
147	cout << "Number of ITS modules = " << nmodules << endl;
148	cout << "\nFilling modules (it takes a while, now)..." << flush;
149	ITS->FillModules(0, 0, nmodules, " ", " ");
150	cout << "DONE!" << endl;
151	}
152
153	// -------------------------------------------------
154	// B. Variables and objects definition and setting
155	// -------------------------------------------------
156
157	Double_t zmax[6] = {16.5, 16.5, 22.2, 29.7, 45.1, 50.8}; // max z measured (cm)
158	Double_t binw[] = {2., 2., 2., 2., 5., 5.}; // bin width for TH1Ds (cm)
159
160	Double_t tot[6]; // total number of channels
161
162	// Histos for plotting hits [1], digits [2] and points [3], and divisor [0]
163	TH1D *hist[4][6];
164
165	// Histograms initialization:
166	// Using the AliITSsegmentation object, here is deduced
167	// the maximum z that can be measurable for each layer, simply
168	// obtaining the globar coordinates of the first digit of the first module
169	// of each layer (that is placed at the maximum z).
170	for (Int_t i = 0; i < 6; i++) {
171	AliITSDetType *det = ITS->DetType(DType(i));
172	AliITSsegmentation *seg = det->GetSegmentationModel();
173	tot[i] = gm->GetNladders(i+1) * gm->GetNdetectors(i+1) * seg->GetNPads();
174	cout.setf(ios::fixed);
175	cout << tot[i] << endl;
176	Text_t name[7];
177	for (Int_t j = 0; j < 4; j++) {
178	name[0] = 'h';
179	name[1] = 'i';
180	name[2] = 's';
181	name[3] = 't';
182	name[4] = '0' + j;
183	name[5] = '0' + i;
184	name[6] = '\0';
185	Int_t nbins = (Int_t)(2. * zmax[i] / binw[i]);
186	hist[j][i] = new TH1D(name, "histo", nbins, -zmax[i], zmax[i]);
187	}
188	}
189	// Here is put the value of the density for each bin of each layer's TH1D
190	for (Int_t i = 0; i < 6; i++) {
191	for (Int_t j = 0; j < hist[0][i]->GetNbinsX(); j++) {
192	Double_t z = hist[0][i]->GetBinLowEdge(j);
193	Double_t bincontent = tot[i] / (zmax[i] * 2.0) * hist[0][i]->GetBinWidth(j);
194	hist[0][i]->Fill(z + 0.1, bincontent / 100.0);
195	z += hist[0][i]->GetBinWidth(j);
196	}
197	}
198
199	// --------------------------------------------------------------
200	// Counting the hits, digits and recpoints contained in the ITS
201	// --------------------------------------------------------------
202
203	for (Int_t L = 0; L < 6; L++) {
204
205	// To avoid two nested loops, are calculated
206	// the ID of the first and last module of the L
207	// (remember the L goest from 0 to 5 (not from 1 to 6)
208	Int_t first, last;
209	first = gm->GetModuleIndex(L + 1, 1, 1);
210	last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
211
212	// Start loop on modules
213	cout << "Examinating layer " << L + 1 << " ... " << flush;
214	for (Int_t ID = first; ID <= last; ID++) {
215
216	// Hits analysis (if requested)
217	if (HITS)
218	GetModuleHits(ITS, ID, hist[2][L]);
219
220	// Digits analysis (if requested)
221	if (DIGS)
222	GetModuleDigits(ITS, ID, DType(L), hist[1][L]);
223
224	// Recpoints analysis (if requested)
225	if (RECS)
226	GetModuleRecPoints(ITS, ID, hist[3][L]);
227	}
228
229	hist[1][L]->Divide(hist[0][L]);
230	hist[2][L]->Divide(hist[0][L]);
231	hist[3][L]->Divide(hist[0][L]);
232	cout << "DONE" << endl;
233	}
234
235	// Write on the screen the total means
236	cout << endl;
237	cout << "******* MEAN OCCUPANCIES *******" << endl;
238	for (Int_t L = 0; L < 6; L++) {
239	Double_t mean[3];
240	cout << " LAYER " << L << ": " << endl;
241	for (Int_t i = 0; i < 3; i++) {
242	mean[i] = hist[i+1][L]->GetSumOfWeights() / hist[i+1][L]->GetNbinsX();
243	Text_t title[50];
244	sprintf(title, " - Layer %d, mean %4.2f - ", L + 1, mean[i]);
245	hist[i+1][L]->SetTitle(title);
246	cout.setf(ios::fixed);
247	cout.precision(3);
248	if (HITS && i == 1) {
249	cout << " hits occupancy = ";
250	if (mean[i] < 10.0) cout << ' ';
251	cout << mean[i] << "%" << endl;
252	}
253	else if (DIGS && i == 0) {
254	cout << " digits occupancy = ";
255	if (mean[i] < 10.0) cout << ' ';
256	cout << mean[i] << "%" << endl;
257	}
258	if (RECS && i == 2) {
259	cout << " recpts occupancy = ";
260	if (mean[i] < 10.0) cout << ' ';
261	cout << mean[i] << "%" << endl;
262	}
263	}
264	cout << "------------------------------------" << endl;
265	}
266
267	// ----------------------
268	// Plots (if requested)
269	// ----------------------
270	if (!PLOT) return 1;
271
272	TCanvas *view[4];
273	if (HITS) {
274	view[2] = new TCanvas("viewH", "Occupancy view (HITS)", 0, 0, 1050, 700);
275	view[2]->Divide(3,2, 0.001, 0.001);
276	}
277	if (DIGS) {
278	view[1] = new TCanvas("viewD", "Occupancy view (DIGITS)", 20, 40, 1050, 700);
279	view[1]->Divide(3,2, 0.001, 0.001);
280	}
281	if (RECS) {
282	view[3] = new TCanvas("viewR", "Occupancy view (RECPOINTS)", 40, 60, 1050, 700);
283	view[3]->Divide(3,2, 0.001, 0.001);
284	}
285
286	for (Int_t L = 0; L < 6; L++) {
287	for (Int_t i = 1; i <= 3; i++) {
288	if (DIGS && i == 1)
289	hist[i][L]->SetFillColor(kBlue);
290	else if (HITS && i == 2)
291	hist[i][L]->SetFillColor(kRed);
292	else if (RECS && i == 3)
293	hist[i][L]->SetFillColor(kGreen);
294
295	if ((HITS && i == 2) \|\| (DIGS && i == 1) \|\| (RECS && i == 3)) {
296	view[i]->cd(L+1);
297	hist[i][L]->SetStats(kFALSE);
298	hist[i][L]->Draw();
299	hist[i][L]->GetXaxis()->SetTitle("zeta");
300	hist[i][L]->GetYaxis()->SetTitle("%");
301	view[i]->Update();
302	}
303	}
304	}
305
306	if (HITS) view[2]->SaveAs(gifH);
307	if (DIGS) view[1]->SaveAs(gifD);
308	if (RECS) view[3]->SaveAs(gifR);
309
310	return 1;
311	}
312
313	void GetModuleDigits(TObject its, Int_t ID, Int_t dtype, TH1D counter) {
314	// First of all, the macro selects the specified module,
315	// then gets the array of recpoints in it and their number.
316	AliITS ITS = (AliITS)its;
317	TTree *TD = gAlice->TreeD();
318	ITS->ResetDigits();
319	TD->GetEvent(ID);
320	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
321	AliITSgeom *gm = ITS->GetITSgeom();
322	AliITSDetType *det = ITS->DetType(dtype);
323	AliITSsegmentation *seg = det->GetSegmentationModel();
324	Int_t digits_num = digits_array->GetEntries();
325
326	// Get the coordinates of the module
327	for (Int_t j = 0; j < digits_num; j++) {
328	Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
329	AliITSdigit digit = (AliITSdigit)digits_array->UncheckedAt(j);
330	Int_t iz=digit->fCoord1; // cell number z
331	Int_t ix=digit->fCoord2; // cell number x
332	// Get local coordinates of the element (microns)
333	seg->GetPadCxz(ix, iz, lx[0], lx[2]);
334	if (dtype != 1) {
335	// !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
336	if (!dtype) {
337	lx[0] -= seg->Dx() / 2.0;
338	lx[2] -= seg->Dz() / 2.0;
339	}
340	lx[0] /= 10000.0; // convert from microns to cm (if not SDD)
341	lx[2] /= 10000.0; // convert from microns to cm (if not SDD)
342	}
343	gm->LtoG(ID, lx, gx);
344	counter->Fill(gx[2]);
345	}
346	}
347
348	/void GetModuleDigits(TObject its, Int_t ID, Int_t dtype, TH1D *counter) {
349	// First of all, the macro selects the specified module,
350	// then gets the array of recpoints in it and their number.
351	AliITS ITS = (AliITS)its;
352	TTree *TD = gAlice->TreeD();
353	ITS->ResetDigits();
354	TD->GetEvent(ID);
355	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
356	AliITSgeom *gm = ITS->GetITSgeom();
357	AliITSDetType *det = ITS->DetType(dtype);
358	AliITSsegmentation *seg = det->GetSegmentationModel();
359	TArrayI ssdone(5000); // used to match p and n side digits of strips
360	TArrayI pair(5000); // as above
361	Int_t digits_num = digits_array->GetEntries();
362
363	// Get the coordinates of the module
364	if (dtype == 2) {
365	for (Int_t j = 0; j < digits_num; j++) {
366	ssdone[j] = 0;
367	pair[j] = 0;
368	}
369	}
370	for (Int_t j = 0; j < digits_num; j++) {
371	Double_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
372	AliITSdigit digit = (AliITSdigit)digits_array->UncheckedAt(j);
373	Int_t iz=digit->fCoord1; // cell number z
374	Int_t ix=digit->fCoord2; // cell number x
375	// Get local coordinates of the element (microns)
376	if(dtype < 2) {
377	seg->GetPadCxz(ix, iz, lx[0], lx[2]);
378	if (dtype == 0) {
379	// !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
380	lx[0] -= seg->Dx() / 2.0;
381	lx[2] -= seg->Dz() / 2.0;
382	lx[0] /= 10000.0; // convert from microns to cm (SPD)
383	lx[2] /= 10000.0; // convert from microns to cm (SPD)
384	}
385	gm->LtoG(ID, lx, gx);
386	counter->Fill(gx[2]);
387	}
388	else {
389	// SSD: if iz==0 ---> N side; if iz==1 P side
390	if (ssdone[j] == 0) {
391	ssdone[j]=1;
392	pair[j]=-1;
393	Bool_t pside = (iz == 1);
394	Bool_t impaired = kTRUE;
395	Int_t pstrip = 0;
396	Int_t nstrip = 0;
397	if (pside) pstrip = ix; else nstrip = ix;
398	for (Int_t k = 0; k < digits_num; k++) {
399	if (ssdone[k] == 0 && impaired) {
400	AliITSdigitSSD sdigit=(AliITSdigitSSD)digits_array->UncheckedAt(k);
401	if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
402	ssdone[k]=2;
403	pair[j]=k;
404	if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
405	impaired=kFALSE;
406	}
407	}
408	}
409	if (!impaired) {
410	seg->GetPadCxz(pstrip, nstrip, lx[0], lx[2]);
411	lx[0] /= 10000.0;
412	lx[2] /= 10000.0;
413	gm->LtoG(ID, lx, gx);
414	counter->Fill(gx[2]);
415	}
416	}
417	}
418	}
419	}*/
420
421	void GetModuleHits (TObject* its, Int_t ID, TH1D *counter) {
422	// First of all, the macro selects the specified module,
423	// then gets the array of hits in it and their number.
424	AliITS ITS = (AliITS) its;
425	AliITSmodule *module = ITS->GetModule(ID);
426	TObjArray *hits_array = module->GetHits();
427	Int_t hits_num = hits_array->GetEntriesFast();
428	// next, fills the counters with the hits coordinates' calcula
429	for (Int_t j = 0; j < hits_num; j++) {
430	AliITShit hit = (AliITShit) hits_array->At(j);
431	Double_t Z = hit->GetZG();
432	if (!hit->StatusEntering()) counter->Fill(Z);
433	}
434	}
435
436	void GetModuleRecPoints (TObject its, Int_t ID, TH1D counter) {
437	// First of all, the macro selects the specified module,
438	// then gets the array of recpoints in it and their number.
439	AliITS ITS = (AliITS) its;
440	AliITSgeom *gm = ITS->GetITSgeom();
441	TTree *TR = gAlice->TreeR();
442	if (!TR \|\| TR->GetEntries() == 0) {
443	cout << "The ITS wasn't clustered..." << endl;
444	return;
445	}
446	ITS->ResetRecPoints();
447	TR->GetEvent(ID);
448	TClonesArray *recs_array = ITS->RecPoints();
449	Int_t recs_num = recs_array->GetEntries();
450
451	for(Int_t j = 0; j < recs_num; j++) {
452	Double_t lx[3] = {0.,0.,0.}, gx[3] = {0.,0.,0.};
453	AliITSRecPoint recp = (AliITSRecPoint)recs_array->At(j);
454	lx[0] = recp->GetX();
455	lx[1] = 0.0;
456	lx[2] = recp->GetZ();
457	gm->LtoG(ID, lx, gx);
458	counter->Fill(gx[2]);
459	}
460	}
461
462	Int_t DType(Int_t layer) {
463	if (layer == 0 \|\| layer == 1)
464	return 0;
465	else if (layer == 2 \|\| layer == 3)
466	return 1;
467	else
468	return 2;
469	}