void testEMCALGeom(const char *filename="galice.root",){ ///////////////////////////////////////////////////////////////////////// // This macro computes the sampling fraction in the EMCAL // // Root > .L ComputeSamplingFraction.C //this loads the macro in memory // Root > ComputeSamplingFraction(); //by default process first event // Root > ComputeSamplingFraction("galice2.root"); //process from the file // galice2.root . //Begin_Html /* */ //End_Html ///////////////////////////////////////////////////////////////////////// if(gAlice){ delete gAlice; gAlice=0; }else{ // Dynamically link some shared libs if(gClassTable->GetID("AliRun") < 0) { gROOT->LoadMacro("loadlibs.C"); loadlibs(); } // end if } // end if gAlice // Connect the Root Galice file containing Geometry, Kine and Hits TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename); if(!file) file = new TFile(filename); // Get AliRun object from file or create it if not on file if(!gAlice) { gAlice = (AliRun*)file->Get("gAlice"); if(gAlice) printf("AliRun object found on file\n"); if(!gAlice) gAlice = new AliRun("gAlice","Alice test program"); } // end if !gAlice // Set event pointer to this event Int_t nparticles = gAlice->GetEvent(0); if (nparticles <= 0){ cout << "No particles found for event " << evNumber; cout << " in file " << filename << endl; return; } // end if nparticles <=0 // Pointer to specific detector hits. AliEMCALHit *emcalHit; // Get pointers to ALL Alice detectors and Hits containers AliEMCALv0 *EMCAL = (AliEMCALv0*) gAlice->GetDetector("EMCAL"); if(!EMCAL) { cout << "EMCAL not Found. Exiting." << endl; return; } // end if !EMCAL AliEMCALGeometry *g = EMCAL->GetGeometry(); Int_t i,in,ix,ieta,iphi,ipre; Float_t eta,phi; cout << "i ieta iphi ipre in eta phi ix" << endl; for(i=0;i<=2*g->GetNEta()*g->GetNPhi();i++){ g->TowerIndexes(i,ieta,iphi,ipre); in = g->TowerIndex(ieta,iphi,ipre); g->EtaPhiFromIndex(i,eta,phi); ix = g->TowerIndexFromEtaPhi(eta,phi); cout << i << " " << ieta << " " << iphi << " " << ipre; cout << " " << in << " " << eta << " " << phi << " " << ix << endl; } // return; // Get pointer to the particle TParticle *part; Int_t ipart; // Create histograms TH1F *hEMCAL = new TH1F("hEMCAL" ,"Energy",100,0.,2.); TH2F *hSfs = new TH2F("hSfs","Sampling Fraction",200,-0.7,0.7, 200,0.0,120.0); TH2F *hSfi = new TH2F("hSfi","Sampling Fraction Normilization", 200,-0.7,0.7,200,0.0,120.0); Int_t track,ntracks = gAlice->TreeH()->GetEntries(); Float_t etot = 0.0; Float_t egam = 0.0; Bool_t first = kTRUE; Float_t eta,phi,x,y,z; // Start loop on tracks in the hits containers for(track=0; trackResetHits(); gAlice->TreeH()->GetEvent(track); for(emcalHit=(AliEMCALHit*)EMCAL->FirstHit(-1);emcalHit; emcalHit=(AliEMCALHit*)EMCAL->NextHit()) { ipart = emcalHit->GetTrack(); part = gAlice->Particle(ipart); if(part->GetFirstMother()!=0) continue; // only those that don't interact. hEMCAL->Fill((Float_t)(emcalHit->GetEnergy())); x = emcalHit->X(); y = emcalHit->Y(); z = emcalHit->Z(); findetaphi(x,y,z,eta,phi); if(first){ egam += part->Energy(); hSfi->Fill(eta,phi,part->Energy()); } // end if first. first = kFALSE; etot += emcalHit->GetEnergy(); hSfs->Fill(eta,phi,emcalHit->GetEnergy()); } // end for tpcHit first = kTRUE; } // end for track hSfs->Divide(hSfi); cout << "Energy of gamma =" << egam << "GeV. Energy in EMCAL =" << etot; cout << "GeV. Sampling Fractions is therefore f=" << etot/egam << endl; //Create a canvas, set the view range, show histograms TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700); c0->Divide(1,3); c0->cd(1); hEMCAL->SetFillColor(42); hEMCAL->Draw(); c0->cd(2); hSfi->Draw("lego"); c0->cd(3); hSfs->Draw("lego"); // c0->SaveAs("analHitsEMCAL.eps"); } void findetaphi(Float_t x,Float_t y,Float_t z,Float_t &eta,Float_t &phi){ // Compute Eta and Phi from x, y, z. Assumes vertex is a zero. Float_t r; r = x*x+y*y; r = TMath::Sqrt(r); phi = 180.*TMath::ATan2(y,x)/TMath::Pi(); if(phi<0.0) phi += 180.; r = TMath::ATan2(r,z); r = TMath::Tan(0.5*r); eta = -TMath::Exp(r); return; }