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                 else
 183                         cout << "Hits scanned..." << endl;
 184                 for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
 185
 186                 // Reads recpoints...
 187                 Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z);
 188                 if (!recs) {
 189                         cout << "No recpoints in module!" << endl;
 190                         continue;
 191                 }
 192                 else
 193                         cout << "Recpoints scanned..." << endl;
 194                 for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
 195
 196                 // Reads digits...
 197                 Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
 198                 if (!digits) {
 199                         cout << "No digits in module!" << endl;
 200                         //continue;
 201                 }
 202                 else
 203                         cout << "Digits scanned..." << endl;
 204                 for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
 205
 206                 // Draws read data...
 207                 // HITS -------> red (2) crosses.
 208                 // DIGITS -----> green (8) boxes.
 209                 // REC-POINTS -> blue (4) St. Andrew's crosses.
 210
 211                 TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
 212                 viewer->cd();
 213
 214                 hdigits->SetMarkerStyle(25);
 215                 hdigits->SetMarkerColor(8);
 216                 hdigits->SetMarkerSize(2);
 217                 hdigits->SetStats(kFALSE);
 218                 hdigits->SetXTitle("Local X (cm)");
 219                 hdigits->SetYTitle("Local Z (cm)");
 220                 hdigits->Draw();
 221
 222                 hhits->SetMarkerStyle(5);
 223                 hhits->SetMarkerColor(2);
 224                 hhits->SetMarkerSize(3);
 225                 hhits->SetStats(kFALSE);
 226                 hhits->Draw("same");
 227
 228                 hrecs->SetMarkerStyle(2);
 229                 hrecs->SetMarkerColor(4);
 230                 hrecs->SetMarkerSize(3);
 231                 hrecs->SetStats(kFALSE);
 232                 hrecs->Draw("same");
 233
 234                 TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
 235                 legend->SetMargin(0.2);
 236                 legend->AddEntry(hhits, "hits","P");
 237                 legend->AddEntry(hrecs, "recpoints","P");
 238                 legend->AddEntry(hdigits, "digits","P");
 239                 legend->SetTextSizePixels(14);
 240                 legend->Draw();
 241
 242                 viewer->Update();
 243
 244
 245                 Text_t fname[250],ans;
 246                 cout << "Do you want to save the current canvas on a file (y/n) ? ";
 247                 cin >> ans;
 248                 if(ans == 'y' || ans == 'Y') {
 249                    cout << "Enter filename: ";
 250                    cin >> fname;
 251                    TString *control = new TString(fname);
 252                    Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif");
 253                    if(!ok){
 254                       cout << "File extension is not recognized. The canvas will be saved as Postscript file";
 255                       strcat(fname,".ps");
 256                    }
 257                    viewer->SaveAs(fname);
 258                 }
 259         }
 260
 261         cout << "Done. Goodbye" << endl;
 262         return;
 263 }
 264
 265 Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {
 266         // First of all, the macro selects the specified module,
 267         // then gets the array of hits in it and their number.
 268         AliITS *ITS = (AliITS*) its;
 269         AliITSmodule *module = ITS->GetModule(ID);
 270         TObjArray *hits_array = module->GetHits();
 271         Int_t hits_num = hits_array->GetEntriesFast();
 272
 273         // Now, if this count returns 0, there's nothing to do,
 274         // while, if it doesn't, the first thing to do is dimensioning
 275         // the coordinate arrays, and then the loop can start.
 276         if (!hits_num)
 277                 return 0;
 278         else {
 279                 if (X) delete [] X;
 280                 if (Y) delete [] Y;
 281                 if (Z) delete [] Z;
 282                 if (St) delete [] St;
 283                 X = new Float_t[hits_num];
 284                 Y = new Float_t[hits_num];
 285                 Z = new Float_t[hits_num];
 286                 St = new Int_t[hits_num];
 287         }
 288
 289         for (Int_t j = 0; j < hits_num; j++) {
 290                 AliITShit *hit = (AliITShit*) hits_array->At(j);
 291                 X[j]  = hit->GetXL();
 292                 Y[j]  = hit->GetYL();
 293                 Z[j]  = hit->GetZL();
 294                 St[j] = hit->StatusEntering();
 295         }
 296         return hits_num;
 297 }
 298
 299 Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z) {
 300
 301         // First of all, the macro selects the specified module,
 302         // then gets the array of recpoints in it and their number.
 303         AliITS *ITS = (AliITS*) its;
 304         TTree *TR = gAlice->TreeR();
 305         ITS->ResetRecPoints();
 306         TR->GetEvent(ID);
 307         TClonesArray *recs_array = ITS->RecPoints();
 308         Int_t recs_num = recs_array->GetEntries();
 309
 310         // Now, if this count returns 0, there's nothing to do,
 311         // while, if it doesn't, the first thing to do is dimensioning
 312         // the coordinate and energy loss arrays, and then the loop can start.
 313         if (!recs_num)
 314                 return 0;
 315         else {
 316                 if (X) delete [] X;
 317                 if (Z) delete [] Z;
 318                 X = new Float_t[recs_num];
 319                 Z = new Float_t[recs_num];
 320         }
 321         for(Int_t j = 0; j < recs_num; j++) {
 322                 AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
 323            X[j] = recp->GetX();
 324                 Z[j] = recp->GetZ();
 325         }
 326         return recs_num;
 327 }
 328
 329 Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {
 330
 331         // First of all, the macro selects the specified module,
 332         // then gets the array of recpoints in it and their number.
 333         AliITS *ITS = (AliITS*)its;
 334         TTree *TD = gAlice->TreeD();
 335         ITS->ResetDigits();
 336         TD->GetEvent(ID);
 337         TClonesArray *digits_array = ITS->DigitsAddress(dtype);
 338         AliITSgeom *gm = ITS->GetITSgeom();
 339         AliITSDetType *det = ITS->DetType(dtype);
 340         AliITSsegmentation *seg = det->GetSegmentationModel();
 341         TArrayI ssdone(5000);  // used to match p and n side digits of strips
 342         TArrayI pair(5000);    // as above
 343         Int_t digits_num = digits_array->GetEntries();
 344         // Now, if this count returns 0, there's nothing to do,
 345         // while, if it doesn't, the first thing to do is dimensioning
 346         // the coordinate and energy loss arrays, and then the loop can start.
 347
 348         if(dtype==2){
 349            Int_t layer, ladder, detec;
 350            gm->GetModuleId(ID,layer,ladder,detec);
 351            seg->SetLayer(layer);
 352         }
 353
 354         if (!digits_num)
 355                 return 0;
 356         else {
 357                 cout << "Digits to scan: " << digits_num << endl;
 358                 if (X) delete [] X;
 359                 if (Z) delete [] Z;
 360                 X = new Float_t[digits_num];
 361                 Z = new Float_t[digits_num];
 362         }
 363
 364         // Get the coordinates of the module
 365         if (dtype == 2) {
 366                 for (Int_t j = 0; j < digits_num; j++) {
 367                         ssdone[j] = 0;
 368                         pair[j] = 0;
 369                 }
 370         }
 371   for (Int_t j = 0; j < digits_num; j++) {
 372                 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                 // ******************************* PARTE CORRETTA ***************************************
 377                 if(dtype < 2) {
 378                         Float_t xx, zz; // aggiunta
 379         seg->DetToLocal(ix, iz, xx, zz);
 380                         X[j] = xx; // aggiunta
 381                         Z[j] = zz; // aggiunta
 382                 }
 383                 // ******************************* FINE PARTE CORRETTA ***************************************
 384     else {
 385                         // SSD: if iz==0 ---> N side; if iz==1 P side
 386       if (ssdone[j] == 0) {
 387                                 ssdone[j]=1;
 388                                 pair[j]=-1;
 389                                 Bool_t pside = (iz == 1);
 390                                 Bool_t impaired = kTRUE;
 391                                 Int_t pstrip = 0;
 392                                 Int_t nstrip = 0;
 393                                 if (pside) pstrip = ix; else nstrip = ix;
 394                                 for (Int_t k = 0; k < digits_num; k++) {
 395                                         if (ssdone[k] == 0 && impaired) {
 396                                                 AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
 397                                                 if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
 398                                                         ssdone[k]=2;
 399                                                         pair[j]=k;
 400                                                         if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
 401                                                         impaired=kFALSE;
 402                                                 }
 403                                         }
 404                                 }
 405         if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
 406                                 X[j] /= 10000.0;  // convert microns to cm
 407                                 Z[j] /= 10000.0;  // convert microns to cm
 408                         }
 409                 }
 410         }
 411         return digits_num;
 412 }
 413