[u/mrichter/AliRoot.git] / ITS / oldmacros / 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 calculations:
     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 and the EPS of the TCanvas where the 
  histograms are plotted. Muliple input files are supported (see the 
  input argument list below).
			  
                                                                        
  It is also possible (and recommended) to compile this macro, 
  to execute it faster. To do this, it is necessary to:
 
  1) 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
  2) load the function via the ".L AliITSOccupancy.C++" statement 
     (only the first time, because the first time the macro must be compiled). 
      From 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
  3) call the function with only its name, e.g.: AliITSOccupancy()

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


#if !defined(__CINT__) || defined(__MAKECINT__)
	#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

TFile* AccessFile(TString inFile="galice.root", TString acctype="R");

Int_t AliITSOccupancy(TString FileHits="galice.root", TString FileDigits="galice.root", Int_t evNum = 0) {       
		 	
  // Open the Root input file containing the AliRun data and 
  // the Root output data file
  TFile *froot = AccessFile(FileHits);
  froot->ls();

  if(!(FileDigits.Data() == FileHits.Data())){
    gAlice->SetTreeDFileName(FileDigits);
  }
  gAlice->GetEvent(evNum);
  // 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();
  }

  view->SaveAs("AliITSOccupancy_digit.gif");
  view->SaveAs("AliITSOccupancy_digit.eps");
  TFile *fOc = new TFile("AliITSOccupancy.root","recreate");
  fOc->cd();
  for (Int_t I = 0; I < 6; I++) above[I]->Write();
  fOc->Close();

  return 1;
}

//-------------------------------------------------------------------
TFile * AccessFile(TString FileName, TString acctype){

  // Function used to open the input file and fetch the AliRun object

  TFile *retfil = 0;
  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(FileName);
  if (file) {file->Close(); delete file; file = 0;}
  if(acctype.Contains("U")){
    file = new TFile(FileName,"update");
  }
  if(acctype.Contains("N") && !file){
    file = new TFile(FileName,"recreate");
  }
  if(!file) file = new TFile(FileName);   // default readonly
  if (!file->IsOpen()) {
	cerr<<"Can't open "<<FileName<<" !" << endl;
	return retfil;
  } 

  // Get AliRun object from file or return if not on file
  if (gAlice) {delete gAlice; gAlice = 0;}
  gAlice = (AliRun*)file->Get("gAlice");
  if (!gAlice) {
	cerr << "AliRun object not found on file"<< endl;
	return retfil;
  } 
  return file;
}
Commit	Line	Data
b3d9d240	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 calculations:
	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 and the EPS of the TCanvas where the
	14	histograms are plotted. Muliple input files are supported (see the
	15	input argument list below).
	16
	17
	18	It is also possible (and recommended) to compile this macro,
	19	to execute it faster. To do this, it is necessary to:
	20
	21	1) type, at the root prompt, the instruction
	22	gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
	23	to adjoin the ALICE include libraries to the main include directory
	24	2) load the function via the ".L AliITSOccupancy.C++" statement
	25	(only the first time, because the first time the macro must be compiled).
	26	From the second time you use the macro, you must only load it via the
	27	".L AliITSOccupancy.C+" instruction (note the number of '+' signs in
	28	each case
	29	3) call the function with only its name, e.g.: AliITSOccupancy()
	30
	31	Author: Alberto Pulvirenti
	32	******************************************************************************/
	33
	34
	35	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	36	#include <iostream.h>
	37	#include <TROOT.h>
	38	#include <TArrayI.h>
	39	#include <TCanvas.h>
	40	#include <TClassTable.h>
	41	#include <TClonesArray.h>
	42	#include <TFile.h>
	43	#include <TObject.h>
	44	#include <TObjArray.h>
	45	#include <TTree.h>
	46	#include <TMath.h>
	47	#include <TString.h>
	48	#include <TH1.h>
	49	#include <AliRun.h>
	50	#include <AliITS.h>
	51	#include <AliITSgeom.h>
	52	#include <AliITSDetType.h>
	53	#include <AliITSRecPoint.h>
	54	#include <AliITSdigit.h>
	55	#include <AliITShit.h>
	56	#include <AliITSmodule.h>
	57	#include <AliITSsegmentation.h>
	58	#include <AliITSsegmentationSPD.h>
	59	#include <AliITSsegmentationSDD.h>
	60	#include <AliITSsegmentationSSD.h>
	61	#endif
	62
	63	TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
	64
