ITS/AliITSReadPlotData.C

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