From 17973466f21c8e1d0cf0a20811c678046918a76d Mon Sep 17 00:00:00 2001 From: barbera Date: Thu, 13 Jun 2002 06:41:00 +0000 Subject: [PATCH] Macro to display hits/digits/recpoints modified to take into account multiple files. This macro should be run in compiled mode. --- ITS/AliITSReadPlotData.C | 478 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 478 insertions(+) create mode 100644 ITS/AliITSReadPlotData.C diff --git a/ITS/AliITSReadPlotData.C b/ITS/AliITSReadPlotData.C new file mode 100644 index 00000000000..9c9e195ea6d --- /dev/null +++ b/ITS/AliITSReadPlotData.C @@ -0,0 +1,478 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St); +Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints); +Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z); + +Int_t AliITSReadPlotData( Int_t evNum = 0, char *filename = "galice.root", TString FileDigits="galice.root", TString FileRec="galice.root", Int_t isfastpoints = 0) { + + /********************************************************************* + * * + * Macro used to read hits, digits and recpoints for a module * + * It draws 3 TH2Fs where stores the 2-D coordinates of these * + * data for a specified module (for which the macro asks when * + * it starts, with some checks to avoid wrong detector coords. * + * * + * Only a little 'experimental' warning: * + * with all the tests I have made, I wasn't able to plot the * + * digits fo the 5th layer... * + * I skipped this problem with an alert to beware th user, while * + * in this situation the macro plots only recpoints and hits * + * * + * Author: Alberto Pulvirenti * + * * + *********************************************************************/ + + + // First of all, here are put some variable declarations + // that are useful in the following part: + Int_t nparticles; // number of particles + // ITS module coordinates [layer = 1, ladder = 2, det = 3] and absolute ID[0] of module [0] + TArrayI ID(4); + Int_t nmodules, dtype; // Total number of modules and module type (SSD, SPD, SDD) + Float_t *x = 0, *y = 0, *z = 0; // Arrays where to store read coords + Bool_t *St = 0; // Status of the track (hits only) + /* + // It's necessary to load the gAlice shared libs + // if they aren't already stored in memory... + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + // Anyway, this macro needs to read a gAlice file, so it + // clears the gAlice object if there is already one in memory... + else { + */ + if(gAlice){ + delete gAlice; + gAlice = 0; + } + // } + + // Now is opened the Root input file containing Geometry, Kine and Hits + // by default its name must be "galice.root". + // When the file is opened, its contens are shown. + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename); + if (!file) file = new TFile(filename); + file->ls(); + + // Then, the macro gets the AliRun object from file. + // If this object is not present, an error occurs + // and the execution is stopped. + // Anyway, this operation needs some time, + // don't worry about an apparent absence of output and actions... + cout << "\nSearching in '" << filename << "' for an AliRun object ... " << flush; + gAlice = (AliRun*)file->Get("gAlice"); + if (gAlice) + cout << "FOUND!" << endl; + else { + cout<<"NOT FOUND! The Macro can't continue!" << endl; + return 0; + } + if(!(FileDigits.Data() == filename)){ + gAlice->SetTreeDFileName(FileDigits); + } + if(!(FileRec.Data() == filename)){ + gAlice->SetTreeRFileName(FileRec); + } + + // Then, the macro selects the event number specified. Default is 0. + nparticles = gAlice->GetEvent(evNum); + cout << "\nNumber of particles = " << nparticles << endl; + if (!nparticles) { + cout << "With no particles I can't do much... Goodbye!" << endl; + return 0; + } + + // The next step is filling the ITS from the AliRun object. + AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); + ITS->InitModules(-1, nmodules); + cout << "Number of ITS modules = " << nmodules << endl; + cout << "\nFilling modules (it takes a while, now)..." << flush; + ITS->FillModules(0, 0, nmodules, " ", " "); + cout << "DONE!" << endl; + AliITSgeom *gm = ITS->GetITSgeom(); + // AliITSDetType *det = ITS->DetType(dtype); + // AliITSsegmentation *seg = det->GetSegmentationModel(); + + for(;;) { + + // Input phase. + // The macro asks if the user wants to put a detector ID[0] + // or prefers to input layer, ladder and detector. + for (Int_t i = 0; i < 4; i++) ID[i] = -1; + Int_t answ; + do { + cout << "\nSelection modes:" << endl; + cout << "1) by layer - ladder - detector numbers" << endl; + cout << "2) by unique detector ID" << endl; + cout << "0) exit macro" << endl; + cout << "\nEnter your choice: "; + cin >> answ; + } while (answ < 0 || answ > 2); + switch (answ) { + case 0: + // input = 0 ---> EXIT + return 0; + break; + case 1: + // input = 1 ---> SELECTION BY COORD + do { + cout << "\nLayer number [1-6, 0 to exit]: "; + cin >> ID[1]; + if (!ID[1]) return 0; + } while (ID[1] < 0 || ID[1] > 6); + + // Detector type: 0 --> SPD, 1 --> SDD, 2 --> SSD. + // Layer 1,2 --> 0 / Layer 3,4 --> 1 / Layer 5,6 --> 2 + dtype = (ID[1] - 1) / 2; + + // Once fixed the layer number, the macro calculates the max number + // for ladder and detector from geometry, and accepts only suitable values. + do { + ID[2] = gm->GetNladders(ID[1]); + cout << "Ladder number [1-" << ID[2] << ", 0 to exit]: "; + cin >> ID[2]; + if (!ID[2]) return 0; + } while (ID[2] < 0 || ID[2] > gm->GetNladders(ID[1])); + do { + ID[3] = gm->GetNdetectors(ID[1]); + cout << "Detector number [1-" << ID[3] << ", 0 to exit]: "; + cin >> ID[3]; + if (!ID[3]) return 0; + } while (ID[3] < 0 || ID[3] > gm->GetNdetectors(ID[1])); + break; + case 2: + // input = 2 ---> SELECTION BY ID[0] + do { + ID[0] = gm->GetIndexMax(); + cout << "\n Detector ID number [0-" << ID[0]-1 << ", -1 to exit]: "; + cin >> ID[0]; + ID[0]++; + if (ID[0] == -1) return 0; + } while (ID[0] < 0 || ID[0] > gm->GetIndexMax()); + break; + }; + + if (ID[0] == -1) + // If ID[0] = -1 the chioce was by coords, so it's necessary to assign the ID: + ID[0] = gm->GetModuleIndex(ID[1], ID[2], ID[3]); + else { + // Else we must get the layer, ladder and detector by the ID; + // ...but first we must remember that the ID goest from 0 to NModules - 1! + ID[0]--; + ID[1] = ITS->GetModule(ID[0])->GetLayer(); + ID[2] = ITS->GetModule(ID[0])->GetLadder(); + ID[3] = ITS->GetModule(ID[0])->GetDet(); + } + + // Defines the histograms inside the `for' cycle, so they are destroyed at the end + // of every read sequqnce, in order to mek another withour segmentation faults + Text_t msg[250]; + Float_t xm = 0.0, zm = 0.0; + dtype = (ID[1] - 1) / 2; + switch (dtype) { + case 0: xm = 1.5; zm = 7.0; break; + case 1: xm = 7.5; zm = 8.0; break; + case 2: xm = 7.5; zm = 4.5; break; + } + sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]); + TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of hits + TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of recpoints + TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits + + cout << endl; + + // Reads hits... + Int_t tmpint = ID[0]; + Int_t hits = GetModuleHits(ITS, tmpint, x, y, z, St); + if (!hits) { + cout << "No hits in module!" << endl; + continue; + } + else + cout << "Hits scanned..." << endl; + for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]); + + // Reads recpoints... + Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z, isfastpoints); + if (!recs) { + cout << "No recpoints in module!" << endl; + continue; + } + else + cout << "Recpoints scanned..." << endl; + for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]); + + // Reads digits... + Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z); + if (!digits) { + cout << "No digits in module!" << endl; + //continue; + } + else + cout << "Digits scanned..." << endl; + for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]); + + // Draws read data... + // HITS -------> red (2) crosses. + // DIGITS -----> green (8) boxes. + // REC-POINTS -> blue (4) St. Andrew's crosses. + + TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800); + viewer->cd(); + + hdigits->SetMarkerStyle(25); + hdigits->SetMarkerColor(8); + hdigits->SetMarkerSize(2); + hdigits->SetStats(kFALSE); + hdigits->SetXTitle("Local X (cm)"); + hdigits->SetYTitle("Local Z (cm)"); + hdigits->Draw(); + + hhits->SetMarkerStyle(5); + hhits->SetMarkerColor(2); + hhits->SetMarkerSize(3); + hhits->SetStats(kFALSE); + hhits->Draw("same"); + + hrecs->SetMarkerStyle(2); + hrecs->SetMarkerColor(4); + hrecs->SetMarkerSize(3); + hrecs->SetStats(kFALSE); + hrecs->Draw("same"); + + TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9); + legend->SetMargin(0.2); + legend->AddEntry(hhits, "hits","P"); + legend->AddEntry(hrecs, "recpoints","P"); + legend->AddEntry(hdigits, "digits","P"); + legend->SetTextSizePixels(14); + legend->Draw(); + + viewer->Update(); + + + Text_t fname[250],ans; + cout << "Do you want to save the current canvas on a file (y/n) ? "; + cin >> ans; + if(ans == 'y' || ans == 'Y') { + cout << "Enter filename: "; + cin >> fname; + TString *control = new TString(fname); + Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif"); + if(!ok){ + cout << "File extension is not recognized. The canvas will be saved as Postscript file"; + strcat(fname,".ps"); + } + viewer->SaveAs(fname); + } + } + + cout << "Done. Goodbye" << endl; + return 0; +} + +Int_t GetModuleHits(TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) { + // First of all, the macro selects the specified module, + // then gets the array of hits in it and their number. + AliITS *ITS = (AliITS*) its; + AliITSmodule *module = ITS->GetModule(ID); + TObjArray *hits_array = module->GetHits(); + Int_t hits_num = hits_array->GetEntriesFast(); + + // Now, if this count returns 0, there's nothing to do, + // while, if it doesn't, the first thing to do is dimensioning + // the coordinate arrays, and then the loop can start. + if (!hits_num) + return 0; + else { + if (X) delete [] X; + if (Y) delete [] Y; + if (Z) delete [] Z; + if (St) delete [] St; + X = new Float_t[hits_num]; + Y = new Float_t[hits_num]; + Z = new Float_t[hits_num]; + St = new Bool_t[hits_num]; + } + + for (Int_t j = 0; j < hits_num; j++) { + AliITShit *hit = (AliITShit*) hits_array->At(j); + X[j] = hit->GetXL(); + Y[j] = hit->GetYL(); + Z[j] = hit->GetZL(); + St[j] = hit->StatusEntering(); + } + return hits_num; +} + +Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints) { + + // First of all, the macro selects the specified module, + // then gets the array of recpoints in it and their number. + AliITS *ITS = (AliITS*) its; + TTree *TR = gAlice->TreeR(); + ITS->ResetRecPoints(); + + TClonesArray *recs_array = ITS->RecPoints(); + + TBranch *branch = 0; + if(TR && recs_array){ + if(isfastpoints==1){ + branch = gAlice->TreeR()->GetBranch("ITSRecPointsF"); + cout<<"using fast points\n"; + } + else { + branch = gAlice->TreeR()->GetBranch("ITSRecPoints"); + } + if(branch)branch->SetAddress(&recs_array); + } + + branch->GetEvent(ID); + + Int_t recs_num = recs_array->GetEntries(); + + + // Now, if this count returns 0, there's nothing to do, + // while, if it doesn't, the first thing to do is dimensioning + // the coordinate and energy loss arrays, and then the loop can start. + if (!recs_num) + return 0; + else { + if (X) delete [] X; + if (Z) delete [] Z; + X = new Float_t[recs_num]; + Z = new Float_t[recs_num]; + } + for(Int_t j = 0; j < recs_num; j++) { + AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j); + X[j] = recp->GetX(); + Z[j] = recp->GetZ(); + } + return recs_num; +} + +Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) { + + // First of all, the macro selects the specified module, + // then gets the array of recpoints in it and their number. + AliITSdigit *digit; + AliITS *ITS = (AliITS*)its; + TTree *TD = gAlice->TreeD(); + ITS->ResetDigits(); + TD->GetEvent(ID); + TClonesArray *digits_array = ITS->DigitsAddress(dtype); + AliITSgeom *gm = ITS->GetITSgeom(); + AliITSDetType *det = ITS->DetType(dtype); + AliITSsegmentation *seg = det->GetSegmentationModel(); + TArrayI ssdone(5000); // used to match p and n side digits of strips + TArrayI pair(5000); // as above + Int_t digits_num = digits_array->GetEntries(); + // Now, if this count returns 0, there's nothing to do, + // while, if it doesn't, the first thing to do is dimensioning + // the coordinate and energy loss arrays, and then the loop can start. + + if(dtype==2){ + Int_t layer, ladder, detec; + gm->GetModuleId(ID,layer,ladder,detec); + seg->SetLayer(layer); + } + + if (!digits_num) + return 0; + else { + cout << "Digits to scan: " << digits_num << endl; + if (X) delete [] X; + if (Z) delete [] Z; + X = new Float_t[digits_num]; + Z = new Float_t[digits_num]; + } + + // Get the coordinates of the module + if (dtype == 2) { + for (Int_t j = 0; j < digits_num; j++) { + ssdone[j] = 0; + pair[j] = 0; + } + } + // AliITSdigit *digit; + for (Int_t j = 0; j < digits_num; j++) { + digit = (AliITSdigit*)digits_array->UncheckedAt(j); + Int_t iz=digit->fCoord1; // cell number z + Int_t ix=digit->fCoord2; // cell number x + // Get local coordinates of the element (microns) + // ******************************* PARTE CORRETTA *************************************** + if(dtype < 2) { + Float_t xx, zz; // aggiunta + seg->DetToLocal(ix, iz, xx, zz); + X[j] = xx; // aggiunta + Z[j] = zz; // aggiunta + } + // ******************************* FINE PARTE CORRETTA *************************************** + else { + // SSD: if iz==0 ---> N side; if iz==1 P side + if (ssdone[j] == 0) { + ssdone[j]=1; + pair[j]=-1; + Bool_t pside = (iz == 1); + Bool_t impaired = kTRUE; + Int_t pstrip = 0; + Int_t nstrip = 0; + if (pside) pstrip = ix; else nstrip = ix; + for (Int_t k = 0; k < digits_num; k++) { + if (ssdone[k] == 0 && impaired) { + AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k); + if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) { + ssdone[k]=2; + pair[j]=k; + if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2; + impaired=kFALSE; + } + } + } + if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]); + X[j] /= 10000.0; // convert microns to cm + Z[j] /= 10000.0; // convert microns to cm + } + } + } + return digits_num; +} + + + + + + + + + + -- 2.31.1