source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/macros/tools/MakeNGDensity.C @ 117

using namespace TMath;

/******************************************************************************
 *
 *  How to calculate the radiance
 *
 *  Rad = radiance
 *  N   = Number of photons on the entrance disc
 *  A   = Area of the entrance disk
 *  ThetaMax = max incident angle
 *   
 *  RadCos = (N*Power(Csc(ThetaMax),2))/(A*Pi)
 *
 *  RadFlat = N/(2.*A*Pi*(1 - Cos(ThetaMax)))
 * 
 ******************************************************************************/


//______________________________________________________________________________
void MakeNGDensity() {

    TTree* t = (TTree*)gDirectory->Get("etree");

    Double_t xmin = 0;
    Double_t xmax = 1130; //size of the ideal focal surface
    const Int_t nbins = 100;
    Double_t xbins[nbins+1];
    Double_t area = Pi()*xmax*xmax/nbins;
    Double_t radius = 0;


    Printf("Area:  %f",area);
    for(Int_t i(0); i<nbins; i++) {
        xbins[i] = radius;
        xbins[i+1] = Sqrt((area+Pi()*radius*radius)/Pi()); 
        radius = xbins[i+1];
    }
    
    TH1D* density = new TH1D("density","density",nbins,0.,xmax);
//    TH1D* density = new TH1D("density","density",nbins,xbins);


    cout << "Entries  "<<etree->Draw("sqrt(fPhotons.fIdealFocalPosX*fPhotons.fIdealFocalPosX+fPhotons.fIdealFocalPosY*fPhotons.fIdealFocalPosY) >> density","fPhotons.CrossedIFS()","") << endl;

    Stat_t x;
    TAxis* ax = density->GetXaxis();

    for( Int_t i(1); i <= density->GetNbinsX(); i++){
       x = density->GetBinContent(i);
       
       area = TMath::Pi()*(ax->GetBinUpEdge(i)*ax->GetBinUpEdge(i)
                           -ax->GetBinLowEdge(i)*ax->GetBinLowEdge(i));

       density->SetBinContent(i,x/area);
        
    }
    
    density->Draw();

    
}


//______________________________________________________________________________
void Dummy() {

    EPhoton *photon;
    TClonesArray *array = new TClonesArray();

    Double_t thetamax = 35*DegToRad();
    // parent initialization with useless values
    ParentPhoton* parent = new ParentPhoton(0, spos, pos, 0., 0., 0., 0, kFluorescence);

    Photon* ph;
    Photon* phCopy;

    TRandom* rndm = EsafRandom::Get();

    Int_t ntrans(0), nhits(0);
    TVector3 dir;
    Double_t alpha, radius;

    for( Int_t i(0); i<nph; i++) {
        if ( (i+1)%(nph/100) == 0 )
            cout << "\rTracing " << (i+1)/(nph/100) << '%' << flush;
            
        // new photon position
        radius = sqrt(rndm->Rndm())*gDiscRadius;
        alpha = rndm->Rndm()*TwoPi();
        pos[0] = radius*Cos(alpha);
        pos[1] = radius*Sin(alpha);
        
        // new incident direction
        dir = CosThetaDist(0,thetamax);
        
        ph = &(parent->GetPhoton());
        phCopy = ph;

        // update the parent
        parent->Build(0, spos, pos, dir.Theta(), dir.Phi(), spec->GetLambda(), 0, kFluorescence);

        ph = os->Transport(ph);

        // reject dead photons and going downward
        if (!ph || ph->dir.Z() < 0.) continue;
        ntrans++;
        
        ifs->HitPosition(ph);
        if (!ph->hitIfs) continue;
        
        points[nhits] = ph->posOnIfs;
        nhits++;
        
    }
    cout << "\rTracing done" << endl;

}

//______________________________________________________________________________
TVector3 FlatDist( Double_t theta1, Double_t theta2) {
    //
    //
    //
    
    TVector3 dir;
    Double_t c1 = Cos(theta1);
    Double_t c2 = Cos(theta2);
    dir.SetMagThetaPhi(1.,ACos(rndm->Uniform(c1,c2)), TwoPi()*rndm->Rndm() );

    return dir;
}

//______________________________________________________________________________
TVector3 CosThetaDist( Double_t theta1, Double_t theta2 ) {
    //
    //
    //

    TVector3 dir;
    Double_t c1 = Cos(2*theta1);
    Double_t c2 = Cos(2*theta2);
    dir.SetMagThetaPhi(1.,ACos(rndm->Uniform(c1,c2))/2., TwoPi()*rndm->Rndm() );

    return dir;
}



//______________________________________________________________________________
void PrintDataFile() {
    //
    //
    //
    Printf("%.4f   %e",density->GetXaxis()->GetBinLowEdge(1), density->GetBinContent(1)); 
    for(Int_t i(1); i<=density->GetNbinsX(); i++) 
        Printf("%.4f   %e",density->GetXaxis()->GetBinUpEdge(i), density->GetBinContent(i)); 

}