ITS/oldmacros/ITSOccupancy.C

   1 /******************************************************************************
   2
   3   "ITSOccupancy.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 ITSOccupancy.C+" instruction (note the number of '+' signs in each case
  41   4) call the function with only its name, es: ITSOccupancy("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 ITSOccupancy(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 }