Small changes
[u/mrichter/AliRoot.git] / EMCAL / testEMCALGeom.C
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 }