]>
Commit | Line | Data |
---|---|---|
1 | void testEMCALGeom(const char *filename="galice.root",){ | |
2 | ///////////////////////////////////////////////////////////////////////// | |
3 | // This macro computes the sampling fraction in the EMCAL | |
4 | // | |
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 | |
8 | // galice2.root . | |
9 | //Begin_Html | |
10 | /* | |
11 | <img src="picts/ComputeSamplingFraction.gif"> | |
12 | */ | |
13 | //End_Html | |
14 | ///////////////////////////////////////////////////////////////////////// | |
15 | if(gAlice){ | |
16 | delete gAlice; | |
17 | gAlice=0; | |
18 | }else{ | |
19 | // Dynamically link some shared libs | |
20 | if(gClassTable->GetID("AliRun") < 0) { | |
21 | gROOT->LoadMacro("loadlibs.C"); | |
22 | loadlibs(); | |
23 | } // end if | |
24 | } // end if gAlice | |
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); | |
28 | ||
29 | // Get AliRun object from file or create it if not on file | |
30 | if(!gAlice) { | |
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"); | |
34 | } // end if !gAlice | |
35 | ||
36 | // Set event pointer to this event | |
37 | Int_t nparticles = gAlice->GetEvent(0); | |
38 | if (nparticles <= 0){ | |
39 | cout << "No particles found for event " << evNumber; | |
40 | cout << " in file " << filename << endl; | |
41 | return; | |
42 | } // end if nparticles <=0 | |
43 | ||
44 | // Pointer to specific detector hits. | |
45 | AliEMCALHit *emcalHit; | |
46 | ||
47 | // Get pointers to ALL Alice detectors and Hits containers | |
48 | AliEMCALv0 *EMCAL = (AliEMCALv0*) gAlice->GetDetector("EMCAL"); | |
49 | if(!EMCAL) { | |
50 | cout << "EMCAL not Found. Exiting." << endl; | |
51 | return; | |
52 | } // end if !EMCAL | |
53 | ||
54 | AliEMCALGeometry *g = EMCAL->GetGeometry(); | |
55 | ||
56 | Int_t i,in,ix,ieta,iphi,ipre; | |
57 | Float_t eta,phi; | |
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; | |
66 | } // | |
67 | return; | |
68 | // Get pointer to the particle | |
69 | TParticle *part; | |
70 | Int_t ipart; | |
71 | ||
72 | // Create histograms | |
73 | TH1F *hEMCAL = new TH1F("hEMCAL" ,"Energy",100,0.,2.); | |
74 | TH2F *hSfs = new TH2F("hSfs","Sampling Fraction",200,-0.7,0.7, | |
75 | 200,0.0,120.0); | |
76 | TH2F *hSfi = new TH2F("hSfi","Sampling Fraction Normilization", | |
77 | 200,-0.7,0.7,200,0.0,120.0); | |
78 | ||
79 | Int_t track,ntracks = gAlice->TreeH()->GetEntries(); | |
80 | Float_t etot = 0.0; | |
81 | Float_t egam = 0.0; | |
82 | Bool_t first = kTRUE; | |
83 | Float_t eta,phi,x,y,z; | |
84 | // Start loop on tracks in the hits containers | |
85 | for(track=0; track<ntracks;track++){ | |
86 | //MI change | |
87 | gAlice->ResetHits(); | |
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())); | |
96 | x = emcalHit->X(); | |
97 | y = emcalHit->Y(); | |
98 | z = emcalHit->Z(); | |
99 | findetaphi(x,y,z,eta,phi); | |
100 | if(first){ | |
101 | egam += part->Energy(); | |
102 | hSfi->Fill(eta,phi,part->Energy()); | |
103 | } // end if first. | |
104 | first = kFALSE; | |
105 | etot += emcalHit->GetEnergy(); | |
106 | hSfs->Fill(eta,phi,emcalHit->GetEnergy()); | |
107 | } // end for tpcHit | |
108 | first = kTRUE; | |
109 | } // end for track | |
110 | hSfs->Divide(hSfi); | |
111 | cout << "Energy of gamma =" << egam << "GeV. Energy in EMCAL =" << etot; | |
112 | cout << "GeV. Sampling Fractions is therefore f=" << etot/egam << endl; | |
113 | ||
114 | //Create a canvas, set the view range, show histograms | |
115 | TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700); | |
116 | c0->Divide(1,3); | |
117 | c0->cd(1); | |
118 | hEMCAL->SetFillColor(42); | |
119 | hEMCAL->Draw(); | |
120 | c0->cd(2); | |
121 | hSfi->Draw("lego"); | |
122 | c0->cd(3); | |
123 | hSfs->Draw("lego"); | |
124 | // c0->SaveAs("analHitsEMCAL.eps"); | |
125 | } | |
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. | |
128 | Float_t r; | |
129 | ||
130 | r = x*x+y*y; | |
131 | r = TMath::Sqrt(r); | |
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); | |
137 | return; | |
138 | } |