Update master to aliroot
[u/mrichter/AliRoot.git] / ITS / UPGRADE / testITSUv1 / GetMaterialBudget.C
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
42 // Chinorat Kobdaj <kobdaj@g.sut.ac.th> , Stefan Rossegger
43 // Revised and adapted to newer AliITSUv1 versions
44 // Mario Sitta <sitta@to.infn.it> - Nov 2014, Feb 2015
45
46
47 enum {
48     kIBModelDummy=0,
49     kIBModel0=1,
50     kIBModel1=2, 
51     kIBModel21=3,
52     kIBModel22=4,
53     kIBModel3=5,
54     kIBModel4=10,
55     kOBModelDummy=6,
56     kOBModel0=7,
57     kOBModel1=8, 
58     kOBModel2=9 
59 };
60
61
62 //______________________________________________________________________
63 void 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 //______________________________________________________________________
113 Double_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 //______________________________________________________________________
186 void 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];
245   Double_t contribX0[maxNumberOfScenarios];
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;
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
270   if (nLay < 3) { // IB has only 6 scenarios
271     nScenarios = nScenariosIB;
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       }
282   } else {
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
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       }
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;
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     }
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   
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
385   c2->SaveAs("Material-details.pdf");
386   
387 }
388
389
390 //______________________________________________________________________
391 void 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 //______________________________________________________________________
461 void 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 //______________________________________________________________________
574 void 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 //______________________________________________________________________
662 void 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 }