[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 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
af8e1c2d	1	/******************************************************************************
	2
	3	"AliITSOccupancy.C"
	4
	5	this macro calculates the mean occupancy of each ITS layer, counting the
b81b287f	6	"fired" digits of each module, and making two overall calculations:
af8e1c2d	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
b81b287f	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).
af8e1c2d	16
b81b287f	17
af8e1c2d	18	It is also possible (and recommended) to compile this macro,
af8e1c2d	19	to execute it faster. To do this, it is necessary to:
b81b287f	20
b81b287f	21	1) type, at the root prompt, the instruction
af8e1c2d	22	gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
af8e1c2d	23	to adjoin the ALICE include libraries to the main include directory
b81b287f	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()
af8e1c2d	30
	31	Author: Alberto Pulvirenti
	32	******************************************************************************/
	33
af8e1c2d	34
b81b287f	35	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
af8e1c2d	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
b81b287f	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
af8e1c2d	83
b81b287f	84	// Histos for plotting occupancy distributions
b81b287f	85	TH1F above[6], below[6];
af8e1c2d	86
b81b287f	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	}
af8e1c2d	116
b81b287f	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];
af8e1c2d	122
b81b287f	123	for (Int_t L = 0; L < 6; L++) {
af8e1c2d	124
b81b287f	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));
af8e1c2d	133
b81b287f	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);
af8e1c2d	140
b81b287f	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	}
af8e1c2d	156
b81b287f	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	}
af8e1c2d	168
b81b287f	169	TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
b81b287f	170	view->Divide(2, 3);
af8e1c2d	171
b81b287f	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;
af8e1c2d	223	}