ITS/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", Bool_t recreate = kFALSE, 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[]  = {20,20,20,28,22,24};             // bins number for TH1Fs
 127
 128         // Histos for plotting occupancy distributions
 129         TH1F *above[6], *below[6];
 130
 131         if (recreate) {
 132                 // This work is done if the "recreate" option is true
 133                 TFile *file = new TFile("totals.root", "RECREATE");
 134                 file->cd();
 135                 for (Int_t lay = 0; lay < 6; lay++) {
 136                         Int_t nlads = gm->GetNladders(lay+1);
 137                         Int_t ndets = gm->GetNdetectors(lay+1);
 138                         Int_t dtype = lay / 2;
 139                         Int_t minID = gm->GetModuleIndex(lay+1, 1, 1);
 140                         Int_t maxID = gm->GetModuleIndex(lay+1, nlads, ndets);
 141                         Text_t ffname[20];
 142                         sprintf(ffname, "h_%d", lay+1);
 143                         below[lay] = new TH1F(ffname, "Total z distribution of digits", nbins[lay], -zmax[lay], zmax[lay]);
 144                         cout << "Evaluating total channels number of layer " << lay+1 << endl;
 145                         for (Int_t I = minID; I <= maxID; I++) {
 146                                 AliITSDetType *det = ITS->DetType(dtype);
 147                                 AliITSsegmentation *seg = det->GetSegmentationModel();
 148                                 Int_t NX = seg->Npx();
 149                                 Int_t NZ = seg->Npz();
 150                                 for (Int_t ix = 0; ix <= NX; ix++) {
 151                                         for (Int_t iz = 0; iz <= NZ; iz++) {
 152                                                 Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
 153                                                 seg->DetToLocal(ix, iz, lx[0], lx[2]);
 154                                                 gm->LtoG(I, lx, gx);
 155                                                 below[lay]->Fill(gx[2]);
 156                                         }
 157                                 }
 158                         }
 159                         // every generated histogram is written to the file
 160                         below[lay]->Write();
 161                         // and every counter histogram is created
 162                         sprintf(ffname, "H_%d", lay+1);
 163                         above[lay] = new TH1F(ffname, "histo", nbins[lay], -zmax[lay], zmax[lay]);
 164                 }
 165         }
 166         else {
 167                 // Thi other work is done if the recreate option is false
 168                 TFile *file = new TFile("totals.root", "READ");
 169                 file->ls();
 170                 for (Int_t i = 0; i < 6; i++) {
 171                         Text_t ffname[20];
 172                         sprintf(ffname, "h_%d", i+1);
 173                         below[i] = (TH1F*)file->Get(ffname);
 174                         sprintf(ffname, "H_%d", i+1);
 175                         above[i] = new TH1F(ffname, "histo", nbins[i], -zmax[i], zmax[i]);
 176                 }
 177         }
 178
 179         // Counting the hits, digits and recpoints contained in the ITS
 180         TTree *TD = gAlice->TreeD();
 181         ITS->ResetDigits();
 182         Float_t mean[6];
 183
 184         for (Int_t L = 0; L < 6; L++) {
 185
 186                 cout << "Layer " << L + 1 << ": " << flush;
 187
 188                 // To avoid two nested loops, are calculated
 189                 // the ID of the first and last module of the L
 190                 // (remember the L goest from 0 to 5 (not from 1 to 6)
 191                 Int_t first, last;
 192                 first = gm->GetModuleIndex(L + 1, 1, 1);
 193                 last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
 194
 195                 // Start loop on modules
 196                 for (Int_t ID = first; ID <= last; ID++) {
 197                         Int_t dtype = L / 2;
 198                         AliITSDetType *det = ITS->DetType(dtype);
 199                         AliITSsegmentation *seg = det->GetSegmentationModel();
 200                         if (dtype == 2) seg->SetLayer(L+1);
 201
 202                         TD->GetEvent(ID);
 203                         TClonesArray *digits_array = ITS->DigitsAddress(dtype);
 204                         Int_t digits_num = digits_array->GetEntries();
 205                         // Get the coordinates of the module
 206                         for (Int_t j = 0; j < digits_num; j++) {
 207                                 Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
 208                                 AliITSdigit *digit = (AliITSdigit*)digits_array->UncheckedAt(j);
 209                                 Int_t iz=digit->fCoord1;  // cell number z
 210                                 Int_t ix=digit->fCoord2;  // cell number x
 211                         // Get local coordinates of the element (microns)
 212                         seg->DetToLocal(ix, iz, lx[0], lx[2]);
 213                                 gm->LtoG(ID, lx, gx);
 214                                 above[L]->Fill(gx[2]);
 215                         }
 216                 }
 217
 218                 above[L]->Divide(above[L], below[L], 100.0, 1.0);
 219                 mean[L] = above[L]->GetSumOfWeights() / above[L]->GetNbinsX();
 220                 cout.setf(ios::fixed);
 221                 cout.precision(2);
 222                 cout << " digits occupancy = " << mean[L] << "%" << endl;
 223         }
 224
 225         TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
 226         view->Divide(2, 3);
 227
 228         for (Int_t I = 1; I < 7; I++) {
 229                 view->cd(I);
 230                 Text_t title[50];
 231                 sprintf(title, "Layer %d: %4.2f%c", I, mean[I-1], '%');
 232                 title[6] = (Text_t)I + '0';
 233                 above[I-1]->SetTitle(title);
 234                 above[I-1]->SetStats(kFALSE);
 235                 above[I-1]->SetXTitle("z (cm)");
 236                 above[I-1]->SetYTitle("%");
 237                 above[I-1]->Draw();
 238                 view->Update();
 239         }
 240         Text_t *filegif  = new Text_t[filelen + 11];
 241         Text_t *fileps  = new Text_t[filelen + 11];
 242         strcpy (filegif, name);
 243         strcpy (fileps, name);
 244         strcat (filegif, "_digs.gif");
 245         strcat (fileps, "_digs.ps");
 246         view->SaveAs(filegif);
 247         view->SaveAs(fileps);
 248
 249
 250         return 1;
 251 }