ITS/ITSReadPlotData.C

   1 Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
   2
   3         /*********************************************************************
   4          *                                                                   *
   5          *  Macro used to read hits, digits and recpoints for a module       *
   6          *  It draws 3 TH2Fs where stores the 2-D coordinates of these       *
   7          *  data for a specified module (for which the macro asks when       *
   8          *  it starts, with some checks to avoid wrong detector coords.      *
   9          *                                                                   *
  10          *  Only a little 'experimental' warning:                            *
  11          *  with all the tests I have made, I wasn't able to plot the        *
  12          *  digits fo the 5th layer...                                       *
  13          *  I skipped this problem with an alert to beware th user, while    *
  14          *  in this situation the macro plots only recpoints and hits        *
  15          *                                                                   *
  16          *  Author: Alberto Pulvirenti                                       *
  17          *                                                                   *
  18          *********************************************************************/
  19
  20
  21         extern Int_t GetModuleHits (TObject* its, Int_t ID[0], Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St);
  22         extern Int_t GetModuleRecPoints (TObject *its, Int_t ID[0], Float_t*& X, Float_t*& Z);
  23         extern Int_t GetModuleDigits(TObject *its, Int_t ID[0], Int_t dtype, Float_t*& X, Float_t*& Z);
  24         extern void  AssignCoords(TArrayI *ID);
  25
  26         // First of all, here are put some variable declarations
  27         // that are useful in the following part:
  28         Int_t nparticles; // number of particles
  29         // ITS module coordinates [layer = 1, ladder = 2, det = 3] and absolute ID[0] of module [0]
  30         TArrayI ID(4);
  31         Int_t nmodules, dtype; // Total number of modules and module type (SSD, SPD, SDD)
  32         Float_t *x = 0, *y = 0, *z = 0; // Arrays where to store read coords
  33         Bool_t *St = 0; // Status of the track (hits only)
  34
  35         // It's necessary to load the gAlice shared libs
  36         // if they aren't already stored in memory...
  37         if (gClassTable->GetID("AliRun") < 0) {
  38         gROOT->LoadMacro("loadlibs.C");
  39                 loadlibs();
  40         }
  41   // Anyway, this macro needs to read a gAlice file, so it
  42   // clears the gAlice object if there is already one in memory...
  43   else {
  44                 if(gAlice){
  45                         delete gAlice;
  46                         gAlice = 0;
  47                 }
  48         }
  49
  50         // Now is opened the Root input file containing Geometry, Kine and Hits
  51   // by default its name must be "galice.root".
  52   // When the file is opened, its contens are shown.
  53         TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  54         if (!file) file = new TFile(filename);
  55                 file->ls();
  56
  57         // Then, the macro gets the AliRun object from file.
  58         // If this object is not present, an error occurs
  59         // and the execution is stopped.
  60         // Anyway, this operation needs some time,
  61         // don't worry about an apparent absence of output and actions...
  62         cout << "\nSearching in '" << filename << "' for an AliRun object ... " << flush;
  63         gAlice = (AliRun*)file->Get("gAlice");
  64         if (gAlice)
  65                 cout << "FOUND!" << endl;
  66         else {
  67                 cout<<"NOT FOUND! The Macro can't continue!" << endl;
  68                 return 0;
  69         }
  70
  71         // Then, the macro selects the event number specified. Default is 0.
  72         nparticles = gAlice->GetEvent(evNum);
  73         cout << "\nNumber of particles   = " << nparticles << endl;
  74         if (!nparticles) {
  75                 cout << "With no particles I can't do much... Goodbye!" << endl;
  76                 return 0;
  77         }
  78
  79         // The next step is filling the ITS from the AliRun object.
  80         AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
  81   ITS->InitModules(-1, nmodules);
  82   cout << "Number of ITS modules = " << nmodules << endl;
  83         cout << "\nFilling modules (it takes a while, now)..." << flush;
  84         ITS->FillModules(0, 0, nmodules, " ", " ");
  85   cout << "DONE!" << endl;
  86         AliITSgeom *gm = ITS->GetITSgeom();
  87         AliITSDetType *det = ITS->DetType(dtype);
  88         AliITSsegmentation *seg = det->GetSegmentationModel();
  89
  90         for(;;) {
  91
  92     // Input phase.
  93     // The macro asks if the user wants to put a detector ID[0]
  94     // or prefers to input layer, ladder and detector.
  95     for (Int_t i = 0; i < 4; i++) ID[i] = -1;
  96     Int_t answ;
  97     do {
  98            cout << "\nSelection modes:" << endl;
  99                 cout << "1) by layer - ladder - detector numbers" << endl;
 100         cout << "2) by unique detector ID" << endl;
 101            cout << "0) exit macro" << endl;
 102                 cout << "\nEnter your choice: ";
 103         cin >> answ;
 104     } while (answ < 0 || answ > 2);
 105     switch (answ) {
 106         case 0:
 107                 // input = 0 ---> EXIT
 108                 return;
 109                 break;
 110         case 1:
 111                 // input = 1 ---> SELECTION BY COORD
 112                 do {
 113                                         cout << "\nLayer number [1-6, 0 to exit]: ";
 114                                         cin >> ID[1];
 115                                         if (!ID[1]) return 0;
 116                                 } while (ID[1] < 0 || ID[1] > 6);
 117
 118                                 // Detector type: 0 --> SPD, 1 --> SDD, 2 --> SSD.
 119                                 // Layer 1,2 --> 0 / Layer 3,4 --> 1 / Layer 5,6 --> 2
 120                                 dtype = ID[1] / 3;
 121
 122                                 // Once fixed the layer number, the macro calculates the max number
 123                                 // for ladder and detector from geometry, and accepts only suitable values.
 124                                 do {
 125                                         ID[2] = gm->GetNladders(ID[1]);
 126                                         cout << "Ladder number [1-" << ID[2] << ", 0 to exit]: ";
 127                                         cin >> ID[2];
 128                                         if (!ID[2]) return 0;
 129                                 } while (ID[2] < 0 || ID[2] > gm->GetNladders(ID[1]));
 130                                 do {
 131                                         ID[3] = gm->GetNdetectors(ID[1]);
 132                                         cout << "Detector number [1-" << ID[3] << ", 0 to exit]: ";
 133                                         cin >> ID[3];
 134                                         if (!ID[3]) return 0;
 135                                 } while (ID[3] < 0 || ID[3] > gm->GetNdetectors(ID[1]));
 136                                 break;
 137         case 2:
 138                 // input = 2 ---> SELECTION BY ID[0]
 139                 do {
 140                                         ID[0] = gm->GetIndexMax();
 141                                         cout << "\n Detector ID number [0-" << ID[0] << ", -1 to exit]: ";
 142                                         cin >> ID[0];
 143                                         if (ID[0] == -1) return 0;
 144                                 } while (ID[0] < 0 || ID[0] > gm->GetIndexMax());
 145           break;
 146     };
 147
 148     if (ID[0] == -1)
 149                 // If ID[0] = -1 the chioce was by coords, so it's necessary to assign the ID:
 150                         ID[0] = gm->GetModuleIndex(ID[1], ID[2], ID[3]);
 151                 else {
 152                         // Else we must get the layer, ladder and detector by the ID;
 153                         // ...but first we must remember that the ID goest from 0 to NModules - 1!
 154                         ID[0]--;
 155                         ID[1] = ITS->GetModule(ID[0])->GetLayer();
 156                         ID[2] = ITS->GetModule(ID[0])->GetLadder();
 157                         ID[3] = ITS->GetModule(ID[0])->GetDet();
 158                 }
 159
 160                 // Defines the histograms inside the `for' cycle, so they are destroyed at the end
 161                 // of every read sequqnce, in order to mek another withour segmentation faults
 162                 Text_t msg[250], xm = 0.0, ym = 0.0;
 163                 switch (dtype) {
 164                         case 0: xm = 1.5; zm = 7.0; break;
 165                         case 1: xm = 7.5; zm = 8.0; break;
 166                         case 2: xm = 7.5; zm = 4.5; break;
 167                 }
 168                 sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]);
 169                 TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of hits
 170                 TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of recpoints
 171                 TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits
 172
 173                 cout << endl;
 174
 175                 // Reads hits...
 176                 Int_t hits = GetModuleHits(ITS, ID[0], x, y, z, St);
 177                 if (!hits) {
 178                         cout << "No hits in module!" << endl;
 179                         continue;
 180                 }
 181                 for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
 182
 183                 // Reads recpoints...
 184                 Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z);
 185                 if (!recs) {
 186                         cout << "No recpoints in module!" << endl;
 187                         continue;
 188                 }
 189                 for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
 190
 191                 // Reads digits...
 192                 Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
 193                 if (!digits) {
 194                         cout << "No digits in module!" << endl;
 195                         //continue;
 196                 }
 197                 for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
 198
 199                 // Draws read data...
 200                 // HITS -------> red (2) crosses.
 201                 // DIGITS -----> green (8) boxes.
 202                 // REC-POINTS -> blue (4) St. Andrew's crosses.
 203
 204                 TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
 205                 viewer->cd();
 206
 207                 hdigits->SetMarkerStyle(25);
 208                 hdigits->SetMarkerColor(8);
 209                 hdigits->SetMarkerSize(2);
 210                 hdigits->SetStats(kFALSE);
 211                 hdigits->SetXTitle("Local X (cm)");
 212                 hdigits->SetYTitle("Local Z (cm)");
 213                 hdigits->Draw();
 214
 215                 hhits->SetMarkerStyle(5);
 216                 hhits->SetMarkerColor(2);
 217                 hhits->SetMarkerSize(3);
 218                 hhits->SetStats(kFALSE);
 219                 hhits->Draw("same");
 220
 221                 hrecs->SetMarkerStyle(2);
 222                 hrecs->SetMarkerColor(4);
 223                 hrecs->SetMarkerSize(3);
 224                 hrecs->SetStats(kFALSE);
 225                 hrecs->Draw("same");
 226
 227                 TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
 228                 legend->SetMargin(0.2);
 229                 legend->AddEntry(hhits, "hits","P");
 230                 legend->AddEntry(hrecs, "recpoints","P");
 231                 legend->AddEntry(hdigits, "digits","P");
 232                 legend->SetTextSizePixels(14);
 233                 legend->Draw();
 234
 235                 viewer->Update();
 236         }
 237
 238         cout << "Done. Goodbye" << endl;
 239         return;
 240 }
 241
 242 Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {
 243         // First of all, the macro selects the specified module,
 244         // then gets the array of hits in it and their number.
 245         AliITS *ITS = (AliITS*) its;
 246         AliITSmodule *module = ITS->GetModule(ID);
 247         TObjArray *hits_array = module->GetHits();
 248         Int_t hits_num = hits_array->GetEntriesFast();
 249
 250         // Now, if this count returns 0, there's nothing to do,
 251         // while, if it doesn't, the first thing to do is dimensioning
 252         // the coordinate arrays, and then the loop can start.
 253         if (!hits_num)
 254                 return 0;
 255         else {
 256                 if (X) delete [] X;
 257                 if (Y) delete [] Y;
 258                 if (Z) delete [] Z;
 259                 if (St) delete [] St;
 260                 X = new Float_t[hits_num];
 261                 Y = new Float_t[hits_num];
 262                 Z = new Float_t[hits_num];
 263                 St = new Int_t[hits_num];
 264         }
 265
 266         for (Int_t j = 0; j < hits_num; j++) {
 267                 AliITShit *hit = (AliITShit*) hits_array->At(j);
 268                 X[j]  = hit->GetXL();
 269                 Y[j]  = hit->GetYL();
 270                 Z[j]  = hit->GetZL();
 271                 St[j] = hit->StatusEntering();
 272         }
 273         return hits_num;
 274 }
 275
 276 Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z) {
 277
 278         // First of all, the macro selects the specified module,
 279         // then gets the array of recpoints in it and their number.
 280         AliITS *ITS = (AliITS*) its;
 281         TTree *TR = gAlice->TreeR();
 282         ITS->ResetRecPoints();
 283         TR->GetEvent(ID);
 284         TClonesArray *recs_array = ITS->RecPoints();
 285         Int_t recs_num = recs_array->GetEntries();
 286
 287         // Now, if this count returns 0, there's nothing to do,
 288         // while, if it doesn't, the first thing to do is dimensioning
 289         // the coordinate and energy loss arrays, and then the loop can start.
 290         if (!recs_num)
 291                 return 0;
 292         else {
 293                 if (X) delete [] X;
 294                 if (Z) delete [] Z;
 295                 X = new Float_t[recs_num];
 296                 Z = new Float_t[recs_num];
 297         }
 298         for(Int_t j = 0; j < recs_num; j++) {
 299                 AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
 300            X[j] = recp->GetX();
 301                 Z[j] = recp->GetZ();
 302         }
 303         return recs_num;
 304 }
 305
 306 Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {
 307
 308         // First of all, the macro selects the specified module,
 309         // then gets the array of recpoints in it and their number.
 310         AliITS *ITS = (AliITS*)its;
 311         TTree *TD = gAlice->TreeD();
 312         ITS->ResetDigits();
 313         TD->GetEvent(ID);
 314         TClonesArray *digits_array = ITS->DigitsAddress(dtype);
 315         AliITSgeom *gm = ITS->GetITSgeom();
 316         AliITSDetType *det = ITS->DetType(dtype);
 317         AliITSsegmentation *seg = det->GetSegmentationModel();
 318         TArrayI ssdone(5000);  // used to match p and n side digits of strips
 319         TArrayI pair(5000);    // as above
 320         Int_t digits_num = digits_array->GetEntries();
 321         // Now, if this count returns 0, there's nothing to do,
 322         // while, if it doesn't, the first thing to do is dimensioning
 323         // the coordinate and energy loss arrays, and then the loop can start.
 324
 325         if (!digits_num)
 326                 return 0;
 327         else {
 328                 if (X) delete [] X;
 329                 if (Z) delete [] Z;
 330                 X = new Float_t[digits_num];
 331                 Z = new Float_t[digits_num];
 332         }
 333
 334         // Get the coordinates of the module
 335         if (dtype == 2) {
 336                 for (Int_t j = 0; j < digits_num; j++) {
 337                         ssdone[j] = 0;
 338                         pair[j] = 0;
 339                 }
 340         }
 341   for (Int_t j = 0; j < digits_num; j++) {
 342         cout << j << endl;
 343                 digit = (AliITSdigit*)digits_array->UncheckedAt(j);
 344                 Int_t iz=digit->fCoord1;  // cell number z
 345                 Int_t ix=digit->fCoord2;  // cell number x
 346     // Get local coordinates of the element (microns)
 347                 if(dtype < 2)
 348         seg->GetPadCxz(ix, iz, X[j], Z[j]);
 349     else {
 350                         // SSD: if iz==0 ---> N side; if iz==1 P side
 351       if (ssdone[j] == 0) {
 352                                 ssdone[j]=1;
 353                                 pair[j]=-1;
 354                                 Bool_t pside = (iz == 1);
 355                                 Bool_t impaired = kTRUE;
 356                                 Int_t pstrip = 0;
 357                                 Int_t nstrip = 0;
 358                                 if (pside) pstrip = ix; else nstrip = ix;
 359                                 for (Int_t k = 0; k < digits_num; k++) {
 360                                         if (ssdone[k] == 0 && impaired) {
 361                                                 AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
 362                                                 if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
 363                                                         ssdone[k]=2;
 364                                                         pair[j]=k;
 365                                                         if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
 366                                                         impaired=kFALSE;
 367                                                 }
 368                                         }
 369                                 }
 370         if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
 371                         }
 372                 }
 373                 if (dtype == 0) {
 374                         // !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
 375                         X[j] = X[j]-seg->Dx() / 2.0;
 376                         Z[j] = Z[j]-seg->Dz() / 2.0;
 377                 }
 378                 if (dtype != 1) {
 379                         X[j] /= 10000.0;
 380                         Z[j] /= 10000.0;
 381                 }
 382         }
 383         return digits_num;
 384 }
 385