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] - 1) / 2;
 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];
 163                 Float_t xm = 0.0, ym = 0.0, zm = 0.0;
 164                 switch (dtype) {
 165                         case 0: xm = 1.5; zm = 7.0; break;
 166                         case 1: xm = 7.5; zm = 8.0; break;
 167                         case 2: xm = 7.5; zm = 4.5; break;
 168                 }
 169                 sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]);
 170                 TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of hits
 171                 TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of recpoints
 172                 TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits
 173
 174                 cout << endl;
 175
 176                 // Reads hits...
 177                 Int_t hits = GetModuleHits(ITS, ID[0], x, y, z, St);
 178                 if (!hits) {
 179                         cout << "No hits in module!" << endl;
 180                         continue;
 181                 }
 182                 for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
 183
 184                 // Reads recpoints...
 185                 Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z);
 186                 if (!recs) {
 187                         cout << "No recpoints in module!" << endl;
 188                         continue;
 189                 }
 190                 for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
 191
 192                 // Reads digits...
 193                 Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
 194                 if (!digits) {
 195                         cout << "No digits in module!" << endl;
 196                         //continue;
 197                 }
 198                 for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
 199
 200                 // Draws read data...
 201                 // HITS -------> red (2) crosses.
 202                 // DIGITS -----> green (8) boxes.
 203                 // REC-POINTS -> blue (4) St. Andrew's crosses.
 204
 205                 TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
 206                 viewer->cd();
 207
 208                 hdigits->SetMarkerStyle(25);
 209                 hdigits->SetMarkerColor(8);
 210                 hdigits->SetMarkerSize(2);
 211                 hdigits->SetStats(kFALSE);
 212                 hdigits->SetXTitle("Local X (cm)");
 213                 hdigits->SetYTitle("Local Z (cm)");
 214                 hdigits->Draw();
 215
 216                 hhits->SetMarkerStyle(5);
 217                 hhits->SetMarkerColor(2);
 218                 hhits->SetMarkerSize(3);
 219                 hhits->SetStats(kFALSE);
 220                 hhits->Draw("same");
 221
 222                 hrecs->SetMarkerStyle(2);
 223                 hrecs->SetMarkerColor(4);
 224                 hrecs->SetMarkerSize(3);
 225                 hrecs->SetStats(kFALSE);
 226                 hrecs->Draw("same");
 227
 228                 TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
 229                 legend->SetMargin(0.2);
 230                 legend->AddEntry(hhits, "hits","P");
 231                 legend->AddEntry(hrecs, "recpoints","P");
 232                 legend->AddEntry(hdigits, "digits","P");
 233                 legend->SetTextSizePixels(14);
 234                 legend->Draw();
 235
 236                 viewer->Update();
 237         }
 238
 239         cout << "Done. Goodbye" << endl;
 240         return;
 241 }
 242
 243 Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {
 244         // First of all, the macro selects the specified module,
 245         // then gets the array of hits in it and their number.
 246         AliITS *ITS = (AliITS*) its;
 247         AliITSmodule *module = ITS->GetModule(ID);
 248         TObjArray *hits_array = module->GetHits();
 249         Int_t hits_num = hits_array->GetEntriesFast();
 250
 251         // Now, if this count returns 0, there's nothing to do,
 252         // while, if it doesn't, the first thing to do is dimensioning
 253         // the coordinate arrays, and then the loop can start.
 254         if (!hits_num)
 255                 return 0;
 256         else {
 257                 if (X) delete [] X;
 258                 if (Y) delete [] Y;
 259                 if (Z) delete [] Z;
 260                 if (St) delete [] St;
 261                 X = new Float_t[hits_num];
 262                 Y = new Float_t[hits_num];
 263                 Z = new Float_t[hits_num];
 264                 St = new Int_t[hits_num];
 265         }
 266
 267         for (Int_t j = 0; j < hits_num; j++) {
 268                 AliITShit *hit = (AliITShit*) hits_array->At(j);
 269                 X[j]  = hit->GetXL();
 270                 Y[j]  = hit->GetYL();
 271                 Z[j]  = hit->GetZL();
 272                 St[j] = hit->StatusEntering();
 273         }
 274         return hits_num;
 275 }
 276
 277 Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z) {
 278
 279         // First of all, the macro selects the specified module,
 280         // then gets the array of recpoints in it and their number.
 281         AliITS *ITS = (AliITS*) its;
 282         TTree *TR = gAlice->TreeR();
 283         ITS->ResetRecPoints();
 284         TR->GetEvent(ID);
 285         TClonesArray *recs_array = ITS->RecPoints();
 286         Int_t recs_num = recs_array->GetEntries();
 287
 288         // Now, if this count returns 0, there's nothing to do,
 289         // while, if it doesn't, the first thing to do is dimensioning
 290         // the coordinate and energy loss arrays, and then the loop can start.
 291         if (!recs_num)
 292                 return 0;
 293         else {
 294                 if (X) delete [] X;
 295                 if (Z) delete [] Z;
 296                 X = new Float_t[recs_num];
 297                 Z = new Float_t[recs_num];
 298         }
 299         for(Int_t j = 0; j < recs_num; j++) {
 300                 AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
 301            X[j] = recp->GetX();
 302                 Z[j] = recp->GetZ();
 303         }
 304         return recs_num;
 305 }
 306
 307 Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {
 308
 309         // First of all, the macro selects the specified module,
 310         // then gets the array of recpoints in it and their number.
 311         AliITS *ITS = (AliITS*)its;
 312         TTree *TD = gAlice->TreeD();
 313         ITS->ResetDigits();
 314         TD->GetEvent(ID);
 315         TClonesArray *digits_array = ITS->DigitsAddress(dtype);
 316         AliITSgeom *gm = ITS->GetITSgeom();
 317         AliITSDetType *det = ITS->DetType(dtype);
 318         AliITSsegmentation *seg = det->GetSegmentationModel();
 319         TArrayI ssdone(5000);  // used to match p and n side digits of strips
 320         TArrayI pair(5000);    // as above
 321         Int_t digits_num = digits_array->GetEntries();
 322         // Now, if this count returns 0, there's nothing to do,
 323         // while, if it doesn't, the first thing to do is dimensioning
 324         // the coordinate and energy loss arrays, and then the loop can start.
 325
 326         if(dtype==2){
 327            Int_t layer, ladder, detec;
 328            gm->GetModuleId(ID,layer,ladder,detec);
 329            seg->SetLayer(layer);
 330         }
 331
 332         if (!digits_num)
 333                 return 0;
 334         else {
 335                 if (X) delete [] X;
 336                 if (Z) delete [] Z;
 337                 X = new Float_t[digits_num];
 338                 Z = new Float_t[digits_num];
 339         }
 340
 341         // Get the coordinates of the module
 342         if (dtype == 2) {
 343                 for (Int_t j = 0; j < digits_num; j++) {
 344                         ssdone[j] = 0;
 345                         pair[j] = 0;
 346                 }
 347         }
 348   for (Int_t j = 0; j < digits_num; j++) {
 349         cout << j << endl;
 350                 digit = (AliITSdigit*)digits_array->UncheckedAt(j);
 351                 Int_t iz=digit->fCoord1;  // cell number z
 352                 Int_t ix=digit->fCoord2;  // cell number x
 353     // Get local coordinates of the element (microns)
 354                 if(dtype < 2)
 355         seg->DetToLocal(ix, iz, X[j], Z[j]);
 356     else {
 357                         // SSD: if iz==0 ---> N side; if iz==1 P side
 358       if (ssdone[j] == 0) {
 359                                 ssdone[j]=1;
 360                                 pair[j]=-1;
 361                                 Bool_t pside = (iz == 1);
 362                                 Bool_t impaired = kTRUE;
 363                                 Int_t pstrip = 0;
 364                                 Int_t nstrip = 0;
 365                                 if (pside) pstrip = ix; else nstrip = ix;
 366                                 for (Int_t k = 0; k < digits_num; k++) {
 367                                         if (ssdone[k] == 0 && impaired) {
 368                                                 AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
 369                                                 if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
 370                                                         ssdone[k]=2;
 371                                                         pair[j]=k;
 372                                                         if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
 373                                                         impaired=kFALSE;
 374                                                 }
 375                                         }
 376                                 }
 377         if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
 378         X[j] /= 10000.0;  // convert microns to cm
 379         Z[j] /= 10000.0;  // convert microns to cm
 380                         }
 381                 }
 382         }
 383         return digits_num;
 384 }
 385