1 #if !defined(__CINT__) || defined(__MAKECINT__)
7 #include<TClassTable.h>
8 #include<TClonesArray.h>
15 #include "AliGenEventHeader.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 #include <AliRunLoader.h>
29 #include <AliITSLoader.h>
30 #include <AliHeader.h>
32 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
33 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
34 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
36 Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+Rec", char *filename="galice.root", Int_t isfastpoints = 0) {
37 /*******************************************************************
38 * This macro displays geometrical information related to the
39 * hits, digits and rec points in ITS.
40 * There are histograms that are not displayed (i.e. energy
41 * deposition) but that are saved onto a file (see below)
45 * Int_t evesel: the selected event number. It's =0 by default
47 * Options: Any combination of:
48 * 1) subdetector name: "SPD", "SDD", "SSD", "All" (default)
49 * 2) Printouts: "Verbose" Almost everything is printed out:
50 * it is wise to redirect the output onto a file
51 * e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
52 * 3) Rec Points option: "Rec" ---> Uses Rec. Points (default)
53 * otherwise ---> uses hits and digits only
55 * .x AliITSGeoPlot(); (All subdetectors; no-verbose; no-recpoints)
56 * .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
58 * filename: It's "galice.root" by default.
59 * isfastpoints: integer. It is set to 0 by defaults. This means that
60 * slow reconstruction is assumed. If fast recpoint are in
61 * in use, isfastpoints must be set =1.
62 * BEWARE: to analyze digits/recpoints produced with aliroot
63 * versions older than v3-07-Release, the input parameters
64 * filename, FileDigits, FileRec and isfastpoints must be left
65 * to their default value
67 * OUTPUT: It produces a root file with a list of histograms
69 * WARNING: spatial information for SSD/DIGITS is obtained by pairing
70 * digits on p and n side originating from the same track, when
71 * possible. This (mis)use of DIGITS is tolerated for debugging
72 * purposes only !!!! The pairing in real life should be done
73 * starting from info really available...
75 * COMPILATION: this macro can be compiled.
76 * 1) You need to set your include path with
77 * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
78 * 3) If you are using root instead of aliroot you need to
79 * execute the macro loadlibs.C in advance
80 * 4) To compile this macro from root (or aliroot):
81 * --- .L AliITSGeoPlot.C++
82 * --- AliITSGeoPlot();
84 * M.Masera 14/05/2001 18:30
85 * Last rev. 09/06/2003 17:00 (Adapted to NewIO) m.m.
86 ********************************************************************/
90 Bool_t All = choice.Contains("All");
91 Bool_t verbose=choice.Contains("Verbose");
92 Bool_t userec=choice.Contains("Rec");
93 Int_t retcode=1; //return code
94 #if !(!defined(__CINT__) || defined(__MAKECINT__))
95 if (gClassTable->GetID("AliRun") < 0) {
96 gROOT->LoadMacro("loadlibs.C");
102 delete gAlice->GetRunLoader();
106 #if !(!defined(__CINT__) || defined(__MAKECINT__))
110 AliRunLoader* rl = AliRunLoader::Open(filename);
112 cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
115 Int_t retval = rl->LoadgAlice();
117 cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
120 gAlice=rl->GetAliRun();
122 retval = rl->LoadHeader();
124 cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
128 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
131 cerr<<"AliITSGeoPlot.C : ITS loader not found"<<endl;
135 ITSloader->LoadHits("read");
136 ITSloader->LoadDigits("read");
137 if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
138 ITSloader->LoadRecPoints("read");
139 rl->GetEvent(evesel);
140 Int_t nparticles = rl->GetHeader()->GetNtrack();
141 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
142 ITS->SetTreeAddress();
144 cout<<" "<<endl<<" "<<endl;
145 cout<<"******* Event processing started *******"<<endl;
146 cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
147 cout<<"will not be processed"<<endl;
148 cout << "Number of particles= " << nparticles <<endl;
152 TTree *TH = ITSloader->TreeH();
153 Stat_t ntracks = TH->GetEntries();
154 if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
158 ITS->InitModules(-1,nmodules);
159 cout<<"Number of ITS modules= "<<nmodules<<endl;
160 cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
161 ITS->FillModules(0,0,nmodules," "," ");
162 cout<<"ITS modules .... DONE!"<<endl;
165 TTree *TD = ITSloader->TreeD();
168 TTree *TR = ITSloader->TreeR();
169 TClonesArray *ITSrec = ITS->RecPoints();
171 if(userec && TR && ITSrec){
173 branch = ITSloader->TreeR()->GetBranch("ITSRecPointsF");
174 cout<<"using fast points\n";
177 branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
179 if(branch)branch->SetAddress(&ITSrec);
182 if(userec && (!TR || !ITSrec || !branch)){
184 cout<<"\n ======================================================= \n";
185 cout<<"WARNING: there are no RECPOINTS on this file ! \n";
186 cout<<"======================================================= \n \n";
190 Int_t mod; //module number
192 Double_t pos[3]; // Global position of the current module
193 Float_t ragdet; // Radius of detector (x y plane)
194 Int_t first; // first module
195 Int_t last; // last module
196 Int_t nrecp; //number of RecPoints for a given module
199 TObjArray histos(26,0); // contains the pointers to the histograms
200 // Book histograms SPD
201 TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
202 TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
203 TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
204 TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
205 TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
206 TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
207 TH2F *rrspd = new TH2F("rrspd","Radii of recpoints - SPD",50,-10.,10.,50,-10.,10.);
208 TH1F *zrspd = new TH1F("zrspd","Z of recpoints - SPD",100,-30.,30.);
209 TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
210 histos.AddLast(rspd); // 0
211 histos.AddLast(rhspd); // 1
212 histos.AddLast(rmspd); // 2
213 histos.AddLast(zspd); // 3
214 histos.AddLast(zhspd); // 4
215 histos.AddLast(zmspd); // 5
216 histos.AddLast(rrspd); // 6
217 histos.AddLast(zrspd); // 7
218 histos.AddLast(enespd); // 8
219 // Book histograms SDD
220 TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
221 TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
222 TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
223 TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
224 TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
225 TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
226 TH2F *rrsdd = new TH2F("rrsdd","Radii of recpoints - SDD",50,-40.,40.,50,-40.,40.);
227 TH1F *zrsdd = new TH1F("zrsdd","Z of recpoints - SDD",100,-40.,40.);
228 TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
229 histos.AddLast(rsdd); // 9
230 histos.AddLast(rhsdd); // 10
231 histos.AddLast(rmsdd); // 11
232 histos.AddLast(zsdd); // 12
233 histos.AddLast(zhsdd); // 13
234 histos.AddLast(zmsdd); // 14
235 histos.AddLast(rrsdd); // 15
236 histos.AddLast(zrsdd); // 16
237 histos.AddLast(enesdd); // 17
238 // Book histogram SSD
239 TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
240 TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
241 TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
242 TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
243 TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
244 TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
245 TH2F *rrssd = new TH2F("rrssd","Radii of recpoints - SSD",50,-50.,50.,50,-50.,50.);
246 TH1F *zrssd = new TH1F("zrssd","Z of recpoints - SSD",100,-70.,70.);
247 TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
248 histos.AddLast(rssd); // 18
249 histos.AddLast(rhssd); // 19
250 histos.AddLast(rmssd); // 20
251 histos.AddLast(zssd); // 21
252 histos.AddLast(zhssd); // 22
253 histos.AddLast(zmssd); // 23
254 histos.AddLast(rrssd); // 24
255 histos.AddLast(zrssd); // 25
256 histos.AddLast(enessd); // 26
258 // Loop on subdetectors
260 AliITSgeom *geom = ITS->GetITSgeom();
261 TString detna; // subdetector name
262 for(Int_t subd=0;subd<3;subd++){
263 if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
264 // Prepare array for the digits
265 TClonesArray *ITSdigits = ITS->DigitsAddress(subd);
266 Bool_t usedigits = kTRUE;
269 cout<<"\n ======================================================= \n";
270 cout<<"WARNING: there are no DIGITS on this file ! \n";
271 cout<<"======================================================= \n \n";
273 // Get segmentation model
274 AliITSDetType *iDetType=ITS->DetType(subd);
275 if(subd==0)detna="SPD";
276 if(subd==1)detna="SDD";
277 if(subd==2)detna="SSD";
278 AliITSsegmentation *seg=(AliITSsegmentation*)iDetType->GetSegmentationModel();
280 first = geom->GetStartDet(subd);
281 last = geom->GetLastDet(subd);
283 cout<<" "<<endl<<"-------------------------------------"<<endl;
284 cout<<"Start processing subdetector "<<detna<<endl;
285 cout<<detna<<" first module "<<first<<endl;
286 cout<<detna<<" last module "<<last<<endl;
287 cout<<" "<<endl<<" "<<endl;
289 for (mod=first; mod<=last; mod++){
290 geom->GetTrans(mod,pos); // position of the module in the MRS
291 ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
292 // The following 2 histos are a check of the geometry
293 TH2F *bidi = (TH2F*)histos.At(2+subd*9);
294 TH1F *uni = (TH1F*)histos.At(5+subd*9);
295 bidi->Fill(pos[0],pos[1]);
298 cout<<"=========================================================\n";
299 cout<<detna<<" module="<<mod<<endl;
300 cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", ";
301 cout<<pos[2]<<" Radius "<<ragdet<<endl;
305 GetHitsCoor(ITS,mod,histos,subd,verbose);
309 ITS->ResetRecPoints();
310 // TR->GetEvent(mod);
311 branch->GetEvent(mod);
312 TH2F *bidi=(TH2F*)histos.At(6+subd*9);
313 TH1F *uni=(TH1F*)histos.At(7+subd*9);
314 nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
320 nbytes += TD->GetEvent(mod);
321 GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
324 } // End of loop on the modules
327 // Plot the histograms
328 TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
331 TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
337 TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
343 TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
348 h2tmp = (TH2F*)histos.At(9*subd);
351 h1tmp=(TH1F*)histos.At(3+subd*9);
354 h2tmp=(TH2F*)histos.At(1+9*subd);
357 h1tmp=(TH1F*)histos.At(4+subd*9);
362 h2tmp=(TH2F*)histos.At(6+9*subd);
365 h1tmp=(TH1F*)histos.At(7+subd*9);
371 h2tmp=(TH2F*)histos.At(2+9*subd);
374 h2tmp=(TH2F*)histos.At(5+9*subd);
378 } // end of loop on subdetectors
379 // Save the histograms
380 TFile *fh = new TFile("AliITSGeoPlot.root","recreate");
381 // The list is written to file as a single entry
382 TList *lihi = new TList();
383 // copy the pointers to the histograms to a TList object.
384 // The histograms concerning recpoints are not copied if
385 // 'userec' is false.
386 for(Int_t i=0;i<histos.GetEntriesFast();i++){
387 if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
388 if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
391 lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
398 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
399 TH2F *h2=(TH2F*)histos.At(1+subd*9);
400 TH1F *h1=(TH1F*)histos.At(4+subd*9);
401 TH1F *ener=(TH1F*)histos.At(8+subd*9);
402 AliITS *ITS= (AliITS*)its;
403 AliITSmodule *modu = ITS->GetModule(mod);
404 TObjArray *fHits = modu->GetHits();
405 Int_t nhits = fHits->GetEntriesFast();
408 cout<<"-------------------------------------------------------"<<endl;
409 cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
411 for (Int_t hit=0;hit<nhits;hit++) {
412 AliITShit *iHit = (AliITShit*) fHits->At(hit);
413 if(!iHit->StatusEntering()){
414 Float_t x=iHit->GetXG();
415 Float_t y=iHit->GetYG();
416 Float_t z=iHit->GetZG();
417 Float_t edep=iHit->GetIonization()*1000000;
422 cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
423 cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= ";
432 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
433 AliITSgeom *geom = (AliITSgeom*)ge;
434 Int_t nrecp = ITSrec->GetEntries();
436 Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
437 Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
439 cout<<"-------------------------------------------------------"<<endl;
440 cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
442 for(Int_t irec=0;irec<nrecp;irec++) {
443 AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
446 geom->LtoG(mod,lc,gc);
448 cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
450 Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
451 cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
453 cout<<"Associated track "<<recp->GetLabel(0)<<endl;
455 h2->Fill(gc[0],gc[1]);
462 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
463 AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
464 AliITSgeom *geom = (AliITSgeom*)ge;
469 geom->GetModuleId(mod,layer,ladder,detec);
470 seg->SetLayer(layer);
472 Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.; //local coord dig.
473 Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
474 Float_t ragdig; // Radius digit
475 TArrayI ssdone(5000); // used to match p and n side digits of strips
476 TArrayI pair(5000); // as above
477 Int_t ndigits = ITSdigits->GetEntries();
481 cout<<"-------------------------------------------------------"<<endl;
482 cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
484 // Get the coordinates of the module
486 for (Int_t digit=0;digit<ndigits;digit++){
491 for (Int_t digit=0;digit<ndigits;digit++) {
492 digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
493 Int_t iz=digs->GetCoord1(); // cell number z
494 Int_t ix=digs->GetCoord2(); // cell number x
495 // Get local coordinates of the element
497 seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
500 // SSD: if iz==0 ---> N side; if iz==1 P side
501 if(ssdone[digit]==0){
504 Bool_t pside=(iz==1);
505 Bool_t impaired=kTRUE;
510 for(Int_t digi2=0;digi2<ndigits;digi2++){
511 if(ssdone[digi2]==0 && impaired){
512 AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
513 if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
516 if(pside)nstrip=dig2->GetCoord2();
517 if(!pside)pstrip=dig2->GetCoord2();
522 if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
525 if(subd<2 || (subd==2 && ssdone[digit]==1)){
526 Int_t coor1=digs->GetCoord1();
527 Int_t coor2=digs->GetCoord2();
528 Int_t tra0=digs->GetTrack(0);
530 cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
531 cout<<coor2<<" track "<<tra0<<" "<<endl;
532 if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z=";
533 cout<<lcoor[2]<<endl;
536 cout<<"(digit not paired)"<<endl;
539 cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2];
541 Int_t dtmp=pair[digit];
542 AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
543 Int_t coor1b=dig2->GetCoord1();
544 Int_t coor2b=dig2->GetCoord2();
545 Int_t tra0b=dig2->GetTrack(0);
546 cout<<"(digit paired with digit #"<<dtmp<<endl;
547 cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
548 cout<<" track "<<tra0b<<")"<<endl;
552 if(subd<2 || (subd==2 && pair[digit]!=-1)){
553 // Global coordinates of the element
554 //SDD and SPD use cm, SSD microns (GetPadCxz)
555 if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
557 geom->LtoG(mod,lcoor,gcoor); // global coord. in cm
558 ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
560 cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
561 cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
564 TH2F *bidi = (TH2F*)histos.At(subd*9);
565 TH1F *uni = (TH1F*)histos.At(3+subd*9);
566 bidi->Fill(gcoor[0],gcoor[1]);
570 } // loop on digits for this module
571 } // if(ndigits>0....