A more efficient way of volume creations and other minor improvements (Mario)
[u/mrichter/AliRoot.git] / ITS / UPGRADE / testITSUv1 / GetMaterialBudget.C
CommitLineData
1fc6eff6 1#ifndef __CINT__
2
3#include "TSystem.h"
4#include "TGeoManager.h"
5#include <AliRunLoader.h>
6#include <AliGeomManager.h>
7#include <AliITSgeom.h>
8#include <AliTracker.h>
9#include "TRandom.h"
10#include "TMath.h"
11#include "TH2F.h"
12#include "TCanvas.h"
13#include "TFile.h"
14#include <AliRun.h>
15#include <TLegend.h>
16
17#include <TGeoVolume.h>
18#include "TGeoMaterial.h"
19#include "TGeoMedium.h"
20#include <TGeoTube.h>
21#include <AliITSUGeomTGeo.h>
22#include <TPaveText.h>
23#include <TText.h>
24
25#endif
26
27/*
28 gSystem->Load("libITSUpgradeBase");
29 gSystem->Load("libITSUpgradeSim");
30
31 .L GetMaterialBudget.C
32
33 DrawMaterialBudget_SPLITTED();
34 DrawMaterialBudget_FromTo();
35
36 GetMaterialBudget_EtaVsPhi(0,0);
37 GetMaterialBudget_inPhi(0,0);
38
39*/
40
41// Original version
94fc91be 42// Chinorat Kobdaj <kobdaj@g.sut.ac.th> , Stefan Rossegger
1fc6eff6 43// Revised and adapted to newer AliITSUv1 versions
e14c51b6 44// Mario Sitta <sitta@to.infn.it> - Nov 2014, Feb 2015
1fc6eff6 45
46
47enum {
48 kIBModelDummy=0,
49 kIBModel0=1,
50 kIBModel1=2,
51 kIBModel21=3,
52 kIBModel22=4,
53 kIBModel3=5,
e14c51b6 54 kIBModel4=10,
1fc6eff6 55 kOBModelDummy=6,
56 kOBModel0=7,
57 kOBModel1=8,
58 kOBModel2=9
59};
60
61
62//______________________________________________________________________
63void PrintMaterialDefs(const char *geofile="geometry.root")
64{
65 //
66 // Print material definition for all ITS materials/mediums
67 // Rewritten - M.Sitta 15 Nov 2014
68 //
69 TGeoManager::Import(geofile);
70
71 TGeoMedium *med;
72 TGeoMaterial *mat;
73 char mediaName[50], matName[50], shortName[50];
74
75 Int_t nMedia = gGeoManager->GetListOfMedia()->GetEntries();
76
77 printf("\n\n\n");
78 printf(" ==== ITSU Material Properties ====\n\n");
79 printf(" A Z d (g/cm3) RadLen (cm) IntLen (cm)\t Name\n");
80
81 // Loop on media, select ITS materials, print their characteristics
82 for (Int_t i = 0; i < nMedia; i++) {
83 med = (TGeoMedium *)(gGeoManager->GetListOfMedia()->At(i));
84 strcpy(mediaName, med->GetName());
85 // The trick here: the name must begin with the string "ITS_"
86 // so the pointer to the first occurrence of the substring as returned
87 // by strstr() must match the beginning of the string
88 if (strstr(mediaName,"ITS_") == mediaName) { // name begins with "ITS_"
89 mat = med->GetMaterial();
90 strcpy(matName, mat->GetName());
91 strncpy(shortName, matName+4, strlen(matName)-5); // get rid of ITS_ , $
92 shortName[strlen(matName)-5] = '\0';
93 printf(" %5.1f %6.1f %8.3f %13.3f %11.3f\t %s\n",
94 mat->GetA(), mat->GetZ(), mat->GetDensity(), mat->GetRadLen(),
95 mat->GetIntLen(), shortName);
96 }
97 }
98
99 printf("\n");
100
101// printf("Values in Corrados estimate");
102// const Int_t nC=6;
103// TString medNC[nC]={"ITS_SI$","ITS_WATER$","ITS_CARBON$","ITS_KAPTON$","ITS_FLEXCABLE$","ITS_GLUE$"};
104// Double_t radl[nC]={9.36,36.1,25,28.6,13.3,44.37};
105// for (Int_t i=0; i<nC; i++) {
106// printf("%8.2lf (cm) \t%s\n", radl[i],medNC[i].Data());
107// }
108
109}
110
111
112//______________________________________________________________________
113Double_t ComputeMaterialBudget(Double_t rmin, Double_t rmax, Double_t etaPos,
114 Double_t phiMax,
115 TH1F* xOverX01d, TH2F* xOverX02d = 0)
116{
117 //
118 // Ancillary function to compute material budget between rmin and rmax
119 // and +/- etaPos in the 0--phiMax range, filling two histograms
120 // Rewritten - M.Sitta 15 Nov 2014
121 //
122 Int_t n = 5e4; //testtracks
123
124 TH1F *nParticles1d = new TH1F("", "", 300, 0, phiMax);
125 TH2F *nParticles2d = new TH2F("", "", 300, -etaPos, etaPos, 300, 0, phiMax);
126
127 Double_t mparam[7]={0.,0.,0.,0.,0.};
128 Double_t point1[3],point2[3];
129
130 for (Int_t it=0; it<n; it++) {
131
132 // PHI VS ETA
133 Double_t phi = gRandom->Rndm()*phiMax;
134 Double_t eta = gRandom->Rndm()*2*etaPos - etaPos;
135 Double_t theta = TMath::ATan(TMath::Exp(-eta))*2;
136 point1[0] = rmin*TMath::Cos(phi);
137 point1[1] = rmin*TMath::Sin(phi);
138 point1[2] = rmin/TMath::Tan(theta);
139 point2[0] = rmax*TMath::Cos(phi);
140 point2[1] = rmax*TMath::Sin(phi);
141 point2[2] = rmax/TMath::Tan(theta);
142
143 AliTracker::MeanMaterialBudget(point1,point2,mparam);
144
145 xOverX01d->Fill(phi, mparam[1]*100);
146 nParticles1d->Fill(phi, 1.);
147
148 if (xOverX02d != 0) {
149 xOverX02d->Fill(eta, phi, mparam[1]*100);
150 nParticles2d->Fill(eta, phi, 1.);
151 }
152
153 if (!(it%10000)) cout<<" : "<<mparam[1]*100<<" phi:"<<phi<<endl;
154
155 }
156
157 // normalization to number of particles in case of phi vs eta
158 Double_t theta = TMath::ATan(TMath::Exp(-etaPos/2))*2;
159 printf("<eta>=%lf -> Sin(theta) %lf\n",etaPos/2,TMath::Sin(theta));
160
161 for (Int_t ix = 1; ix<=nParticles1d->GetNbinsX(); ix++) {
162 if (nParticles1d->GetBinContent(ix) > 0)
163 xOverX01d->SetBinContent(ix,xOverX01d->GetBinContent(ix)/nParticles1d->GetBinContent(ix)*TMath::Sin(theta));
164 }
165
166 if (xOverX02d) {
167 for (Int_t ix = 1; ix<=nParticles2d->GetNbinsX(); ix++) {
168 for (Int_t iy = 1; iy<=nParticles2d->GetNbinsY(); iy++) {
169 if (nParticles2d->GetBinContent(ix,iy) > 0)
170 xOverX02d->SetBinContent(ix,iy,xOverX02d->GetBinContent(ix,iy)/nParticles2d->GetBinContent(ix,iy));
171 }
172 }
173 }
174
175 Double_t mean = 0;
176 for (Int_t ix = 1; ix<=xOverX01d->GetNbinsX(); ix++)
177 mean+=xOverX01d->GetBinContent(ix);
178
179 mean /= xOverX01d->GetNbinsX();
180
181 return mean;
182}
183
184
185//______________________________________________________________________
186void DrawMaterialBudget_Splitted(Int_t nLay = 0, Double_t rmin = 1.,
187 Double_t rmax = 5.)
188{
189 //
190 // Function to factorize the percentage of each medium in the
191 // total material budget between rmin and rmax for layer nLay
192 // Cleaned up - M.Sitta 15 Nov 2014
193 //
194 TCanvas *c2 = new TCanvas("c2","c2");
195
196 Double_t etaPos = 1.;
197
198 Bool_t firstPlot = 1;
199
200 // Check parameters
201 if (nLay < 0 || nLay > 6) {
202 printf("ERROR! Wrong layer number %d\n",nLay);
203 return;
204 }
205
206 if (rmin < 0. || rmax < 0. || rmax < rmin) {
207 printf("ERROR! Wrong radial values rmin = %f rmax = %f\n",rmin,rmax);
208 return;
209 }
210
211 // In current version, build levels have different meaning for IB and OB
212 const Int_t nScenariosIB = 6, nScenariosOB = 7;
213 TString strDescripOB[nScenariosOB] = {"Copper", "Aluminum", "Glue", "Water",
214 "Kapton", "Carbon" , "Silicon" };
215 // Color codes
216 // Water kBlack
217 // Kapton kYellow
218 // Carbon kRed
219 // Glue kBlue
220 // Aluminum kCyan
221 // Si-Sensor kGreen
222 // Copper kGray
223 Int_t colorsOB[nScenariosOB] = {kGray , kCyan, kBlue, kBlack,
224 kYellow, kRed , kGreen };
225
226// TString strDescripIB[nScenariosIB] = {"Carbon Structure", "Water",
227// "Cooling Pipe Walls and ColdPlate",
228// "Glue", "Flex Cable", "Pixel Chip"};
229 TString strDescripIB[nScenariosIB] = {"Flex cable", "Glue",
230 "Carbon structure", "Water",
231 "Cooling walls", "Pixel Chip"};
232 // Color codes
233 // Flex Cable kBlack
234 // Glue kYellow
235 // Carbon Structure kRed
236 // Water kBlue
237 // Cooling Walls kCyan
238 // Si-Sensor kGreen
239 Int_t colorsIB[nScenariosIB] = {kCyan, kBlue, kBlack, kYellow, kRed, kGreen};
240
241 // Setup
242 const Int_t maxNumberOfScenarios = nScenariosOB;
243
244 Double_t meanX0[maxNumberOfScenarios];
e14c51b6 245 Double_t contribX0[maxNumberOfScenarios];
1fc6eff6 246 TH1F *l[maxNumberOfScenarios+1];
247
248 TString strDescrip[maxNumberOfScenarios];
249 Int_t colors[maxNumberOfScenarios];
250
251 // Choose which scenario based on Layer number and Model
252 Int_t nScenarios;
253
254 TGeoManager::Import(Form("geometry_0.root"));
255 AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
256 Int_t nLad = gm->GetNStaves(nLay);
257
258 Float_t phiMax = TMath::TwoPi()/nLad*2;
259
260 char title[30];
261 Int_t model, buildlevel;
e14c51b6 262
263 strcpy(title,
264 gGeoManager->GetVolume(Form("ITSULayer%d",nLay))->GetTitle());
265 if (strlen(title) == 0) // No Model info -> old model schema
266 model = -1;
267 else // New model schema -> extract model
268 sscanf(title, "Model %d Build level %d", &model, &buildlevel);
269
1fc6eff6 270 if (nLay < 3) { // IB has only 6 scenarios
271 nScenarios = nScenariosIB;
e14c51b6 272 if (model == kIBModel4) // New model -> same name/colors as OB
273 for (Int_t i=0; i<nScenarios; i++) {
274 strDescrip[i] = strDescripOB[i+1];
275 colors[i] = colorsOB[i+1];
276 }
277 else // Old model -> old names/colors
278 for (Int_t i=0; i<nScenarios; i++) {
279 strDescrip[i] = strDescripIB[i];
280 colors[i] = colorsIB[i];
281 }
1fc6eff6 282 } else {
e14c51b6 283 if (model == kOBModel2) { // New OB
284 nScenarios = nScenariosOB;
285 for (Int_t i=0; i<nScenarios; i++) {
286 strDescrip[i] = strDescripOB[i];
287 colors[i] = colorsOB[i];
288 }
289 } else { // Old OB has only 6 scenarios like IB
1fc6eff6 290 nScenarios = nScenariosIB;
291 for (Int_t i=0; i<nScenarios; i++) {
292 strDescrip[i] = strDescripOB[i+1];
293 colors[i] = colorsOB[i+1];
294 }
1fc6eff6 295 }
296 } // if (nLay < 3)
297
298 delete gGeoManager;
299
300 for (Int_t i=0; i<nScenarios; i++) {
301
302 printf(" -> Loading geometry_%d.root .....\n",i);
303 TGeoManager::Import(Form("geometry_%d.root",i));
304
305 strcpy(title,
306 gGeoManager->GetVolume(Form("ITSULayer%d",nLay))->GetTitle());
307 if (strlen(title) != 0) {
308 sscanf(title, "Model %d Build level %d", &model, &buildlevel);
309 if (i != buildlevel)
310 printf("WARNING! Possible mismatch: file geometry_%d.root created with Build level %d\n",i,buildlevel);
311 }
312
313 TH1F *xOverX01d = new TH1F("", "", 300, 0, phiMax);
314
315 Double_t mean = ComputeMaterialBudget(rmin, rmax, etaPos, phiMax,
316 xOverX01d);
317 meanX0[i] = mean;
318 cout<<"Mean X/X0: " << meanX0[i] << " %" << endl;
319
320 xOverX01d->GetXaxis()->SetTitle("#phi (rad)");
321 xOverX01d->GetYaxis()->SetTitle(Form("X/X_{0} (%%) at #eta=0 "));
322 xOverX01d->SetFillColor(colors[i]-7);
323 if (i==0) xOverX01d->SetFillColor(colors[i]);
324 if ( (nScenarios == nScenariosIB && i==2) ||
325 (nScenarios == nScenariosOB && i==3) )
326 xOverX01d->SetFillColor(colors[i]);
327 xOverX01d->SetStats(0);
328
329 l[i] = new TH1F("","",1,0,phiMax);
330 l[i]->SetBinContent(1,mean);
331 l[i]->SetStats(0);
332 l[i]->SetLineColor(colors[i]-7);
333 if (i==0) l[i]->SetLineColor(kGray);
334 if ( (nScenarios == nScenariosIB && i==2) ||
335 (nScenarios == nScenariosOB && i==3) ) l[i]->SetLineColor(12);
336
337 c2->cd();
338 if (firstPlot) {
339 xOverX01d->SetMinimum(0.);
340 xOverX01d->DrawCopy();
341 firstPlot=0;
342 } else
343 xOverX01d->DrawCopy("same");
344
345 delete gGeoManager;
346 }
347
348 // Build meaningful legend
349 TLegend *leg = new TLegend(0.76,0.77,0.99,0.99,"");
350 leg->SetFillColor(0);
351
352 for (Int_t i=0; i<nScenarios; i++) {
353 // contribution in percent
354 Double_t contr = 0;
e14c51b6 355 if (i == nScenarios-1) {
356 contr = (meanX0[i])/meanX0[0]*100;
357 contribX0[i] = meanX0[i];
358 } else {
359 contr = (meanX0[i]-meanX0[i+1])/meanX0[0]*100;
360 contribX0[i] = meanX0[i]-meanX0[i+1];
361 }
1fc6eff6 362
363 strDescrip[i].Append(Form(" (%3.1lf%%)",contr));
364 leg->AddEntry(l[i],strDescrip[i].Data(),"l");
365 }
366
367 TPaveText *pt = new TPaveText(0.76,0.70,0.99,0.77,"brNDC");
368 pt->SetBorderSize(1); // no shadow
369 pt->SetTextFont(12);
370 pt->SetFillColor(0);
371 pt->AddText(Form("Mean X/X0 = %4.3lf%%",meanX0[0]));
372 pt->Draw();
373
374 leg->Draw();
375
376 l[nScenarios]=(TH1F*)l[0]->Clone();
377 l[nScenarios]->SetLineColor(1);
378 l[nScenarios]->SetLineWidth(2);
379 l[nScenarios]->DrawCopy("same");
380
e14c51b6 381 printf("\n X/X0 contributions:\n");
382 for (Int_t i=0; i<nScenarios; i++)
383 printf("%s\t\t%f\n",strDescrip[i].Data(),contribX0[i]);
384
1fc6eff6 385 c2->SaveAs("Material-details.pdf");
386
387}
388
389
390//______________________________________________________________________
391void DrawMaterialBudget_FromTo(Int_t nLay = 0, Double_t rmin = 1.,
392 Double_t rmax = 5., Bool_t only2D = 1)
393{
394 //
395 // Function to compute material budget between rmin and rmax
396 // for layer nLay. If only2D is false a 1D histogram is also plotted
397 // Cleaned up and simplified - M.Sitta 18 Nov 2014
398 //
399
400 Double_t etaPos = 0.25;
401
402 // Check parameters
403 if (nLay < 0 || nLay > 6) {
404 printf("ERROR! Wrong layer number %d\n",nLay);
405 return;
406 }
407
408 if (rmin < 0. || rmax < 0. || rmax < rmin) {
409 printf("ERROR! Wrong radial values rmin = %f rmax = %f\n",rmin,rmax);
410 return;
411 }
412
413
414 TGeoManager::Import("geometry.root");
415 AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
416 Int_t nLad = gm->GetNStaves(nLay);
417
418 Float_t phiMax = TMath::TwoPi()/nLad*2;
419
420 TH2F *xOverX0 = new TH2F("", "", 300, -1, 1., 300, 0, phiMax);
421 TH1F *xOverX01d = new TH1F("", "", 300, 0, phiMax);
422
423 Double_t mean = ComputeMaterialBudget(rmin, rmax, -1, phiMax,
424 xOverX01d, xOverX0);
425 cout<<"Mean X/X0: " << mean << " %" << endl;
426
427 TH1F *l = new TH1F("", "", 1, 0, phiMax);
428 l->SetBinContent(1, mean);
429 l->SetLineColor(2);
430 l->SetStats(0);
431 l->SetLineWidth(2);
432
433 // Draw the histograms
434 xOverX0->SetTitle(0);
435 xOverX0->GetXaxis()->SetTitle("pseudorapidity #eta ");
436 xOverX0->GetYaxis()->SetTitle("#phi (rad)");
437 xOverX0->GetZaxis()->SetTitle("X/X_{0} (%)");
438 xOverX0->SetStats(0);
439
440 xOverX01d->SetTitle(Form("X/X_{0} (%%) within r=(%2.1lf-%2.1lf)cm average over #eta=(%1.1lf,%1.1lf)", rmin, rmax, -1., 1.));
441 xOverX01d->GetXaxis()->SetTitle("#phi (rad)");
442 xOverX01d->GetYaxis()->SetTitle("X/X_{0} (%) average over #eta=(-1,1)");
443 xOverX01d->SetStats(0);
444
445 if (!only2D) {
446 TCanvas *c1 = new TCanvas("c1","c1",800,800);
447 c1->Divide(1,2);
448 c1->cd(1); xOverX0->Draw("colz");
449 c1->cd(2); xOverX01d->DrawCopy(); l->DrawCopy("same");
450 c1->SaveAs("Material-1D.pdf");
451 } else {
452 TCanvas *c1 = new TCanvas("c1","c1");
453 xOverX0->Draw("colz");
454 c1->SaveAs("Material-2D.pdf");
455 }
456
457}
458
459
460//______________________________________________________________________
461void ComputeGeometricalParameters(const char *geofile, Double_t *rmin,
462 Double_t *rmax, Double_t *zPos,
463 Int_t *nlad, Bool_t fromGDML = 0)
464{
465 //
466 // Ancillary function to retrieve various geometrical parameters
467 // from a geometry file. The caller must ensure that all vectors
468 // have proper dimensions: since this function is designed to be used
469 // internally, no check is made on parameters!
470 // Rewritten - M.Sitta 20 Nov 2014
471 //
472
473 if (fromGDML) // from GDML
474 TGeoManager::Import(geofile);
475// gGeoManager->ViewLeaves(true);
476// gGeoManager->GetTopVolume()->Draw("ogl");
477 else { // from AliRoot simulation
478
479// gAlice=NULL;
480// AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
481// runLoader->LoadgAlice();
482
483// gAlice = runLoader->GetAliRun();
484// AliGeomManager::LoadGeometry(geofile);
485 TGeoManager::Import(geofile);
486 AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
487
488 TGeoVolume *pipeV = gGeoManager->GetVolume("IP_PIPE");
489 if (!pipeV) {
490 printf("Pipe volume %s is not in the geometry\n", "IP_PIPE");
491 return;
492 } else {
493 printf("%s ","IP_PIPE");
494 TGeoTube *t = (TGeoTube*)(pipeV->GetShape());
495 printf(" r = (%3.2lf %3.2lf) ", t->GetRmin(), t->GetRmax());
496 printf(" z = %3.2lf\n", t->GetDz());
497
498 rmin[0] = t->GetRmin();
499 rmax[0] = t->GetRmax();
500 nlad[0] = 0;
501 zPos[0] = 0;
502 }
503
504 TGeoVolume *itsV = gGeoManager->GetVolume(gm->GetITSVolPattern());
505 if (!itsV) {
506 printf("ITS volume %s is not in the geometry\n", gm->GetITSVolPattern());
507 return;
508 }
509 //
510 // Loop on all ITSV nodes, count Layer volumes by checking names
511 Int_t nNodes = itsV->GetNodes()->GetEntries();
512 Int_t numberOfLayers = 0;
513 for (Int_t j=0; j<nNodes; j++) {
514 if ( strstr(itsV->GetNodes()->At(j)->GetName(),
515 gm->GetITSWrapVolPattern()) ) {
516
517 TGeoNode *wrapN = (TGeoNode*)(itsV->GetNodes()->At(j));
518 TGeoVolume *wrapV = wrapN->GetVolume();
519 Int_t nNodesWrap = wrapV->GetNodes()->GetEntries();
520
521 for (Int_t k=0; k<nNodesWrap; k++) {
522 if ( strstr(wrapV->GetNodes()->At(k)->GetName(),
523 gm->GetITSLayerPattern()) ) {
524
525 Int_t l = numberOfLayers;
526 numberOfLayers++;
527
528 char laynam[30];
529 snprintf(laynam, 30, "%s%d", gm->GetITSLayerPattern(), l);
530 TGeoVolume* volLy = gGeoManager->GetVolume(laynam);
531 printf("%s\t",volLy->GetName());
532
533 // Layers are Assemblies not Tubes, so:
534 // for rmax we assume the maximum extension of the assembly
535 // for rmin: for OB (no Turbo) we take the stave height and
536 // we subtract this from rmax; for IB (Turbo) we compute the
537 // position of the most internal corner
538 TGeoBBox *b = (TGeoBBox*)(volLy->GetShape());
539 TGeoBBox *s = (TGeoBBox*)(gGeoManager->GetVolume(Form("%s%d",gm->GetITSStavePattern(),l))->GetShape());
540 Double_t minr;
541
542 Double_t loc[3], mas[3];
543 if (l < 3) {
544 Double_t loc[3], mas[3];
545 loc[0] = s->GetDX();
546 loc[1] = s->GetDY();
547 loc[2] = 0;
548 volLy->GetNode(Form("%s%d_0",gm->GetITSStavePattern(),l))->GetMatrix()->LocalToMaster(loc,mas);
549 minr = TMath::Sqrt(mas[0]*mas[0] + mas[1]*mas[1]);
550 zPos[l+1] = 0;
551 } else {
552 minr = b->GetDX() - 2*s->GetDX();
553 TGeoBBox *c = (TGeoBBox*)(gGeoManager->GetVolume(Form("%s%d",gm->GetITSChipPattern(),l))->GetShape());
554 zPos[l+1] = c->GetDZ();
555 }
556 rmin[l+1] = minr;
557 rmax[l+1] = b->GetDX();
558 nlad[l+1] = gm->GetNStaves(l);
559
560 printf(" r = (%5.2lf , %5.2lf) ", rmin[l+1], rmax[l+1]);
561 printf(" z = +/- %5.2lf ", b->GetDZ());
562 printf(" #lad = %d \n", nlad[l+1]);
563
564 }
565 }
566 }
567 }
568 }
569
570}
571
572
573//______________________________________________________________________
574void DrawMaterialBudget_inPhi(const char *geofile = "geometry.root",
575 Int_t lay = -1, Bool_t fromGDML = 0)
576{
577 //
578 // Function to compute material budget as a function of phi
579 // for layer lay (all layers and pipe if lay=-1). geofile is the
580 // geometry file name. If fromGDML is true, use a GDML file, else
581 // a standard geometry.root file.
582 // Cleaned up and simplified - M.Sitta 18 Nov 2014
583 //
584
585 Double_t rmin[8], rmax[8], zPos[8];
586 Int_t nlad[8];
587
588 ComputeGeometricalParameters(geofile, rmin, rmax, zPos, nlad, fromGDML);
589
590 Int_t n = 100000; //testtracks
591 // Test Get material ?
592 Double_t mparam[7]={0.,0.,0.,0.,0.};
593 Double_t point1[3],point2[3];
594
595 TH2F *xOvsPhi[8];
596 Int_t ls = 0, le = 8;
597 if (lay != -1) {
598 ls = lay+1;
599 le = lay+2;
600 }
601
602 for(Int_t i=ls; i<le; i++) {
603 Double_t phiMax;
604 if (i == 0) // PIPE doesn't have staves
605 phiMax = TMath::TwoPi()/4.;
606 else
607 phiMax = TMath::TwoPi()/nlad[i]*5.; // show approx 5 ladders
608
609// xOvsPhi[i] = new TH2F("", "", 100, 0, phiMax, 100, 0.2, 0.8);
610 if (i < 4)
611 xOvsPhi[i] = new TH2F("", "", 100, 0, phiMax, 100, 0.0, 1.2);
612 else
613 xOvsPhi[i] = new TH2F("", "", 100, 0, phiMax, 100, 0.2, 2.2);
614
615 for (Int_t it=0; it<n; it++) {
616 Double_t phi = phiMax*gRandom->Rndm();
617 Double_t z = zPos[i];
618 point1[0] = rmin[i]*TMath::Cos(phi);
619 point1[1] = rmin[i]*TMath::Sin(phi);
620 point1[2] = z;
621 point2[0] = rmax[i]*TMath::Cos(phi);
622 point2[1] = rmax[i]*TMath::Sin(phi);
623 point2[2] = z;
624 AliTracker::MeanMaterialBudget(point1,point2,mparam);
625 if (mparam[1] < 100) // don't save fakes due to errors
626 xOvsPhi[i]->Fill(phi,mparam[1]*100); //fxOverX0Layer
627 if (!(it%10000)) cout << "layer" << i-1 << " : " << mparam[1]*100
628 << " r=(" <<rmin[i] << "," << rmax[i] << ")"
629 << " phi=" << phi <<endl;
630 }
631
632 }
633
634 if (lay==-1) {
635 TCanvas *c1 = new TCanvas("c1","Material Budget",700,800);
636 c1->Divide(2,4);
637 for(Int_t i=0; i<8; i++) {
638 c1->cd(i+1);
639 if (i>0)
640 xOvsPhi[i]->SetTitle(Form("Layer %d",i-1));
641 else
642 xOvsPhi[i]->SetTitle(Form("Beampipe"));
643 xOvsPhi[i]->GetXaxis()->SetTitle("phi (rad)");
644 xOvsPhi[i]->GetYaxis()->SetTitle("X/X_{0} (%)");
645 xOvsPhi[i]->SetStats(0);
646 xOvsPhi[i]->Draw("col");
647 }
648 } else {
649 TCanvas *c1 = new TCanvas("c1","Material Budget");
650 Int_t i = lay+1;
651 xOvsPhi[i]->SetTitle(Form("Layer %d",lay));
652 xOvsPhi[i]->GetXaxis()->SetTitle("phi (rad)");
653 xOvsPhi[i]->GetYaxis()->SetTitle("X/X_{0} (%)");
654 xOvsPhi[i]->SetStats(1111111);
655 xOvsPhi[i]->Draw("col");
656 }
657
658}
659
660
661//______________________________________________________________________
662void DrawMaterialBudget_EtaVsPhi(const char *geofile = "geometry.root",
663 Int_t lay = -1, Bool_t fromGDML = 0)
664{
665 //
666 // Function to compute material budget as a function of eta and phi
667 // in eta +/- 1 for layer lay (all layers and pipe if lay=-1).
668 // geofile is the geometry file name. If fromGDML is true, use a
669 // GDML file, else a standard geometry.root file.
670 // Cleaned up and simplified - M.Sitta 20 Nov 2014
671 //
672
673 Double_t rmin[10], rmax[10], zPos[10]; // zPos not used but needed for call
674 Int_t nlad[10];
675
676 ComputeGeometricalParameters(geofile, rmin, rmax, zPos, nlad, fromGDML);
677
678 rmin[8] = rmin[1];
679 rmax[8] = rmax[3];
680 rmin[9] = rmin[4];
681 rmax[9] = rmax[7];
682
683 Int_t n = 100000; //testtracks
684 // Test Get material ?
685 Double_t mparam[7]={0.,0.,0.,0.,0.};
686 Double_t point1[3],point2[3];
687
688 TH2F *xOverX0[10];
689 TH2F *nParticles[10];
690 Int_t ls = 0, le = 10;
691 if (lay != -1) {
692 ls = lay+1;
693 le = lay+2;
694 }
695
696 for(Int_t i=ls; i<le; i++) {
697 Double_t phiMax;
698 if (i == 0) // PIPE doesn't have staves
699 phiMax = TMath::TwoPi()/4.;
700 else if (i == 8) // Mean value of layer 0 to 2
701 phiMax = TMath::TwoPi()/(0.5*(nlad[1]+nlad[3]))*5.;
702 else if (i == 9) // Mean value of layer 3 to 6
703 phiMax = TMath::TwoPi()/(0.5*(nlad[4]+nlad[7]))*5.;
704 else
705 phiMax = TMath::TwoPi()/nlad[i]*5.; // show approx 5 ladders
706
707 xOverX0[i] = new TH2F("", "", 100, -1., 1., 100, 0., phiMax);
708 nParticles[i] = new TH2F("", "", 100, -1., 1., 100, 0., phiMax);
709
710 for (Int_t it=0; it<n; it++) {
711 Double_t phi = phiMax*gRandom->Rndm();
712 Double_t eta = gRandom->Rndm()*2 - 1.; // +/- 1 eta
713 Double_t theta = TMath::ATan(TMath::Exp(-eta))*2;
714 point1[0] = rmin[i]*TMath::Cos(phi);
715 point1[1] = rmin[i]*TMath::Sin(phi);
716 point1[2] = rmin[i]/TMath::Tan(theta);
717 point2[0] = rmax[i]*TMath::Cos(phi);
718 point2[1] = rmax[i]*TMath::Sin(phi);
719 point2[2] = rmax[i]/TMath::Tan(theta);
720
721 AliTracker::MeanMaterialBudget(point1,point2,mparam);
722 if (mparam[1] < 100) { // don't save fakes due to errors
723 xOverX0[i]->Fill(eta,phi,mparam[1]*100); //fxOverX0Layer
724 nParticles[i]->Fill(eta,phi,1.);
725 }
726 if (!(it%10000)) cout << "layer" << i-1 << " : " << mparam[1]*100
727 << " r=(" <<rmin[i] << "," << rmax[i] << ")"
728 << " phi=" << phi << " theta=" << theta <<endl;
729 }
730
731 // normalization to number of particles
732 for (Int_t ix=1; ix<=nParticles[i]->GetNbinsX(); ix++) {
733 for (Int_t iy=1; iy<=nParticles[i]->GetNbinsY(); iy++) {
734 if (nParticles[i]->GetBinContent(ix,iy) > 0)
735 xOverX0[i]->SetBinContent(ix,iy,xOverX0[i]->GetBinContent(ix,iy)/nParticles[i]->GetBinContent(ix,iy));
736 }
737 }
738 }
739
740 if (lay==-1) {
741 TCanvas *c1 = new TCanvas("c1","Material Budget",750,900);
742 c1->Divide(2,5);
743 for(Int_t i=0; i<10; i++) {
744 c1->cd(i+1);
745 if (i==0)
746 xOverX0[i]->SetTitle(Form("Beampipe"));
747 else if (i==8)
748 xOverX0[i]->SetTitle(Form("Layer 0 to 2"));
749 else if (i==9)
750 xOverX0[i]->SetTitle(Form("Layer 3 to 6"));
751 else
752 xOverX0[i]->SetTitle(Form("Layer %d",i-1));
753
754
755 xOverX0[i]->SetStats(0);
756 xOverX0[i]->GetXaxis()->SetTitle("pseudorapidity #eta ");
757 xOverX0[i]->GetYaxis()->SetTitle("#phi (rad)");
758 xOverX0[i]->GetZaxis()->SetTitle("X/X_{0} (%)");
759 xOverX0[i]->Draw("colz");
760 }
761 } else {
762 TCanvas *c1 = new TCanvas("c1","Material Budget");
763 Int_t i = lay+1;
764 xOverX0[i]->SetTitle(Form("Layer %d",lay));
765 xOverX0[i]->SetStats(0);
766 xOverX0[i]->GetXaxis()->SetTitle("pseudorapidity #eta ");
767 xOverX0[i]->GetYaxis()->SetTitle("#phi (rad)");
768 xOverX0[i]->GetZaxis()->SetTitle("X/X_{0} (%)");
769 xOverX0[i]->Draw("colz");
770 }
771
772}