Coding rule violations corrected.
[u/mrichter/AliRoot.git] / ITS / AliITSGeoPlot.C
CommitLineData
813e97e8 1#if !defined(__CINT__) || defined(__MAKECINT__)
88cb7938 2#include<Riostream.h>
af8e1c2d 3#include<TROOT.h>
4#include<TArrayI.h>
813e97e8 5#include<TBranch.h>
af8e1c2d 6#include<TCanvas.h>
7#include<TClassTable.h>
8#include<TClonesArray.h>
9#include<TFile.h>
10#include<TH1.h>
11#include<TH2.h>
a6c9fb45 12#include <TInterpreter.h>
af8e1c2d 13#include<TObject.h>
14#include<TObjArray.h>
15#include<TTree.h>
979e3647 16#include "AliGenEventHeader.h"
af8e1c2d 17#include <AliRun.h>
18#include <AliITS.h>
19#include <AliITSgeom.h>
20#include <AliITSDetType.h>
21#include <AliITSRecPoint.h>
a6c9fb45 22#include <AliITSclusterV2.h>
af8e1c2d 23#include <AliITSdigit.h>
a6c9fb45 24#include <AliITSdigitSSD.h>
af8e1c2d 25#include <AliITShit.h>
26#include <AliITSmodule.h>
27#include <AliITSsegmentation.h>
28#include <AliITSsegmentationSPD.h>
29#include <AliITSsegmentationSDD.h>
30#include <AliITSsegmentationSSD.h>
88cb7938 31#include <AliRunLoader.h>
32#include <AliITSLoader.h>
33#include <AliHeader.h>
af8e1c2d 34#endif
813e97e8 35void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
36Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
a6c9fb45 37Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
813e97e8 38void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
39
88cb7938 40Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+Rec", char *filename="galice.root", Int_t isfastpoints = 0) {
af8e1c2d 41 /*******************************************************************
42 * This macro displays geometrical information related to the
a6c9fb45 43 * hits, digits and rec points (or V2 clusters) in ITS.
af8e1c2d 44 * There are histograms that are not displayed (i.e. energy
45 * deposition) but that are saved onto a file (see below)
46 *
813e97e8 47 * INPUT arguments:
48 *
49 * Int_t evesel: the selected event number. It's =0 by default
50 *
af8e1c2d 51 * Options: Any combination of:
52 * 1) subdetector name: "SPD", "SDD", "SSD", "All" (default)
53 * 2) Printouts: "Verbose" Almost everything is printed out:
54 * it is wise to redirect the output onto a file
813e97e8 55 * e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
af8e1c2d 56 * 3) Rec Points option: "Rec" ---> Uses Rec. Points (default)
a6c9fb45 57 *
58 * 4) ClustersV2 option: "ClustersV2" ---->Uses ClustersV2
59 * otherwise ---->uses hits and digits only
af8e1c2d 60 * Examples:
813e97e8 61 * .x AliITSGeoPlot(); (All subdetectors; no-verbose; no-recpoints)
62 * .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
63 *
88cb7938 64 * filename: It's "galice.root" by default.
813e97e8 65 * isfastpoints: integer. It is set to 0 by defaults. This means that
66 * slow reconstruction is assumed. If fast recpoint are in
67 * in use, isfastpoints must be set =1.
af8e1c2d 68 *
69 * OUTPUT: It produces a root file with a list of histograms
813e97e8 70 *
af8e1c2d 71 * WARNING: spatial information for SSD/DIGITS is obtained by pairing
72 * digits on p and n side originating from the same track, when
73 * possible. This (mis)use of DIGITS is tolerated for debugging
74 * purposes only !!!! The pairing in real life should be done
75 * starting from info really available...
76 *
77 * COMPILATION: this macro can be compiled.
813e97e8 78 * 1) You need to set your include path with
af8e1c2d 79 * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
80 * 3) If you are using root instead of aliroot you need to
81 * execute the macro loadlibs.C in advance
82 * 4) To compile this macro from root (or aliroot):
813e97e8 83 * --- .L AliITSGeoPlot.C++
84 * --- AliITSGeoPlot();
af8e1c2d 85 *
86 * M.Masera 14/05/2001 18:30
a6c9fb45 87 * Last rev. 31/05/2004 14:00 (Clusters V2 added) E.C.
af8e1c2d 88 ********************************************************************/
89
af8e1c2d 90 //Options
91 TString choice(opt);
92 Bool_t All = choice.Contains("All");
93 Bool_t verbose=choice.Contains("Verbose");
94 Bool_t userec=choice.Contains("Rec");
a6c9fb45 95 Bool_t useclustersv2=choice.Contains("ClustersV2");
af8e1c2d 96 Int_t retcode=1; //return code
a6c9fb45 97
af8e1c2d 98 if (gClassTable->GetID("AliRun") < 0) {
a6c9fb45 99 gInterpreter->ExecuteMacro("loadlibs.C");
af8e1c2d 100 }
a6c9fb45 101 else {
af8e1c2d 102 if(gAlice){
88cb7938 103 delete gAlice->GetRunLoader();
af8e1c2d 104 delete gAlice;
105 gAlice=0;
106 }
af8e1c2d 107 }
a6c9fb45 108
109 AliRunLoader* rl = AliRunLoader::Open(filename);
110 if (rl == 0x0){
111 cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
112 return -1;
113 }
114 Int_t retval = rl->LoadgAlice();
115 if (retval){
116 cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
117 return -1;
118 }
119 gAlice=rl->GetAliRun();
88cb7938 120
a6c9fb45 121 retval = rl->LoadHeader();
122 if (retval){
123 cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
124 return -1;
125 }
88cb7938 126
a6c9fb45 127 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
88cb7938 128
a6c9fb45 129 if(!ITSloader){
130 cerr<<"AliITSGeoPlot.C : ITS loader not found"<<endl;
131 return -1;
132 }
88cb7938 133
a6c9fb45 134 ITSloader->LoadHits("read");
135 ITSloader->LoadDigits("read");
136 if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
137 ITSloader->LoadRecPoints("read");
138 rl->GetEvent(evesel);
139 Int_t nparticles = rl->GetHeader()->GetNtrack();
140 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
141 ITS->SetTreeAddress();
af8e1c2d 142 if(verbose) {
143 cout<<" "<<endl<<" "<<endl;
144 cout<<"******* Event processing started *******"<<endl;
145 cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
146 cout<<"will not be processed"<<endl;
147 cout << "Number of particles= " << nparticles <<endl;
148 }
149
150 // HITS
88cb7938 151 TTree *TH = ITSloader->TreeH();
af8e1c2d 152 Stat_t ntracks = TH->GetEntries();
153 if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
154
155 // ITS
af8e1c2d 156 Int_t nmodules;
157 ITS->InitModules(-1,nmodules);
158 cout<<"Number of ITS modules= "<<nmodules<<endl;
159 cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
160 ITS->FillModules(0,0,nmodules," "," ");
161 cout<<"ITS modules .... DONE!"<<endl;
162
163 // DIGITS
88cb7938 164 TTree *TD = ITSloader->TreeD();
af8e1c2d 165
a6c9fb45 166 //RECPOINTS (V2 clusters)
88cb7938 167 TTree *TR = ITSloader->TreeR();
af8e1c2d 168 TClonesArray *ITSrec = ITS->RecPoints();
813e97e8 169 TBranch *branch = 0;
170 if(userec && TR && ITSrec){
171 if(isfastpoints==1){
88cb7938 172 branch = ITSloader->TreeR()->GetBranch("ITSRecPointsF");
813e97e8 173 cout<<"using fast points\n";
174 }
175 else {
88cb7938 176 branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
813e97e8 177 }
178 if(branch)branch->SetAddress(&ITSrec);
179 }
180
181 if(userec && (!TR || !ITSrec || !branch)){
af8e1c2d 182 userec = kFALSE;
813e97e8 183 cout<<"\n ======================================================= \n";
184 cout<<"WARNING: there are no RECPOINTS on this file ! \n";
185 cout<<"======================================================= \n \n";
af8e1c2d 186 }
a6c9fb45 187 if(useclustersv2 && TR && ITSrec){
188
189 branch = ITSloader->TreeR()->GetBranch("Clusters");
190 if(branch)branch->SetAddress(&ITSrec);
191 }
192
193 if(useclustersv2 && (!TR || !ITSrec || !branch)){
194 useclustersv2 = kFALSE;
195 cout<<"\n ======================================================= \n";
196 cout<<"WARNING: there are no CLUSTERSV2 on this file ! \n";
197 cout<<"======================================================= \n \n";
198 }
199
813e97e8 200
af8e1c2d 201 //local variables
202 Int_t mod; //module number
203 Int_t nbytes = 0;
204 Double_t pos[3]; // Global position of the current module
205 Float_t ragdet; // Radius of detector (x y plane)
206 Int_t first; // first module
207 Int_t last; // last module
208 Int_t nrecp; //number of RecPoints for a given module
209
210 //List of histograms
211 TObjArray histos(26,0); // contains the pointers to the histograms
212 // Book histograms SPD
213 TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
214 TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
215 TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
216 TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
217 TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
218 TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
a6c9fb45 219
220 Char_t title1[50]="";
221 Char_t title2[50]="";
222 if(userec){
223 sprintf(title1,"Radii of recpoints - %s","SPD");
224 sprintf(title2,"Z of recpoints - %s","SPD");
225 }
226 if(useclustersv2){
227 sprintf(title1,"Radii of clustersV2 - %s","SPD");
228 sprintf(title2,"Z of clustersV2 - %s","SPD");
229 }
230 TH2F *rrspd = new TH2F("rrspd",title1,50,-10.,10.,50,-10.,10.);
231 TH1F *zrspd = new TH1F("zrspd",title2,100,-30.,30.);
af8e1c2d 232 TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
233 histos.AddLast(rspd); // 0
234 histos.AddLast(rhspd); // 1
235 histos.AddLast(rmspd); // 2
236 histos.AddLast(zspd); // 3
237 histos.AddLast(zhspd); // 4
238 histos.AddLast(zmspd); // 5
239 histos.AddLast(rrspd); // 6
240 histos.AddLast(zrspd); // 7
241 histos.AddLast(enespd); // 8
242 // Book histograms SDD
243 TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
244 TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
245 TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
246 TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
247 TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
248 TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
a6c9fb45 249 Char_t title3[50];
250 Char_t title4[50];
251 if(userec){
252 sprintf(title3,"Radii of recpoints - %s","SDD");
253 sprintf(title4,"Z of recpoints - %s","SDD");
254 }
255 if(useclustersv2){
256 sprintf(title3,"Radii of clustersV2 - %s","SDD");
257 sprintf(title4,"Z of clustersV2 - %s","SDD");
258 }
259 TH2F *rrsdd = new TH2F("rrsdd",title3,50,-40.,40.,50,-40.,40.);
260 TH1F *zrsdd = new TH1F("zrsdd",title4,100,-40.,40.);
af8e1c2d 261 TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
262 histos.AddLast(rsdd); // 9
263 histos.AddLast(rhsdd); // 10
264 histos.AddLast(rmsdd); // 11
265 histos.AddLast(zsdd); // 12
266 histos.AddLast(zhsdd); // 13
267 histos.AddLast(zmsdd); // 14
268 histos.AddLast(rrsdd); // 15
269 histos.AddLast(zrsdd); // 16
270 histos.AddLast(enesdd); // 17
271 // Book histogram SSD
272 TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
273 TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
274 TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
275 TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
276 TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
277 TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
a6c9fb45 278 Char_t title5[50];
279 Char_t title6[50];
280 if(userec){
281 sprintf(title5,"Radii of recpoints - %s","SSD");
282 sprintf(title6,"Z of recpoints - %s","SSD");
283 }
284 if(useclustersv2){
285 sprintf(title5,"Radii of clustersV2 - %s","SSD");
286 sprintf(title6,"Z of clustersV2 - %s","SSD");
287 }
288
289 TH2F *rrssd = new TH2F("rrssd",title5,50,-50.,50.,50,-50.,50.);
290 TH1F *zrssd = new TH1F("zrssd",title6,100,-70.,70.);
af8e1c2d 291 TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
292 histos.AddLast(rssd); // 18
293 histos.AddLast(rhssd); // 19
294 histos.AddLast(rmssd); // 20
295 histos.AddLast(zssd); // 21
296 histos.AddLast(zhssd); // 22
297 histos.AddLast(zmssd); // 23
298 histos.AddLast(rrssd); // 24
299 histos.AddLast(zrssd); // 25
300 histos.AddLast(enessd); // 26
301 //
302 // Loop on subdetectors
303 //
304 AliITSgeom *geom = ITS->GetITSgeom();
305 TString detna; // subdetector name
306 for(Int_t subd=0;subd<3;subd++){
307 if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
308 // Prepare array for the digits
309 TClonesArray *ITSdigits = ITS->DigitsAddress(subd);
310 Bool_t usedigits = kTRUE;
311 if(!ITSdigits){
312 usedigits = kFALSE;
813e97e8 313 cout<<"\n ======================================================= \n";
314 cout<<"WARNING: there are no DIGITS on this file ! \n";
315 cout<<"======================================================= \n \n";
af8e1c2d 316 }
317 // Get segmentation model
318 AliITSDetType *iDetType=ITS->DetType(subd);
319 if(subd==0)detna="SPD";
320 if(subd==1)detna="SDD";
321 if(subd==2)detna="SSD";
322 AliITSsegmentation *seg=(AliITSsegmentation*)iDetType->GetSegmentationModel();
323 // Loop on modules
324 first = geom->GetStartDet(subd);
325 last = geom->GetLastDet(subd);
326 if(verbose){
327 cout<<" "<<endl<<"-------------------------------------"<<endl;
328 cout<<"Start processing subdetector "<<detna<<endl;
329 cout<<detna<<" first module "<<first<<endl;
330 cout<<detna<<" last module "<<last<<endl;
331 cout<<" "<<endl<<" "<<endl;
332 }
333 for (mod=first; mod<=last; mod++){
334 geom->GetTrans(mod,pos); // position of the module in the MRS
335 ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
336 // The following 2 histos are a check of the geometry
337 TH2F *bidi = (TH2F*)histos.At(2+subd*9);
338 TH1F *uni = (TH1F*)histos.At(5+subd*9);
339 bidi->Fill(pos[0],pos[1]);
340 uni->Fill(pos[2]);
341 if(verbose){
813e97e8 342 cout<<"=========================================================\n";
af8e1c2d 343 cout<<detna<<" module="<<mod<<endl;
813e97e8 344 cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", ";
345 cout<<pos[2]<<" Radius "<<ragdet<<endl;
af8e1c2d 346 }
347
348 // Hits
349 GetHitsCoor(ITS,mod,histos,subd,verbose);
350
351 //RecPoints
352 if(userec){
353 ITS->ResetRecPoints();
813e97e8 354 branch->GetEvent(mod);
af8e1c2d 355 TH2F *bidi=(TH2F*)histos.At(6+subd*9);
356 TH1F *uni=(TH1F*)histos.At(7+subd*9);
357 nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
358 }
a6c9fb45 359 if(useclustersv2){
360 ITS->ResetRecPoints();
361 branch->GetEvent(mod);
362 TH2F *bidi=(TH2F*)histos.At(6+subd*9);
363 TH1F *uni=(TH1F*)histos.At(7+subd*9);
364 nrecp=GetClusCoor(geom,ITSrec,mod,bidi,uni,verbose);
365 }
af8e1c2d 366
367 // Digits
368 if(usedigits){
369 ITS->ResetDigits();
370 nbytes += TD->GetEvent(mod);
371 GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
372 }
373
374 } // End of loop on the modules
375 TH1F *h1tmp;
376 TH2F *h2tmp;
377 // Plot the histograms
378 TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
379 if(subd==0){
380 // Prepare canvas 1
381 TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
382 c1->Divide(2,3);
383 current=c1;
384 }
385 if(subd==1){
386 // Prepare canvas 2
387 TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
388 c2->Divide(2,3);
389 current=c2;
390 }
391 if(subd==2){
392 // Prepare canvas 3
393 TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
394 c3->Divide(2,3);
395 current=c3;
396 }
397 current->cd(1);
398 h2tmp = (TH2F*)histos.At(9*subd);
399 h2tmp->Draw();
400 current->cd(2);
401 h1tmp=(TH1F*)histos.At(3+subd*9);
402 h1tmp->Draw();
403 current->cd(3);
404 h2tmp=(TH2F*)histos.At(1+9*subd);
405 h2tmp->Draw();
406 current->cd(4);
407 h1tmp=(TH1F*)histos.At(4+subd*9);
408 h1tmp->Draw();
409
a6c9fb45 410 if(userec || useclustersv2){
af8e1c2d 411 current->cd(5);
412 h2tmp=(TH2F*)histos.At(6+9*subd);
413 h2tmp->Draw();
414 current->cd(6);
415 h1tmp=(TH1F*)histos.At(7+subd*9);
416 h1tmp->Draw();
417 }
418
419 else {
420 current->cd(5);
421 h2tmp=(TH2F*)histos.At(2+9*subd);
422 h2tmp->Draw();
423 current->cd(6);
424 h2tmp=(TH2F*)histos.At(5+9*subd);
425 h2tmp->Draw();
426 }
427 } // if(All.....
428 } // end of loop on subdetectors
429 // Save the histograms
813e97e8 430 TFile *fh = new TFile("AliITSGeoPlot.root","recreate");
af8e1c2d 431 // The list is written to file as a single entry
432 TList *lihi = new TList();
433 // copy the pointers to the histograms to a TList object.
434 // The histograms concerning recpoints are not copied if
435 // 'userec' is false.
436 for(Int_t i=0;i<histos.GetEntriesFast();i++){
437 if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
438 if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
439 }
440 }
441 lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
442 fh->Close();
443
444 return retcode;
445}
446
447
448void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
449 TH2F *h2=(TH2F*)histos.At(1+subd*9);
450 TH1F *h1=(TH1F*)histos.At(4+subd*9);
451 TH1F *ener=(TH1F*)histos.At(8+subd*9);
452 AliITS *ITS= (AliITS*)its;
453 AliITSmodule *modu = ITS->GetModule(mod);
454 TObjArray *fHits = modu->GetHits();
455 Int_t nhits = fHits->GetEntriesFast();
456 if(nhits>0){
457 if(verb){
458 cout<<"-------------------------------------------------------"<<endl;
459 cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
460 }
461 for (Int_t hit=0;hit<nhits;hit++) {
462 AliITShit *iHit = (AliITShit*) fHits->At(hit);
463 if(!iHit->StatusEntering()){
464 Float_t x=iHit->GetXG();
465 Float_t y=iHit->GetYG();
466 Float_t z=iHit->GetZG();
467 Float_t edep=iHit->GetIonization()*1000000;
468 h2->Fill(x,y);
469 h1->Fill(z);
470 ener->Fill(edep);
471 if(verb){
472 cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
813e97e8 473 cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= ";
474 cout<<edep<<endl;
af8e1c2d 475 }
476 }
477 }
478 }
479}
480
481
a6c9fb45 482Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
483 AliITSgeom *geom = (AliITSgeom*)ge;
484 Int_t nrecp = ITSrec->GetEntries();
485 if(nrecp>0){
486 Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
487 Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
488 if(verb){
489 cout<<"-------------------------------------------------------"<<endl;
490 cout<<"Number of CLUSTERS for module "<<mod<<": "<<nrecp<<endl;
491 }
492 for(Int_t irec=0;irec<nrecp;irec++) {
493 AliITSclusterV2 *recp = (AliITSclusterV2*)ITSrec->At(irec);
494 Double_t rot[9];
495 geom->GetRotMatrix(mod,rot);
496 Int_t lay,lad,det;
497 geom->GetModuleId(mod,lay,lad,det);
498 Float_t tx,ty,tz;
499 geom->GetTrans(lay,lad,det,tx,ty,tz);
500
501 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
502 Double_t phi1=TMath::Pi()/2+alpha;
503 if(lay==1) phi1+=TMath::Pi();
504
505 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
506 Float_t r=tx*cp+ty*sp;
507 gc[0]= r*cp - recp->GetY()*sp;
508 gc[1]= r*sp + recp->GetY()*cp;
509 gc[2]=recp->GetZ();
510
511 if(verb){
512 Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
513 cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
514 cout<<r<<endl;
515 cout<<"Associated track "<<recp->GetLabel(0)<<endl;
516 }
517 h2->Fill(gc[0],gc[1]);
518 h1->Fill(gc[2]);
519 }
520 }
521 return nrecp;
522}
af8e1c2d 523Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
524 AliITSgeom *geom = (AliITSgeom*)ge;
525 Int_t nrecp = ITSrec->GetEntries();
526 if(nrecp>0){
527 Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
528 Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
529 if(verb){
530 cout<<"-------------------------------------------------------"<<endl;
531 cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
532 }
533 for(Int_t irec=0;irec<nrecp;irec++) {
534 AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
535 lc[0]=recp->GetX();
536 lc[2]=recp->GetZ();
537 geom->LtoG(mod,lc,gc);
538 if(verb){
813e97e8 539 cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
540 cout<<lc[2]<<endl;
af8e1c2d 541 Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
813e97e8 542 cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
543 cout<<r<<endl;
544 cout<<"Associated track "<<recp->GetLabel(0)<<endl;
af8e1c2d 545 }
546 h2->Fill(gc[0],gc[1]);
547 h1->Fill(gc[2]);
548 }
549 }
550 return nrecp;
551}
552
553void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
554 AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
555 AliITSgeom *geom = (AliITSgeom*)ge;
556 Int_t layer;
557 Int_t ladder;
558 Int_t detec;
559 if(subd==2){
560 geom->GetModuleId(mod,layer,ladder,detec);
561 seg->SetLayer(layer);
562 }
563 Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.; //local coord dig.
564 Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
565 Float_t ragdig; // Radius digit
566 TArrayI ssdone(5000); // used to match p and n side digits of strips
567 TArrayI pair(5000); // as above
568 Int_t ndigits = ITSdigits->GetEntries();
569 AliITSdigit *digs;
570 if(ndigits){
571 if(verbose){
572 cout<<"-------------------------------------------------------"<<endl;
573 cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
574 }
575 // Get the coordinates of the module
576 if(subd==2){
577 for (Int_t digit=0;digit<ndigits;digit++){
578 ssdone[digit]=0;
579 pair[digit]=0;
580 }
581 }
582 for (Int_t digit=0;digit<ndigits;digit++) {
583 digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
0be9848a 584 Int_t iz=digs->GetCoord1(); // cell number z
585 Int_t ix=digs->GetCoord2(); // cell number x
af8e1c2d 586 // Get local coordinates of the element
587 if(subd<2){
588 seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
589 }
590 else{
591 // SSD: if iz==0 ---> N side; if iz==1 P side
592 if(ssdone[digit]==0){
593 ssdone[digit]=1;
594 pair[digit]=-1;
595 Bool_t pside=(iz==1);
596 Bool_t impaired=kTRUE;
597 Int_t pstrip=0;
598 Int_t nstrip=0;
599 if(pside)pstrip=ix;
600 if(!pside)nstrip=ix;
601 for(Int_t digi2=0;digi2<ndigits;digi2++){
602 if(ssdone[digi2]==0 && impaired){
603 AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
0be9848a 604 if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
af8e1c2d 605 ssdone[digi2]=2;
606 pair[digit]=digi2;
0be9848a 607 if(pside)nstrip=dig2->GetCoord2();
608 if(!pside)pstrip=dig2->GetCoord2();
af8e1c2d 609 impaired=kFALSE;
610 }
611 }
612 }
613 if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
614 }
615 }
616 if(subd<2 || (subd==2 && ssdone[digit]==1)){
0be9848a 617 Int_t coor1=digs->GetCoord1();
618 Int_t coor2=digs->GetCoord2();
af8e1c2d 619 Int_t tra0=digs->GetTrack(0);
620 if(verbose){
813e97e8 621 cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
622 cout<<coor2<<" track "<<tra0<<" "<<endl;
623 if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z=";
624 cout<<lcoor[2]<<endl;
af8e1c2d 625 if(subd==2){
626 if(pair[digit]==-1){
627 cout<<"(digit not paired)"<<endl;
628 }
629 else {
813e97e8 630 cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2];
631 cout<<endl;
af8e1c2d 632 Int_t dtmp=pair[digit];
633 AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
0be9848a 634 Int_t coor1b=dig2->GetCoord1();
635 Int_t coor2b=dig2->GetCoord2();
af8e1c2d 636 Int_t tra0b=dig2->GetTrack(0);
637 cout<<"(digit paired with digit #"<<dtmp<<endl;
813e97e8 638 cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
639 cout<<" track "<<tra0b<<")"<<endl;
af8e1c2d 640 }
641 }
642 }
643 if(subd<2 || (subd==2 && pair[digit]!=-1)){
644 // Global coordinates of the element
645 //SDD and SPD use cm, SSD microns (GetPadCxz)
646 if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
647 lcoor[1]=0.;
648 geom->LtoG(mod,lcoor,gcoor); // global coord. in cm
649 ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
544f2385 650 if(verbose){
a6c9fb45 651 cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
652 cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
544f2385 653 }
af8e1c2d 654 //Fill histograms
655 TH2F *bidi = (TH2F*)histos.At(subd*9);
656 TH1F *uni = (TH1F*)histos.At(3+subd*9);
657 bidi->Fill(gcoor[0],gcoor[1]);
658 uni->Fill(gcoor[2]);
659 }
660 }
661 } // loop on digits for this module
662 } // if(ndigits>0....
663}