ITS/ITSOccupancy.C

   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 }