1 void testEMCALGeom(const char *filename="galice.root",){
2 /////////////////////////////////////////////////////////////////////////
3 // This macro computes the sampling fraction in the EMCAL
5 // Root > .L ComputeSamplingFraction.C //this loads the macro in memory
6 // Root > ComputeSamplingFraction(); //by default process first event
7 // Root > ComputeSamplingFraction("galice2.root"); //process from the file
11 <img src="picts/ComputeSamplingFraction.gif">
14 /////////////////////////////////////////////////////////////////////////
19 // Dynamically link some shared libs
20 if(gClassTable->GetID("AliRun") < 0) {
21 gROOT->LoadMacro("loadlibs.C");
25 // Connect the Root Galice file containing Geometry, Kine and Hits
26 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
27 if(!file) file = new TFile(filename);
29 // Get AliRun object from file or create it if not on file
31 gAlice = (AliRun*)file->Get("gAlice");
32 if(gAlice) printf("AliRun object found on file\n");
33 if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
36 // Set event pointer to this event
37 Int_t nparticles = gAlice->GetEvent(0);
39 cout << "No particles found for event " << evNumber;
40 cout << " in file " << filename << endl;
42 } // end if nparticles <=0
44 // Pointer to specific detector hits.
45 AliEMCALHit *emcalHit;
47 // Get pointers to ALL Alice detectors and Hits containers
48 AliEMCALv0 *EMCAL = (AliEMCALv0*) gAlice->GetDetector("EMCAL");
50 cout << "EMCAL not Found. Exiting." << endl;
54 AliEMCALGeometry *g = EMCAL->GetGeometry();
56 Int_t i,in,ix,ieta,iphi,ipre;
58 cout << "i ieta iphi ipre in eta phi ix" << endl;
59 for(i=0;i<=2*g->GetNEta()*g->GetNPhi();i++){
60 g->TowerIndexes(i,ieta,iphi,ipre);
61 in = g->TowerIndex(ieta,iphi,ipre);
62 g->EtaPhiFromIndex(i,eta,phi);
63 ix = g->TowerIndexFromEtaPhi(eta,phi);
64 cout << i << " " << ieta << " " << iphi << " " << ipre;
65 cout << " " << in << " " << eta << " " << phi << " " << ix << endl;
68 // Get pointer to the particle
73 TH1F *hEMCAL = new TH1F("hEMCAL" ,"Energy",100,0.,2.);
74 TH2F *hSfs = new TH2F("hSfs","Sampling Fraction",200,-0.7,0.7,
76 TH2F *hSfi = new TH2F("hSfi","Sampling Fraction Normilization",
77 200,-0.7,0.7,200,0.0,120.0);
79 Int_t track,ntracks = gAlice->TreeH()->GetEntries();
83 Float_t eta,phi,x,y,z;
84 // Start loop on tracks in the hits containers
85 for(track=0; track<ntracks;track++){
88 gAlice->TreeH()->GetEvent(track);
89 for(emcalHit=(AliEMCALHit*)EMCAL->FirstHit(-1);emcalHit;
90 emcalHit=(AliEMCALHit*)EMCAL->NextHit()) {
91 ipart = emcalHit->GetTrack();
92 part = gAlice->Particle(ipart);
93 if(part->GetFirstMother()!=0) continue;
94 // only those that don't interact.
95 hEMCAL->Fill((Float_t)(emcalHit->GetEnergy()));
99 findetaphi(x,y,z,eta,phi);
101 egam += part->Energy();
102 hSfi->Fill(eta,phi,part->Energy());
105 etot += emcalHit->GetEnergy();
106 hSfs->Fill(eta,phi,emcalHit->GetEnergy());
111 cout << "Energy of gamma =" << egam << "GeV. Energy in EMCAL =" << etot;
112 cout << "GeV. Sampling Fractions is therefore f=" << etot/egam << endl;
114 //Create a canvas, set the view range, show histograms
115 TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700);
118 hEMCAL->SetFillColor(42);
124 // c0->SaveAs("analHitsEMCAL.eps");
126 void findetaphi(Float_t x,Float_t y,Float_t z,Float_t &eta,Float_t &phi){
127 // Compute Eta and Phi from x, y, z. Assumes vertex is a zero.
132 phi = 180.*TMath::ATan2(y,x)/TMath::Pi();
133 if(phi<0.0) phi += 180.;
134 r = TMath::ATan2(r,z);
135 r = TMath::Tan(0.5*r);
136 eta = -TMath::Exp(r);