65	Int_t AliITSOccupancy(TString FileHits="galice.root", TString FileDigits="galice.root", Int_t evNum = 0) {
66
67	// Open the Root input file containing the AliRun data and
68	// the Root output data file
69	TFile *froot = AccessFile(FileHits);
70	froot->ls();
71
72	if(!(FileDigits.Data() == FileHits.Data())){
73	gAlice->SetTreeDFileName(FileDigits);
74	}
75	gAlice->GetEvent(evNum);
76	// Initialize the ITS and its geometry
77	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
78	AliITSgeom *gm = ITS->GetITSgeom();
79
80	// Variables and objects definition and setting
81	Float_t zmax[] = {20.0,20.0,25.0,35.0,49.5,54.0}; // z edges for TH1F (cm)
82	Int_t nbins[] = {40,40,50,70,99,108}; // bins number for TH1Fs
83
84	// Histos for plotting occupancy distributions
85	TH1F above[6], below[6];
86
87	for (Int_t lay = 0; lay < 6; lay++) {
88	Int_t nlads = gm->GetNladders(lay+1);
89	Int_t ndets = gm->GetNdetectors(lay+1);
90	Int_t dtype = lay / 2;
91	Int_t minID = gm->GetModuleIndex(lay+1, 1, 1);
92	Int_t maxID = gm->GetModuleIndex(lay+1, nlads, ndets);
93	Text_t ffname[20];
94	sprintf(ffname, "h_%d", lay+1);
95	below[lay] = new TH1F(ffname, "Total z distribution of digits", nbins[lay], -zmax[lay], zmax[lay]);
96	cout << "Evaluating total channels number of layer " << lay+1 << endl;
97	for (Int_t I = minID; I <= maxID; I++) {
98	AliITSDetType *det = ITS->DetType(dtype);
99	AliITSsegmentation *seg = det->GetSegmentationModel();
100	Int_t NX = seg->Npx();
101	if(lay == 2 \|\| lay == 3) NX = 192;
102	Int_t NZ = seg->Npz();
103	// cout << "ID: " << I << ", Layer: " << lay+1 << ", NX = " << NX << ", NZ = " << NZ << endl;
104	for (Int_t ix = 0; ix <= NX; ix++) {
105	for (Int_t iz = 0; iz <= NZ; iz++) {
106	Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
107	seg->DetToLocal(ix, iz, lx[0], lx[2]);
108	gm->LtoG(I, lx, gx);
109	below[lay]->Fill(gx[2]);
110	}
111	}
112	}
113	sprintf(ffname, "H_%d", lay+1);
114	above[lay] = new TH1F(ffname, "histo", nbins[lay], -zmax[lay], zmax[lay]);
115	}
116
117	// Counting the hits, digits and recpoints contained in the ITS
118	TTree *TD = gAlice->TreeD();
119	ITS->ResetDigits();
120	Float_t mean[6];
121	Float_t average[6];
122
123	for (Int_t L = 0; L < 6; L++) {
124
125	cout << "Layer " << L + 1 << ": " << flush;
126
127	// To avoid two nested loops, are calculated
128	// the ID of the first and last module of the L
129	// (remember the L goest from 0 to 5 (not from 1 to 6)
130	Int_t first, last;
131	first = gm->GetModuleIndex(L + 1, 1, 1);
132	last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
133
134	// Start loop on modules
135	for (Int_t ID = first; ID <= last; ID++) {
136	Int_t dtype = L / 2;
137	AliITSDetType *det = ITS->DetType(dtype);
138	AliITSsegmentation *seg = det->GetSegmentationModel();
139	if (dtype == 2) seg->SetLayer(L+1);
140
141	TD->GetEvent(ID);
142	TClonesArray *digits_array = ITS->DigitsAddress(dtype);
143	Int_t digits_num = digits_array->GetEntries();
144	// Get the coordinates of the module
145	for (Int_t j = 0; j < digits_num; j++) {
146	Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
147	AliITSdigit digit = (AliITSdigit)digits_array->UncheckedAt(j);
148	Int_t iz=digit->fCoord1; // cell number z
149	Int_t ix=digit->fCoord2; // cell number x
150	// Get local coordinates of the element (microns)
151	seg->DetToLocal(ix, iz, lx[0], lx[2]);
152	gm->LtoG(ID, lx, gx);
153	above[L]->Fill(gx[2]);
154	}
155	}
156
157	Float_t occupied = above[L]->GetEntries();
158	Float_t total = below[L]->GetEntries();
159	cout << "Entries: " << occupied << ", " << total << endl;
160	average[L] = 100.*occupied/total;
161	above[L]->Divide(above[L], below[L], 100.0, 1.0);
162	mean[L] = above[L]->GetSumOfWeights() / above[L]->GetNbinsX();
163	cout.setf(ios::fixed);
164	cout.precision(2);
165	cout << " digits occupancy = " << mean[L] << "%" << endl;
166	cout << " average digits occupancy = " << average[L] << "%" << endl;
167	}
168
169	TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
170	view->Divide(2, 3);
171
172	for (Int_t I = 1; I < 7; I++) {
173	view->cd(I);
174	Text_t title[50];
175	sprintf(title, "Layer %d: %4.2f%c", I, mean[I-1], '%');
176	title[6] = (Text_t)I + '0';
177	above[I-1]->SetTitle(title);
178	above[I-1]->SetStats(kFALSE);
179	above[I-1]->SetXTitle("z (cm)");
180	above[I-1]->SetYTitle("%");
181	above[I-1]->Draw();
182	view->Update();
183	}
184
185	view->SaveAs("AliITSOccupancy_digit.gif");
186	view->SaveAs("AliITSOccupancy_digit.eps");
187	TFile *fOc = new TFile("AliITSOccupancy.root","recreate");
188	fOc->cd();
189	for (Int_t I = 0; I < 6; I++) above[I]->Write();
190	fOc->Close();
191
192	return 1;
193	}
194
195	//-------------------------------------------------------------------
196	TFile * AccessFile(TString FileName, TString acctype){
197
198	// Function used to open the input file and fetch the AliRun object
199
200	TFile *retfil = 0;
201	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(FileName);
202	if (file) {file->Close(); delete file; file = 0;}
203	if(acctype.Contains("U")){
204	file = new TFile(FileName,"update");
205	}
206	if(acctype.Contains("N") && !file){
207	file = new TFile(FileName,"recreate");
208	}
209	if(!file) file = new TFile(FileName); // default readonly
210	if (!file->IsOpen()) {
211	cerr<<"Can't open "<<FileName<<" !" << endl;
212	return retfil;
213	}
214
215	// Get AliRun object from file or return if not on file
216	if (gAlice) {delete gAlice; gAlice = 0;}
217	gAlice = (AliRun*)file->Get("gAlice");
218	if (!gAlice) {
219	cerr << "AliRun object not found on file"<< endl;
220	return retfil;
221	}
222	return file;
223	}