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