source: JEM-EUSO/esaf_cc_at_lal/macros/G4/Efficiency.C @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 4.2 KB
Line 
1
2int angles[3]={323,357,391};
3
4TList l_graphs;
5TList th_graphs;
6
7void MakeAll(const char* fname){
8        l_graphs.SetOwner();
9        th_graphs.SetOwner();
10        l_graphs.Clear();
11        th_graphs.Clear();
12
13    for ( int i(0); i<3; i++ ){
14                Efficiency(fname,0,32,angles[i],50000);
15    }
16
17    TFile f(Form("%s_eff.root",fname),"RECREATE");
18    l_graphs.Write("l_graphs",TObject::kSingleKey);
19    th_graphs.Write("th_graphs",TObject::kSingleKey);
20    f.Close();
21
22}
23
24void Efficiency(const char* fname,int theta1,int theta2,int l,int nphot=50000){
25    int th=theta1;
26    double x_max,y_max;
27
28    while(th<theta2){
29        TFile* f=new TFile(Form("%sa%d.l%d.n%i.root",fname,th,l,nphot),"READ");
30        if (  !f || f->IsZombie() ) {
31            cerr<<"Error opening file."<<endl;
32            return;
33        }
34
35        TTree* tr=(TTree*)f->Get("etree");
36        if ( !tr ) {
37            cerr<<"Error getting tree."<<endl;
38            return;
39        }
40
41        double xn=600.;
42        if(th>10) xn=1200.;
43        TH2F* h1=new TH2F(Form("eff%d",th),Form("eff%d",th),1000,-200.,200.,1000,-200.,xn);
44        tr->Draw(Form("detector.fPhotons.fIdealFocalPosX:detector.fPhotons.fIdealFocalPosY>>eff%d",th),"detector.fPhotons.fStatusBits&0x1","goff");
45
46        int bmax_x, bmax_y, bmax_z;
47        h1->GetMaximumBin(bmax_x, bmax_y, bmax_z);
48
49        x_max = h1->GetXaxis()->GetBinCenter(bmax_x);
50        y_max = h1->GetYaxis()->GetBinCenter(bmax_y);
51
52        h1->GetXaxis()->SetRangeUser(x_max-7.,x_max+7.);
53        h1->GetYaxis()->SetRangeUser(y_max-7.,y_max+7.);
54
55//        TCanvas *c = new TCanvas();
56//        c->Print(Form("%s_Analysis_%d_%d.eps[",fname,th,l));
57//        h1->SetMarkerStyle(8);
58//        h1->SetMarkerSize(0.6);
59//        h1->Draw();
60//        c->Print(Form("%s_Analysis_%d_%d.eps",fname,th,l));
61
62        double box=0.;
63        while( box<=5. ){
64            box+=2.5;
65            double max=0.;
66            float x[5]={max-(box*0.5),max-(box*0.5),max+(box*0.5),max+(box*0.5),max-(box*0.5)};
67            float y[5]={y_max-(box*0.5),y_max+(box*0.5),y_max+(box*0.5),y_max-(box*0.5),y_max-(box*0.5)};
68            TCutG *cut = new TCutG(Form("cut_%fmm",box),5,x,y);
69            cut->SetLineWidth(1);
70            cut->SetLineColor(2);
71
72            double e_total = cut->Integral(h1);
73            double n_fp = h1 ->GetEntries();
74//            cout<<"e_total"<<"      "<<e_total<<endl;
75//            cout<<"e_total"<<"      "<<n_fp<<endl;
76//            cut->Draw("same");
77//            c->Print(Form("%s_Analysis_%d_%d.eps",fname,th,l));
78            //cut->Print();
79
80            TGraph* th_fpeff =    GetGraph(&th_graphs, box, Form("#lambda=%i. FP efficiency",l));
81            TGraph* th_eff =      GetGraph(&th_graphs, box, Form("#lambda=%i. Total efficiency",l));
82            TGraph* th_efficacy = GetGraph(&th_graphs, box, Form("#lambda=%i. Efficacy",l));
83            TGraph* l_fpeff =     GetGraph(&l_graphs , box, Form("#theta=%i. FP efficiency",th));
84            TGraph* l_eff =       GetGraph(&l_graphs , box, Form("#theta=%i. Total efficiency",th));
85            TGraph* l_efficacy =  GetGraph(&l_graphs , box, Form("#theta=%i. Efficacy",th));
86
87            double ttotalefficiency=e_total/nphot;
88            double tefficiency=e_total/n_fp;
89            double tefficacy=ttotalefficiency*TMath::Cos(th*TMath::DegToRad())*1.250*1.250*TMath::Pi();
90
91            th_fpeff->SetPoint(th_fpeff->GetN(), th, tefficiency);
92            th_eff->SetPoint(th_eff->GetN(), th, ttotalefficiency);
93            th_efficacy->SetPoint(th_efficacy->GetN(), th, tefficacy );
94            l_fpeff->SetPoint(l_fpeff->GetN(), l, tefficiency);
95            l_eff->SetPoint(l_eff->GetN(), l, ttotalefficiency);
96            l_efficacy->SetPoint(l_efficacy->GetN(), l, tefficacy);
97        }
98
99//        c->Print(Form("%s_Analysis_%d_%d.eps]",fname,th,l));
100        th+=5;
101
102//              delete c;
103                delete h1;
104    }
105}
106
107TGraph* GetGraph(TList* list, Double_t box, const char* title){
108    TString name=Form("%s. Box %.2f X %.2f cm S = %.2fcm^{2}",title,box/10.,box/10.,(box)*(box)/100.);
109    TGraph* gr=(TGraph*)list->FindObject(name.Data());
110    if( gr ) return gr;
111    else {
112        TGraph *gr = new TGraph();
113        gr->SetName(name.Data());
114        gr->SetTitle(name.Data());
115        list->Add(gr);
116    }
117    return gr;
118}
Note: See TracBrowser for help on using the repository browser.