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 }