source: JEM-EUSO/esaf_cc_at_lal/macros/ETreeEventViewer.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: 100.9 KB
Line 
1// Author: A.Thea, D.Naumov   Aug,  9 2004
2
3///*****************************************************************************
4//
5// EventViewer macro:
6//
7// ETreeEventViewer is an interface for the classes of eventviewer library made to
8// display multiple events: 3D shower viewer, some histos related to the shower.
9// With time we will implement also detector info for many events
10//
11// See also: EEventViewer macro.
12//
13//*****************************************************************************/
14
15Bool_t axis = kTRUE;
16
17
18ETreePainter* treepainter = NULL;
19TTree *etree = NULL;
20TTree *stree = NULL;
21TChain *chain = NULL;
22Bool_t SAVED = false;  // to know if stree has been already saved
23EGViewer* pad_shower = 0;
24EGViewer* pad_atmo_photons = 0;
25EGViewer* pad_p_on_pupil = 0;
26
27const char *defaultfilename= "";
28const char *defaultStree= "";
29
30
31
32// !!!!!!!!!!!!!!!!! WARNING : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33//
34//      - GammaA-B   :  - should modify ETreePainter (X1 pb) when building Stree
35//
36//      - other than Zeta_A, Zeta_Mono, Gamma_A : change _air_ into _back_   (scattered ckov)
37//
38// !!!!!!!!!!!!!!!!! WARNING : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
39
40
41
42TCut cut1 = "H_max_last_bin && X_max_last_bin &&L_max_last_bin && H_fluomax_last_bin &&X_fluomax_last_bin && L_fluomax_last_bin";
43TCut cut2 = "X1>0.002";  // because some pb below (shower generation)
44TCut cut3 = "";
45TCut cut4OR = "";
46TCut cut5OR = "";
47TCut drawcut;
48
49
50
51//__________________________________________________________________________________________________________________________
52void ETreeEventViewer() {
53   
54   bar = new TControlBar("vertical", "ETree Event Viewer",10,10);
55   
56   bar->AddButton("How to run        ",                 "help()","Instructions for running this macro"); 
57   bar->AddButton("Open File         ",                 "openfile()","Open a new file"); 
58   bar->AddButton("Open Stree        ",                 "openstree()","Open a stree file"); 
59   bar->AddButton("Browser           ",                 "TBrowser b", "Open a Browser"); 
60   bar->AddButton("STREE Viewer",                       "simpletreeviewer()","Simple Results Tree Viewer"); 
61   bar->AddButton("Cut summary",                        "CutSummary()",""); 
62   bar->AddButton("Shower - Truth Analysis",            "ShowerAnalysis()","Shower analysis");
63   bar->AddButton("Photons at creation Analysis",       "AtmoPhotonsAnalysis()","Analysis of Photons created by shower");
64   bar->AddButton("Photons on pupil analysis",          "PhotonsOnPupilAnalysis()","Show Tree Atmo Histos");
65   bar->AddButton("Quit Root        ",                  ".q","Quit root");
66   bar->Show();
67   
68   gROOT->SaveContext();
69   Bool_t def = OpenDefaultStree();
70   if(!def) openstree();
71   
72   drawcut += cut1;
73   drawcut += cut2;
74   drawcut += cut3;
75   drawcut = drawcut || cut4OR;
76   drawcut = drawcut || cut5OR;
77}
78
79//__________________________________________________________________________________________________________________________
80void ShowerAnalysis() {
81   bar = new TControlBar("vertical", "Shower and truth analysis",10,10);
82   
83   bar->AddButton("Truth Histograms",                   "treetruthhistos()","Show Tree  Truth Histos");
84   bar->AddButton("Shower energy - xmax correlations",  "treeshowercorrel()","Show Tree shower Correlations");
85   bar->AddButton("Shower theta correlations",          "shower_theta_correl()","Show Tree shower Correlations");
86   bar->Show();
87   
88   gROOT->SaveContext();
89   if(!pad_shower) pad_shower = new EGViewer();
90   else pad_shower->Clear();
91   pad_shower->SetName("Shower Analysis");
92   pad_shower->SetTitle("Shower Analysis");
93}
94
95//__________________________________________________________________________________________________________________________-
96void AtmoPhotonsAnalysis() {
97   bar = new TControlBar("vertical", "Photons at creation analysis",10,10);
98   
99   bar->AddButton("Electrons / Photons correlations",   "phot_elec_correl()","");
100   bar->AddButton("Crash ground",                       "crash_ground()","");
101   bar->AddButton("Xmax shift",                         "Xmax_shift()","");
102   bar->AddButton("photons Histo",                      "treeatmohistos()","");
103   bar->AddButton("photons / theta correlations",       "atmo_theta_correl()","");
104   bar->AddButton("photons / energy correlations",      "atmo_energy_correl()","");
105   bar->Show();
106
107   gROOT->SaveContext();
108   if(!pad_atmo_photons) pad_atmo_photons = new EGViewer();
109   else pad_atmo_photons->Clear();
110   pad_atmo_photons->SetName("Photons in atmosphere Analysis");
111   pad_atmo_photons->SetTitle("Photons in atmosphere Analysis");
112}
113
114//__________________________________________________________________________________________________________________________-
115void PhotonsOnPupilAnalysis() {
116   bar = new TControlBar("vertical", "Photons on pupil analysis",10,10);
117   
118   bar->AddButton("Photons on pupil Histograms",        "treepupilhistos()","");
119   bar->AddButton("Pupil / theta Correlations",         "pupil_theta_correl()","");
120   bar->AddButton("Pupil / energy Correlations",        "pupil_energy_correl()","");
121   bar->AddButton("Transmissions",                      "transmission()","");
122   bar->AddButton("Miscellaneous",                      "Miscellaneous()","");
123   bar->Show();
124
125   gROOT->SaveContext();
126   if(!pad_p_on_pupil) pad_p_on_pupil = new EGViewer();
127   else pad_p_on_pupil->Clear();
128   pad_p_on_pupil->SetName("Photons on pupil Analysis");
129   pad_p_on_pupil->SetTitle("Photons on pupil Analysis");
130}
131
132//__________________________________________________________________________________________________________________________-
133void help() {
134   //
135 
136   new TCanvas("chelp","Help to run demos",200,10,700,600);
137
138   welcome = new TPaveText(.1,.8,.9,.97);
139   welcome->AddText("Welcome to ESAF Tree Event Viewer");
140   welcome->SetTextFont(32);
141   welcome->SetTextColor(4);
142   welcome->SetFillColor(24);
143   welcome->Draw();
144
145   hdemo = new TPaveText(.05,.05,.95,.7);
146   hdemo->SetTextAlign(12);
147   hdemo->SetTextFont(52);
148   hdemo->AddText("- Demo for visualization of the shower development in atmosphere");
149   hdemo->AddText("  and signal observed by the EUSO detector");
150   hdemo->AddText("- Click left mouse button to execute one demo");
151   hdemo->AddText("- Click left mouse button on Animation On/Off to start/stop animation");
152   hdemo->AddText(" ");
153   hdemo->SetAllWith("....","color",2);
154   hdemo->SetAllWith("....","font",72);
155   hdemo->SetAllWith("....","size",0.03);
156
157   hdemo->Draw();
158}
159
160//__________________________________________________________________________________________________________________________-
161void CutSummary() {
162    //
163    if(!stree) Printf("NO STREE : open stree before");
164   
165    TCanvas *c = 0;
166    if (( c = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c_all") )) {
167        c->cd(kTruthPadID);
168        gPad->Clear();
169    } 
170    else {
171        c = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c_truth");
172        if (!c) c = new TCanvas("c_truth","c_truth");
173
174        c->Clear();
175        c->cd();
176    }
177    TPaveText* hdemo = new TPaveText(.05,.05,.95,.95);
178    hdemo->SetTextAlign(12);
179    hdemo->SetTextFont(52);
180    hdemo->SetTextColor(4);
181    TText* temptext = 0;
182   
183    hdemo->AddText("");
184    temptext = hdemo->AddText("Summarize cut effects on nb of events:");
185    temptext->SetTextColor(1);
186   
187    hdemo->AddText("");
188    temptext = hdemo->AddText(Form("Cut1 =  %s ",cut1.GetTitle()));
189    temptext->SetTextColor(4);
190    stree->Draw("Energy>>httemp",cut1,"goff");
191    TH1F* hhtemp = (TH1F*)gROOT->FindObject("httemp");
192    temptext = hdemo->AddText(Form("        Nb of events                  %g",hhtemp->GetEntries()));
193    temptext->SetTextColor(2);
194   
195    hdemo->AddText("");
196    temptext = hdemo->AddText(Form("Cut2 =  %s ",cut2.GetTitle()));
197    temptext->SetTextColor(4);
198    stree->Draw("Energy>>httemp",cut2,"goff");
199    hhtemp = (TH1F*)gROOT->FindObject("httemp");
200    temptext = hdemo->AddText(Form("        Nb of events      %g",hhtemp->GetEntries()));
201    temptext->SetTextColor(2);
202    stree->Draw("Energy>>httemp",cut1+cut2,"goff");
203    hhtemp = (TH1F*)gROOT->FindObject("httemp");
204    temptext = hdemo->AddText(Form("        Cumulated nb      %g",hhtemp->GetEntries()));
205    temptext->SetTextColor(2);
206       
207    hdemo->AddText("");
208    temptext = hdemo->AddText(Form("Cut3 =  %s ",cut3.GetTitle()));
209    temptext->SetTextColor(4);
210    stree->Draw("Energy>>httemp",cut3,"goff");
211    hhtemp = (TH1F*)gROOT->FindObject("httemp");
212    temptext = hdemo->AddText(Form("        Nb of events      %g",hhtemp->GetEntries()));
213    temptext->SetTextColor(2);
214    stree->Draw("Energy>>httemp",cut1+cut2+cut3,"goff");
215    hhtemp = (TH1F*)gROOT->FindObject("httemp");
216    temptext = hdemo->AddText(Form("        Cumulated nb      %g",hhtemp->GetEntries()));
217    temptext->SetTextColor(2);
218       
219    hdemo->AddText("");
220    temptext = hdemo->AddText(Form("Cut4OR =  %s ",cut4OR.GetTitle()));
221    temptext->SetTextColor(4);
222    stree->Draw("Energy>>httemp",cut4OR,"goff");
223    hhtemp = (TH1F*)gROOT->FindObject("httemp");
224    temptext = hdemo->AddText(Form("        Nb of events      %g",hhtemp->GetEntries()));
225    temptext->SetTextColor(2);
226    stree->Draw("Energy>>httemp",cut1+cut2+cut3||cut4OR,"goff");
227    hhtemp = (TH1F*)gROOT->FindObject("httemp");
228    temptext = hdemo->AddText(Form("        Cumulated nb      %g",hhtemp->GetEntries()));
229    temptext->SetTextColor(2);
230       
231    hdemo->AddText("");
232    temptext = hdemo->AddText(Form("Cut5OR =  %s ",cut5OR.GetTitle()));
233    temptext->SetTextColor(4);
234    stree->Draw("Energy>>httemp",cut5OR,"goff");
235    hhtemp = (TH1F*)gROOT->FindObject("httemp");
236    temptext = hdemo->AddText(Form("        Nb of events      %g",hhtemp->GetEntries()));
237    temptext->SetTextColor(2);
238    stree->Draw("Energy>>httemp",cut1+cut2+cut3||cut4OR||cut5OR,"goff");
239    hhtemp = (TH1F*)gROOT->FindObject("httemp");
240    temptext = hdemo->AddText(Form("        Cumulated nb      %g",hhtemp->GetEntries()));
241    temptext->SetTextColor(2);
242       
243    hdemo->AddText("");
244    hdemo->Draw();
245
246}
247
248//__________________________________________________________________________________________________________________________-
249void gettree() {
250    etree = (TTree*)gDirectory->Get("etree");
251    if( etree ) {
252        Printf("etree with %d events found.", (Int_t)etree->GetEntries());
253        if(treepainter) delete treepainter;
254        treepainter = new ETreePainter( etree );
255    }
256   
257    stree = (TTree*)gDirectory->Get("stree"); 
258    if(!stree) {
259      Warning("gettree()","Fail to generate stree");
260      return;
261    }
262    savesimpletree();
263}
264
265//__________________________________________________________________________________________________________________________-
266void getstree() {
267    stree = (TTree*)gDirectory->Get("stree");
268    if( stree ) {
269        Printf("%d events found.", (Int_t)stree->GetEntries());
270        if(treepainter) delete treepainter;
271        treepainter = new ETreePainter( stree );
272        SAVED = true;
273    }
274    else Warning("getstree()","stree not found");
275}
276
277//__________________________________________________________________________________________________________________________-
278Bool_t checktree( Option_t *opt = "default" ) {
279    TString option(opt);
280    Bool_t rtn = true;
281    if((opt == "etree") && !etree ) {
282        Warning("I do not find etree");
283        rtn = false;
284    }
285    else if((opt == "stree") && !stree ) {
286        Warning("I do not find stree");
287        rtn = false;
288        //help();
289   
290    }
291    else if((opt == "chain") && !chain ) {
292        Warning("I do not find chain");
293        rtn = false;
294    }
295    else if((opt == "default")&&!etree&&!chain&&!stree) {
296        Warning("checktree()","NO TREE FOUND (no chain, etree or stree)");
297        rtn = false;
298    }
299   
300    return rtn;
301}
302
303//__________________________________________________________________________________________________________________________-
304void getchain() {
305  //    chain = (TChain*)gDirectory->Get("etree");
306    if(chain) {
307        Printf("chain etree with %d events found.\n", (Int_t)chain->GetEntries());
308        if(treepainter)
309           delete treepainter;
310        treepainter = new ETreePainter( chain );
311    }
312   
313    stree = (TTree*)gDirectory->Get("stree"); 
314    if(!stree) {
315      Warning("getchain()","Fail to generate stree");
316      return;
317    }
318    savesimpletree();
319}
320
321//__________________________________________________________________________________________________________________________-
322void openfile() {
323    TFile *f = 0;
324    delete f;
325    TGFileInfo* info = new TGFileInfo();
326    while ( info->fFilename == 0 || !f || f->IsZombie() ) {
327        TGFileDialog* filedialog = new TGFileDialog(gClient->GetRoot(), 0, kFDOpen, info);
328        if( info->fFilename ) 
329        f = new TFile( info->fFilename );
330    }
331
332    if( f && !f->IsZombie() ) {
333        Printf("%s successfully opened.", info->fFilename);
334        gettree();
335        }
336    else 
337        Printf("Unable to open %s",info->fFilename);   
338}
339
340//__________________________________________________________________________________________________________________________-
341void openstree() {
342    TFile *f = 0;
343    delete f;
344    TGFileInfo* info = new TGFileInfo();
345    while ( info->fFilename == 0 || !f || f->IsZombie() ) {
346        TGFileDialog* filedialog = new TGFileDialog(gClient->GetRoot(), 0, kFDOpen, info);
347        if( info->fFilename ) 
348        f = new TFile( info->fFilename );
349    }
350
351    if( f && !f->IsZombie() ) {
352        Printf("%s successfully opened.", info->fFilename);
353        Warning("openstree()","SHOWER DRAWINGS NOT USABLE in this configuration");
354        getstree();
355    }
356    else 
357        Printf("Unable to open %s",info->fFilename);   
358}
359
360//__________________________________________________________________________________________________________________________-
361Bool_t OpenDefaultStree() {
362    TFile *f = new TFile(defaultStree);
363
364    if( f && !f->IsZombie() ) {
365        Printf("Default Stree %s successfully opened. ", defaultStree);
366        getstree();
367        return true;
368    }
369    else {
370        Printf("Unable to open default Stree %s",defaultStree);
371        return false;
372    }
373}
374
375//__________________________________________________________________________________________________________________________-
376void treetruthhistos() {
377    Bool_t flag = true;
378    if(!chain&&!etree&&!stree) flag = checktree();     
379    if(!flag) {
380        Warning("treetruthhistos()","No Tree found");
381        return;
382    }
383   
384    Double_t min = 0;
385    Double_t max = 0;
386    Double_t textpos = 0;
387    Double_t par[3];
388    TLatex *t = 0;
389
390   
391    if(!pad_shower) pad_shower = new EGViewer();
392    TH1* histo = 0;
393   
394    TCanvas *c = pad_shower->FindCanvas("c_treetruth_histo");
395    if(!c ) c = pad_shower->AddTab("c_treetruth_histo","c_treetruth_histo");
396    TCanvas *c2 = pad_shower->FindCanvas("c_treetruth_histo2");
397    if(!c2 ) c2 = pad_shower->AddTab("c_treetruth_histo2","c_treetruth_histo2");
398
399    c->Clear();
400    c->cd();
401
402    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
403    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
404    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
405    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
406
407    pad1->Draw();
408    pad2->Draw();
409    pad3->Draw();
410    pad4->Draw();
411
412    pad1->cd();
413    treepainter->DrawHistos("Energy",drawcut);
414    histo = (TH1*)gPad->GetPrimitive("h_Energy");
415    histo->SetTitle("Energy distribution in log scale");
416    histo->SetXTitle("log10(Energy) (MeV)");
417
418    pad2->cd();
419    pad2->SetLogy(1);
420    treepainter->DrawHistos("X1",drawcut);
421    histo = (TH1*)gPad->GetPrimitive("h_X1");
422    histo->SetTitle("Depth of first interaction point");
423    histo->SetXTitle("X1 (g/cm^{2})");
424    min = 0;
425    max = 0;
426    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
427        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
428        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
429        if(fabs(min*max) > 0) break;
430    }
431    mypol = new TF1("mypol_X1","expo");
432    mypol->SetLineWidth(1);
433    mypol->SetLineColor(2);
434    textpos = histo->GetBinLowEdge(1) + 0.05*(histo->GetBinLowEdge(histo->GetXaxis()->GetLast()) - histo->GetBinLowEdge(1));
435    histo->Fit("mypol_X1","WOQ");
436    mypol->GetParameters(par);
437    t = new TLatex();
438    char txt[200];
439    t->SetTextFont(62); 
440    t->SetTextSize(0.06);
441    t->SetTextColor(2);
442    t->SetTextColor(1);
443    sprintf(txt,"%.2g  exp(%.2g X_{1})",exp(par[0]),par[1]);
444    t->DrawLatex(textpos,pow(10,log10(histo->GetMaximum())*0.15),txt);
445
446    pad3->cd();
447    treepainter->DrawHistos("Phi",drawcut);
448    histo = (TH1*)gPad->GetPrimitive("h_Phi");
449    histo->SetTitle("Azimuth angle");
450    histo->SetXTitle("phi (deg)");
451
452    pad4->cd();
453    treepainter->DrawHistos("Theta",drawcut);
454    histo = (TH1*)gPad->GetPrimitive("h_Theta");
455    histo->SetTitle("Local #theta");
456    histo->SetXTitle("#theta (deg)");
457    min = 0;
458    max = 0;
459    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
460        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
461        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
462        if(fabs(min*max) > 0) break;
463    }
464    mypol = new TF1("mypol_Theta","[0]*sin(2*(x*3.14/180))");
465    mypol->SetLineWidth(1);
466    mypol->SetLineColor(2);
467    textpos = min + 0.3*(max-min);
468    histo->Fit("mypol_Theta","WOQ");
469    mypol->GetParameters(par);
470    t = new TLatex();
471    char txt[200];
472    t->SetTextFont(62); 
473    t->SetTextSize(0.06);
474    t->SetTextColor(2);
475    t->SetTextColor(1);
476    sprintf(txt,"%.2f sin(2x) fit",par[0]);
477    t->DrawText(textpos,histo->GetMaximum()*0.84,txt);
478   
479
480/////////////////////////////
481    c2->Clear();
482    c2->cd();
483
484    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
485    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
486    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
487    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
488
489    pad11->Draw();
490    pad21->Draw();
491    pad31->Draw();
492    pad41->Draw();
493
494    pad11->cd();
495    treepainter->DrawHistos("Hmax",drawcut);
496    histo = (TH1*)gPad->GetPrimitive("h_Hmax");
497    histo->SetTitle("Height of maximum");
498    histo->SetXTitle("Hmax (km)");
499   
500    pad21->cd();
501    pad21->SetLogy(1);
502    treepainter->DrawHistos("Xmax",drawcut);
503    histo = (TH1*)gPad->GetPrimitive("h_Xmax");
504    histo->SetTitle("Depth of maximum");
505    histo->SetXTitle("Xmax (g/cm^{2})");
506    min = histo->GetBinCenter(histo->GetMaximumBin()) + 10;
507    max = histo->GetBinLowEdge(histo->GetXaxis()->GetLast()) + histo->GetBinWidth(histo->GetXaxis()->GetLast());
508    mypol = new TF1("mypol_Xmax","expo",min,max);   // fit the Xmax distribution tail (10 g/cm2 after maximum)
509    mypol->SetLineWidth(1);
510    mypol->SetLineColor(2);
511    textpos = histo->GetBinLowEdge(1) + 0.05*(histo->GetBinLowEdge(histo->GetXaxis()->GetLast()) - histo->GetBinLowEdge(1));
512    histo->Fit("mypol_Xmax","WOQR");
513    mypol->GetParameters(par);
514    t = new TLatex();
515    char txt[200];
516    t->SetTextFont(62); 
517    t->SetTextSize(0.06);
518    t->SetTextColor(2);
519    t->SetTextColor(1);
520    sprintf(txt,"%.2g  exp(%.2g X_{max})",exp(par[0]),par[1]);
521    t->DrawLatex(textpos,pow(10,log10(histo->GetMaximum())*0.15),txt);
522   
523    pad31->cd();
524    treepainter->DrawHistos("Xmax_X1",drawcut);
525    pad31->SetLogy(1);
526    histo = (TH1*)gPad->GetPrimitive("h_Xmax_X1");
527    histo->SetTitle("Depth of maximum from depth of first interaction");
528    histo->SetXTitle("Xmax - X1 (g/cm^{2})");
529   
530    pad41->cd();
531    pad41->SetLogy(1);
532    treepainter->DrawHistos("Nemax",drawcut);
533    histo = (TH1*)gPad->GetPrimitive("h_Nemax");
534    histo->SetTitle("Number of electrons at maximum");
535    histo->SetXTitle("Nemax");
536}
537
538//__________________________________________________________________________________________________________________________-
539void treeshowercorrel() {
540   
541    Bool_t flag = true;
542    if(!chain&&!etree&&!stree) flag = checktree();     
543    if(!flag) {
544        Warning("treeshowercorrel()","No Tree found");
545        return;
546    }
547    Double_t min = 0;
548    Double_t max = 0;
549    Double_t textpos = 0;
550    Double_t par[3];
551    TLatex *t = 0;
552    Int_t nbin = 0;
553    TF1* fitfunc = 0;
554    if(!pad_shower) pad_shower = new EGViewer();
555
556    TCanvas *c = pad_shower->FindCanvas("c_treeshower_correl");
557    if(!c ) c = pad_shower->AddTab("c_treeshower_correl","c_treeshower_correl");
558   
559    TCanvas *c2 = pad_shower->FindCanvas("c_treeshower_correl_F");
560    if(!c2 ) c2 = pad_shower->AddTab("c_treeshower_correl_F","c_treeshower_correl_F");
561
562
563    // *************** PLOTS FOR BEHAVIOR ********************
564   
565    c->Clear();
566    c->cd();
567
568    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
569    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
570    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
571    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
572
573    pad1->Draw();
574    pad2->Draw();
575    pad3->Draw();
576    pad4->Draw();
577   
578    pad1->cd();
579    treepainter->DrawHistos("Hmax_Theta",drawcut);
580    histo = (TH1*)gPad->GetPrimitive("h_Hmax_Theta");
581    histo->SetTitle("Shower Height of maximum vs/ local #theta");
582    histo->SetXTitle("#theta (deg)");
583    histo->SetYTitle("Hmax (km)");
584    histo->SetStats(0);
585    min = 0;
586    max = 0;
587    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
588        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
589        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
590        if(fabs(min*max) > 0) break;
591    }
592    /*
593    mypol = new TF1("mypol_Hmax_Theta",my_logcos,min,max,2);
594    mypol->SetLineWidth(1);
595    mypol->SetLineColor(2);
596    mypol->SetParameter(1.,-10.);
597    mypol->SetParLimits(0,0.01,10.);
598    mypol->SetParLimits(1,-10.,-0.1);
599    textpos = min + 0.1*(max-min);
600    histo->Fit("mypol_Hmax_Theta","WROQ");
601    mypol->GetParameters(par);
602    t = new TLatex();
603    char txt[200];
604    t->SetTextFont(62); 
605    t->SetTextSize(0.06);
606    t->SetTextColor(2);
607    t->SetTextColor(1);
608    sprintf(txt,"%.3g  -  %.3g * log(cos#theta)",par[0],-par[1]);
609    t->DrawLatex(textpos,histo->GetMaximum()*0.64,txt);
610    */
611   
612    pad2->cd();
613    treepainter->DrawHistos("Xmax_Energy",drawcut);
614    profXmax = (TProfile*)gPad->GetPrimitive("h_Xmax_Energy");
615    profXmax->SetTitle("Shower Depth of maximum vs/ Energy");
616    profXmax->SetXTitle("log10(Energy) (MeV)");
617    profXmax->SetYTitle("Xmax (g/cm^{2})");
618    profXmax->SetStats(0);
619    min = 0;
620    max = 0;
621    for(Int_t m=0; m<profXmax->GetNbinsX(); m++) {
622        if(profXmax->GetBinContent(m+1) && fabs(min) == 0) min = profXmax->GetBinLowEdge(m+1);
623        if(profXmax->GetBinContent(profXmax->GetNbinsX() - m) && fabs(max) == 0) max = profXmax->GetBinLowEdge(profXmax->GetNbinsX() - m);
624        if(fabs(min*max) > 0) break;
625    }
626    textpos = min + 0.03*(max-min);
627    /*
628    TF1* mypol = new TF1("mypol_Xmax_Energy","pol1",min,max);
629    mypol->SetLineWidth(1);
630    mypol->SetLineColor(2);
631    profXmax->Fit("mypol_Xmax_Energy","WRCOQ");
632    mypol->GetParameters(par);
633    t = new TLatex();
634    char txt[200];
635    t->SetTextFont(62); 
636    t->SetTextSize(0.06);
637    t->SetTextColor(2);
638    t->SetTextColor(1);
639    sprintf(txt,"Xmax = %.3g  +  %.3g * log_{10}(E/MeV) g/cm^{2}",par[0],par[1]);
640    t->DrawLatex(textpos,profXmax->GetMaximum()*0.54,txt);
641    */
642
643    pad3->cd();
644    treepainter->DrawHistos("Nemax_Energy",drawcut);
645    profNemax = (TProfile*)gPad->GetPrimitive("h_Nemax_Energy");
646    profNemax->SetTitle("Number of electrons at maximum vs/ Energy");
647    profNemax->SetXTitle("log10(Energy) (MeV)");
648    profNemax->SetYTitle("log10(Nemax)");
649    profNemax->SetStats(0);
650    min = 0;
651    max = 0;
652    for(Int_t m=0; m<profNemax->GetNbinsX(); m++) {
653        if(profNemax->GetBinContent(m+1) && fabs(min) == 0) min = profNemax->GetBinLowEdge(m+1);
654        if(profNemax->GetBinContent(profNemax->GetNbinsX() - m) && fabs(max) == 0) max = profNemax->GetBinLowEdge(profNemax->GetNbinsX() - m);
655        if(fabs(min*max) > 0) break;
656    }
657    mypol = new TF1("mypol_Nemax_Energy",my_linear,min,max,2);  // pente fixee a 1
658    mypol->SetLineWidth(1);
659    mypol->SetLineColor(2);
660    mypol->SetParameter(1.,1.);
661    mypol->SetParLimits(0,-10.,10);
662    mypol->SetParLimits(1,1.,1.);
663    textpos = min + 0.1*(max-min);
664    profNemax->Fit("mypol_Nemax_Energy","WRCOQ");
665    mypol->GetParameters(par);
666    t = new TLatex();
667    char txt[200];
668    t->SetTextFont(62); 
669    t->SetTextSize(0.06);
670    t->SetTextColor(2);
671    t->SetTextColor(1);
672    sprintf(txt,"Ne_{max} = %.3e  *  (E/MeV)",pow(10,par[0]));
673    t->DrawLatex(textpos,profNemax->GetMaximum()*0.54,txt);
674
675    pad4->cd();
676    treepainter->DrawHistos("Ne2_Energy",drawcut);
677    profNe2 = (TProfile*)gPad->GetPrimitive("h_Ne2_Energy");
678    profNe2->SetTitle("Ne Longitudinal profile first HALF vs/ Energy");
679    profNe2->SetXTitle("log10(Energy) (MeV)");
680    profNe2->SetYTitle("log10(Ne Longitudinal profile first HALF) (g/cm^{2})");
681    profNe2->SetStats(0);
682    min = 0;
683    max = 0;
684    for(Int_t m=0; m<profNe2->GetNbinsX(); m++) {
685        if(profNe2->GetBinContent(m+1) && fabs(min) == 0) min = profNe2->GetBinLowEdge(m+1);
686        if(profNe2->GetBinContent(profNe2->GetNbinsX() - m) && fabs(max) == 0) max = profNe2->GetBinLowEdge(profNe2->GetNbinsX() - m);
687        if(fabs(min*max) > 0) break;
688    }
689    mypol = new TF1("mypol_Ne2_Energy",my_linear,min,max,2);  // pente fixee a 1
690    mypol->SetLineWidth(1);
691    mypol->SetLineColor(2);
692    mypol->SetParameter(1.,1.);
693    mypol->SetParLimits(0,-10.,10);
694    mypol->SetParLimits(1,1.,1.);
695
696    textpos = min + 0.1*(max-min);
697    profNe2->Fit("mypol_Ne2_Energy","WRCOQ");
698    mypol->GetParameters(par);
699    t = new TLatex();
700    char txt[200];
701    t->SetTextFont(62); 
702    t->SetTextSize(0.06);
703    t->SetTextColor(2);
704    t->SetTextColor(1);
705    sprintf(txt,"\"Ne_{1/2}\" = %.3e  *  (E/MeV)  g/cm^{2}",pow(10,par[0]));
706    t->DrawLatex(textpos,profNe2->GetMaximum()*0.54,txt);
707   
708   
709    // *************** PLOTS FOR SIGMA  ****************************************
710
711    c2->Clear();
712    c2->cd();
713    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
714    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
715    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
716    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
717
718    pad11->Draw();
719    pad21->Draw();
720    pad31->Draw();
721    pad41->Draw();
722    TProfile* prof = 0;
723   
724    pad11->cd();
725    prof = profXmax;
726    nbin = profXmax->GetNbinsX();
727    min = profXmax->GetBinLowEdge(1);
728    max =  profXmax->GetBinLowEdge(nbin) + profXmax->GetBinWidth(nbin);
729    resolXmax = new TProfile("resolXmax","Normalized sigma of Xmax VS Energy",nbin,min,max,1e-6,1e6);
730    for(Int_t i=0; i<nbin; i++) {
731        if(profXmax->GetBinContent(i+1)) resolXmax->Fill(profXmax->GetBinCenter(i+1),profXmax->GetBinError(i+1)/profXmax->GetBinContent(i+1));  // sigma/mean
732    }
733    resolXmax->SetMarkerStyle(27);
734    resolXmax->SetXTitle("log10(Energy) (MeV)");
735    resolXmax->Draw();
736   
737    pad21->cd();
738    treepainter->DrawHistos("Nemax_Energy_linear",drawcut);
739    prof = (TProfile*)gDirectory->Get("h_Nemax_Energy_linear");
740    nbin = prof->GetNbinsX();
741    min = prof->GetBinLowEdge(1);
742    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
743    resolNemax = new TProfile("resolNemax","Normalized sigma of Nemax VS Energy",nbin,min,max,1e-6,1e6);
744    for(Int_t i=0; i<nbin; i++) {
745        if(prof->GetBinContent(i+1)) resolNemax->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
746    }
747    resolNemax->SetXTitle("log10(Energy) (MeV)");
748    resolNemax->SetMarkerStyle(27);
749    resolNemax->Draw();
750 
751    pad31->cd();
752    treepainter->DrawHistos("Ne2_Energy_linear",drawcut);
753    prof = (TProfile*)gDirectory->Get("h_Ne2_Energy_linear");
754    nbin = prof->GetNbinsX();
755    min = prof->GetBinLowEdge(1);
756    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
757    resolNe2 = new TProfile("resolNe2","#splitline{Normalized sigma of \"Ne longitudinal}{profile first HALF\" VS Energy}",nbin,min,max,1e-6,1e6);
758    for(Int_t i=0; i<nbin; i++) {
759        if(prof->GetBinContent(i+1)) resolNe2->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
760    }
761    gPad->SetTopMargin(0.2);
762    resolNe2->SetXTitle("log10(Energy) (MeV)");
763    resolNe2->SetMarkerStyle(27);
764    resolNe2->SetTitleOffset(2);
765    resolNe2->Draw();
766}
767
768//__________________________________________________________________________________________________________________________-
769void shower_theta_correl() {
770    Double_t min = 0;
771    Double_t max = 0;
772    Double_t textpos = 0;
773    Double_t par[3];
774    TLatex *t = 0;
775    Int_t nbin = 0;
776    TF1* fitfunc = 0;
777    TF1* mypol = 0;   
778
779    Bool_t flag = true;
780    if(!chain&&!etree&&!stree) flag = checktree();     
781    if(!flag) {
782        Warning("treeatmoEnergycorrel()","No Tree found");
783        return;
784    }
785   
786    TCanvas *cc = pad_shower->FindCanvas("c_treeshower_correltheta");
787    if(!cc ) cc = pad_shower->AddTab("c_treeshower_correltheta","c_treeshower_correltheta");
788
789
790    // *************** PLOTS FOR BEHAVIOR ********************
791
792    cc->Clear();
793    cc->cd();
794
795    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
796    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
797    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
798    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
799
800    pad11->Draw();
801    pad21->Draw();
802    pad31->Draw();
803    pad41->Draw();
804
805    pad11->cd();
806    treepainter->DrawHistos("Ne_Theta",drawcut);
807    histo = (TH1*)gPad->GetPrimitive("h_Ne_Theta");
808    histo->SetTitle("Electrons profile integration vs/ theta -- energy normalized");
809    histo->SetYTitle("Electron profile integration    g/(MeV.cm^{2})");
810    histo->SetXTitle("#theta (deg)");
811    histo->SetStats(0);
812       
813    pad21->cd();
814    treepainter->DrawHistos("Ne2_Theta",drawcut);
815    histo = (TH1*)gPad->GetPrimitive("h_Ne2_Theta");
816    histo->SetTitle("HALF electrons profile integration vs/ theta -- energy normalized");
817    histo->SetYTitle("HALF lectron profile integration    g/(MeV.cm^{2})");
818    histo->SetXTitle("#theta (deg)");
819    histo->SetStats(0);
820   
821    pad31->cd();
822    treepainter->DrawHistos("Nemax_Theta",drawcut);
823    prof_widthTheta = (TProfile*)gPad->GetPrimitive("h_Nemax_Theta");
824    prof_widthTheta->SetTitle("Nemax vs/ #theta -- energy normalized");
825    prof_widthTheta->SetXTitle("#theta (deg)");
826    prof_widthTheta->SetYTitle("Nemax  (Mev^{-1})");
827    prof_widthTheta->SetStats(0);
828   
829    pad41->cd();
830    treepainter->DrawHistos("Widthshow_Hmax",drawcut);
831    prof_widthTheta = (TProfile*)gPad->GetPrimitive("h_Widthshow_Hmax");
832    prof_widthTheta->SetTitle("Shower width vs/ Hmax");
833    prof_widthTheta->SetXTitle("Hmax  (km)");
834    prof_widthTheta->SetYTitle("Width  (km)");
835    prof_widthTheta->SetStats(0);
836}
837
838//__________________________________________________________________________________________________________________________-
839void phot_elec_correl() {
840    //
841    //
842    //
843   
844    Double_t min = 0;
845    Double_t max = 0;
846    Double_t textpos = 0;
847    Double_t par[3];
848    TLatex *t = 0;
849    Int_t nbin = 0;
850    TF1* fitfunc = 0;
851
852    TCanvas *c = pad_atmo_photons->FindCanvas("c_treeshowphot_correl");
853    if(!c ) c = pad_atmo_photons->AddTab("c_treeshowphot_correl","c_treeshowphot_correl");
854    TCanvas *c2 = pad_atmo_photons->FindCanvas("c_treeshowphot_correl_F");
855    if(!c2 ) c2 = pad_atmo_photons->AddTab("c_treeshowphot_correl_F","c_treeshowphot_correl_F");
856
857   
858    // *************** PLOTS FOR BEHAVIOR ********************
859
860    c->Clear();
861    c->cd();
862    Bool_t flag = true;
863    if(!chain&&!etree&&!stree) flag = checktree();     
864    if(!flag) {
865        delete c; c= 0;
866        Warning("treeshowphotcorrel()","No Tree found");
867        return;
868    }
869
870    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
871    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
872    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
873    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
874
875    pad1->Draw();
876    pad2->Draw();
877    pad3->Draw();
878    pad4->Draw();
879
880    pad1->cd();
881    treepainter->DrawHistos("Nph2fluo_Ne2",drawcut);
882    prof_Nph2fluoNe = (TProfile*)gPad->GetPrimitive("h_Nph2fluo_Ne2");
883    prof_Nph2fluoNe->SetTitle("HALF Fluo vs/ integration of Ne profile");
884    prof_Nph2fluoNe->SetXTitle("(log10) HALF integration of Ne profile (g/cm^{2})");
885    prof_Nph2fluoNe->SetYTitle("(log10) HALF fluo photons generated");
886    prof_Nph2fluoNe->SetStats(0);
887    min = 0;
888    max = 0;
889    for(Int_t m=0; m<prof_Nph2fluoNe->GetNbinsX(); m++) {
890        if(prof_Nph2fluoNe->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nph2fluoNe->GetBinLowEdge(m+1);
891        if(prof_Nph2fluoNe->GetBinContent(prof_Nph2fluoNe->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nph2fluoNe->GetBinLowEdge(prof_Nph2fluoNe->GetNbinsX() - m);
892        if(fabs(min*max) > 0) break;
893    }
894    mypol = new TF1("mypol_Nph2fluo_Ne2",my_linear,min,max,2);  // pente fixee a 1
895    mypol->SetLineWidth(1);
896    mypol->SetLineColor(2);
897    mypol->SetParameter(1.,1.);
898    mypol->SetParLimits(0,-10.,10.);
899    mypol->SetParLimits(1,1.,1.);
900    textpos = min + 0.1*(max-min);
901    prof_Nph2fluoNe->Fit("mypol_Nph2fluo_Ne2","WRCOQ");
902    mypol->GetParameters(par);
903    t = new TLatex();
904    char txt[200];
905    t->SetTextFont(62); 
906    t->SetTextSize(0.06);
907    t->SetTextColor(2);
908    t->SetTextColor(1);
909    sprintf(txt,"Nfluo_{1/2} = %.1f  *  Ne_{1/2}",pow(10,par[0]));
910    t->DrawLatex(textpos,prof_Nph2fluoNe->GetMaximum()*0.54,txt);
911
912    pad2->cd();
913    treepainter->DrawHistos("Nph2ckov_Ne2",drawcut);
914    prof_Nph2ckovNe = (TProfile*)gPad->GetPrimitive("h_Nph2ckov_Ne2");
915    prof_Nph2ckovNe->SetTitle("HALF Cerenkov vs/ HALF integration of Ne profile");
916    prof_Nph2ckovNe->SetXTitle("(log10) HALF integration Ne profile (g/cm^{2})");
917    prof_Nph2ckovNe->SetYTitle("(log10) HALF cerenkov photons generated");
918    prof_Nph2ckovNe->SetStats(0);
919    min = 0;
920    max = 0;
921    for(Int_t m=0; m<prof_Nph2ckovNe->GetNbinsX(); m++) {
922        if(prof_Nph2ckovNe->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nph2ckovNe->GetBinLowEdge(m+1);
923        if(prof_Nph2ckovNe->GetBinContent(prof_Nph2ckovNe->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nph2ckovNe->GetBinLowEdge(prof_Nph2ckovNe->GetNbinsX() - m);
924        if(fabs(min*max) > 0) break;
925    }
926    mypol = new TF1("mypol_Nph2ckov_Ne2",my_linear,min,max,2);  // pente fixee a 1
927    mypol->SetLineWidth(1);
928    mypol->SetLineColor(2);
929    mypol->SetParameter(1.,1.);
930    mypol->SetParLimits(0,-10.,10.);
931    mypol->SetParLimits(1,1.,1.);
932    textpos = min + 0.1*(max-min);
933    prof_Nph2ckovNe->Fit("mypol_Nph2ckov_Ne2","WRCOQ");
934    mypol->GetParameters(par);
935    t = new TLatex();
936    char txt[200];
937    t->SetTextFont(62); 
938    t->SetTextSize(0.06);
939    t->SetTextColor(2);
940    t->SetTextColor(1);
941    sprintf(txt,"Nckov_{1/2} = %.1f  *  Ne_{1/2}",pow(10,par[0]));
942    t->DrawLatex(textpos,prof_Nph2ckovNe->GetMaximum()*0.54,txt);
943
944    pad3->cd();
945    treepainter->DrawHistos("Nphmaxfluo_Nemax",drawcut);
946    prof_Nphmax = (TProfile*)gPad->GetPrimitive("h_Nphmaxfluo_Nemax");
947    prof_Nphmax->SetTitle("Maximum of fluo profile vs/ Maximum of Ne profile");
948    prof_Nphmax->SetXTitle("log10(Nemax)");
949    prof_Nphmax->SetYTitle("(log10) Maximum of fluo profile");
950    prof_Nphmax->SetStats(0);
951    min = 0;
952    max = 0;
953    for(Int_t m=0; m<prof_Nphmax->GetNbinsX(); m++) {
954        if(prof_Nphmax->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nphmax->GetBinLowEdge(m+1);
955        if(prof_Nphmax->GetBinContent(prof_Nphmax->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nphmax->GetBinLowEdge(prof_Nphmax->GetNbinsX() - m);
956        if(fabs(min*max) > 0) break;
957    }
958    mypol = new TF1("mypol_Nphmaxfluo_Nemax",my_linear,min,max,2);  // pente fixee a 1
959    mypol->SetLineWidth(1);
960    mypol->SetLineColor(2);
961    mypol->SetParameter(1.,1.);
962    mypol->SetParLimits(0,-10.,10.);
963    mypol->SetParLimits(1,1.,1.);
964    textpos = min + 0.1*(max-min);
965    prof_Nphmax->Fit("mypol_Nphmaxfluo_Nemax","WRCOQ");
966    mypol->GetParameters(par);
967    t = new TLatex();
968    char txt[200];
969    t->SetTextFont(62); 
970    t->SetTextSize(0.06);
971    t->SetTextColor(2);
972    t->SetTextColor(1);
973    sprintf(txt,"Nfluo_{max} = %.1f  *  Ne_{max}  cm^{2}/g",pow(10,par[0]));
974    t->DrawLatex(textpos,prof_Nphmax->GetMaximum()*0.54,txt);
975
976    pad4->cd();
977    treepainter->DrawHistos("Nph2fluo_Nemax",drawcut);
978    prof_NphfluoNemax = (TProfile*)gPad->GetPrimitive("h_Nph2fluo_Nemax");
979    prof_NphfluoNemax->SetTitle("HALF Fluo vs/ Maximum of Ne profile");
980    prof_NphfluoNemax->SetXTitle("log10(Nemax)");
981    prof_NphfluoNemax->SetYTitle("(log10) HALF fluo generated");
982    prof_NphfluoNemax->SetStats(0);
983    min = 0;
984    max = 0;
985    for(Int_t m=0; m<prof_NphfluoNemax->GetNbinsX(); m++) {
986        if(prof_NphfluoNemax->GetBinContent(m+1) && fabs(min) == 0) min = prof_NphfluoNemax->GetBinLowEdge(m+1);
987        if(prof_NphfluoNemax->GetBinContent(prof_NphfluoNemax->GetNbinsX() - m) && fabs(max) == 0) max = prof_NphfluoNemax->GetBinLowEdge(prof_NphfluoNemax->GetNbinsX() - m);
988        if(fabs(min*max) > 0) break;
989    }
990    mypol = new TF1("mypol_Nph2fluo_Nemax",my_linear,min,max,2);  // pente fixee a 1
991    mypol->SetLineWidth(1);
992    mypol->SetLineColor(2);
993    mypol->SetParameter(1.,1.);
994    mypol->SetParLimits(0,-10.,10.);
995    mypol->SetParLimits(1,1.,1.);
996    textpos = min + 0.1*(max-min);
997    prof_NphfluoNemax->Fit("mypol_Nph2fluo_Nemax","WRCOQ");
998    mypol->GetParameters(par);
999    t = new TLatex();
1000    char txt[200];
1001    t->SetTextFont(62); 
1002    t->SetTextSize(0.06);
1003    t->SetTextColor(2);
1004    t->SetTextColor(1);
1005    sprintf(txt,"Nfluo_{1/2} = %.1f  *  Ne_{max}",pow(10,par[0]));
1006    t->DrawLatex(textpos,prof_NphfluoNemax->GetMaximum()*0.54,txt);
1007
1008   
1009    // *************** PLOTS FOR SIGMA  ***********************************
1010   
1011    c2->Clear();
1012    c2->cd();
1013    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1014    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1015    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1016    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1017
1018    pad11->Draw();
1019    pad21->Draw();
1020    pad31->Draw();
1021    pad41->Draw();
1022    TProfile* prof = 0;
1023   
1024    pad11->cd();
1025    treepainter->DrawHistos("Nph2fluo_Ne2_linear",drawcut);
1026    prof = (TProfile*)gDirectory->Get("h_Nph2fluo_Ne2_linear");
1027    nbin = prof->GetNbinsX();
1028    min = prof->GetBinLowEdge(1);
1029    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1030    resolNph2FluoNe = new TProfile("resolNph2FluoNe","#splitline{Normalized sigma of Fluo first HALF VS/}{Ne first HALF profile}",nbin,min,max,1e-6,1e6);
1031    for(Int_t i=0; i<nbin; i++)
1032        if(prof->GetBinContent(i+1)) resolNph2FluoNe->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1033    resolNph2FluoNe->SetMarkerStyle(27);
1034    resolNph2FluoNe->SetXTitle("(log10) HALF integration of Ne profile (g/cm^{2})");
1035    resolNph2FluoNe->Draw();
1036    gPad->SetTopMargin(0.18);
1037   
1038    pad21->cd();
1039    treepainter->DrawHistos("Nph2ckov_Ne2_linear",drawcut);
1040    prof = (TProfile*)gDirectory->Get("h_Nph2ckov_Ne2_linear");
1041    nbin = prof->GetNbinsX();
1042    min = prof->GetBinLowEdge(1);
1043    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1044    resolNph2CkovNe = new TProfile("resolNph2CkovNe","#splitline{Normalized sigma of Cerenkov first HALF VS/}{Ne first HALF profile}",nbin,min,max,1e-6,1e6);
1045    for(Int_t i=0; i<nbin; i++) {
1046        if(prof->GetBinContent(i+1)) resolNph2CkovNe->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1047    }
1048    resolNph2CkovNe->SetMarkerStyle(27);
1049    resolNph2CkovNe->SetXTitle("(log10) HALF integration of Ne profile (g/cm^{2})");
1050    resolNph2CkovNe->Draw();
1051    gPad->SetTopMargin(0.18);
1052 
1053    pad31->cd();
1054    treepainter->DrawHistos("Nphmaxfluo_Nemax_linear",drawcut);
1055    prof = (TProfile*)gDirectory->Get("h_Nphmaxfluo_Nemax_linear");
1056    nbin = prof->GetNbinsX();
1057    min = prof->GetBinLowEdge(1);
1058    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1059    resolNmaxFluoNemax = new TProfile("resolNmaxFluoNemax","#splitline{Normalized sigma of max Fluo VS/}{Max of electrons}",nbin,min,max,1e-6,1e6);
1060    for(Int_t i=0; i<nbin; i++) {
1061        if(prof->GetBinContent(i+1)) resolNmaxFluoNemax->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1062    }
1063    resolNmaxFluoNemax->SetMarkerStyle(27);
1064    resolNmaxFluoNemax->SetXTitle("log10(Nemax)");
1065    resolNmaxFluoNemax->Draw();
1066    gPad->SetTopMargin(0.18);
1067     
1068    pad41->cd();
1069    treepainter->DrawHistos("Nph2fluo_Nemax_linear",drawcut);
1070    prof = (TProfile*)gDirectory->Get("h_Nph2fluo_Nemax_linear");
1071    nbin = prof->GetNbinsX();
1072    min = prof->GetBinLowEdge(1);
1073    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1074    resolNphFluoNemax = new TProfile("resolNphFluoNemax","#splitline{Normalized sigma of Fluo first HALF  VS/}{Max of electrons}",nbin,min,max,1e-6,1e6);
1075    for(Int_t i=0; i<nbin; i++) {
1076        if(prof->GetBinContent(i+1)) resolNphFluoNemax->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1077    }
1078    resolNphFluoNemax->SetMarkerStyle(27);
1079    resolNphFluoNemax->SetXTitle("log10(Nemax)");
1080    resolNphFluoNemax->Draw();
1081    gPad->SetTopMargin(0.18);
1082}
1083
1084//__________________________________________________________________________________________________________________________-
1085void crash_ground() {
1086   
1087    Double_t min = 0;
1088    Double_t max = 0;
1089    Double_t textpos = 0;
1090    Double_t par[3];
1091    TLatex *t = 0;
1092    Int_t nbin = 0;
1093    TF1* fitfunc = 0;
1094   
1095    TCanvas *cc = pad_atmo_photons->FindCanvas("c_treeshower_crashGround");
1096    if(!cc ) cc = pad_atmo_photons->AddTab("c_treeshower_crashGround","c_treeshower_crashGround");
1097
1098    cc->Clear();
1099    cc->cd();
1100
1101    TPad *pad10 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1102    TPad *pad20 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1103    TPad *pad30 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1104    TPad *pad40 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1105
1106    pad10->Draw();
1107    pad20->Draw();
1108    pad30->Draw();
1109    pad40->Draw();
1110
1111    pad10->cd();
1112    treepainter->DrawHistos("Ne2_Over_Ne_Hmax",drawcut);
1113    histo = (TH1*)gPad->GetPrimitive("h_Ne2_Over_Ne_Hmax");
1114    histo->SetTitle("Shower crash effect on electrons profile");
1115    histo->SetXTitle("Local altitude of shower maximum (km)");
1116    histo->SetYTitle("Ratio ( first HALF / FULL profile )");
1117    histo->GetYaxis()->SetTitleOffset(1.5);
1118    gPad->SetLeftMargin(0.13);
1119    histo->SetStats(0);
1120   
1121    pad20->cd();
1122    treepainter->DrawHistos("Nph2_Over_Nph_Hmax_fluo",drawcut);
1123    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Over_Nph_Hmax_fluo");
1124    histo->SetTitle("Shower crash effect on fluo photons nb");
1125    histo->SetXTitle("Local altitude of shower maximum (km)");
1126    histo->SetYTitle("Ratio ( first HALF / FULL profiles )");
1127    histo->SetStats(0);
1128   
1129    pad30->cd();
1130    treepainter->DrawHistos("Width2_Over_Width_Hmax_shower",drawcut);
1131    histo = (TH1*)gPad->GetPrimitive("h_Width2_Over_Width_Hmax_shower");
1132    histo->SetTitle("Crash effect on shower Width");
1133    histo->SetXTitle("Local altitude of shower maximum (km)");
1134    histo->SetYTitle("Ratio ( first HALF width / FULL width )");
1135    histo->SetStats(0);
1136   
1137    pad40->cd();
1138    treepainter->DrawHistos("Width2_Over_Width_Hmax_fluo",drawcut);
1139    histo = (TH1*)gPad->GetPrimitive("h_Width2_Over_Width_Hmax_fluo");
1140    histo->SetTitle("Shower crash effect on fluo profile Width");
1141    histo->SetXTitle("Local altitude of shower maximum (km)");
1142    histo->SetYTitle("Ratio ( first HALF width / FULL width )");
1143    histo->SetStats(0);
1144   
1145    /* for CKOV : not really used so far
1146    treepainter->DrawHistos("Nph2_Over_Nph_Hmax_cer",drawcut);
1147    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Over_Nph_Hmax_cer");
1148    histo->SetTitle("Shower crash effect on cerenkov photons nb");
1149    histo->SetXTitle("Local altitude of shower maximum (km)");
1150    histo->SetYTitle("Ratio ( first HALF / FULL profiles )");
1151    histo->SetStats(0);
1152    */
1153   
1154//////////////////////////////////////////
1155   
1156    TCanvas *c = pad_atmo_photons->FindCanvas("c_treeshower_crashGround2");
1157    if(!c ) c = pad_atmo_photons->AddTab("c_treeshower_crashGround2","c_treeshower_crashGround2");
1158
1159    c->Clear();
1160    c->cd();
1161
1162    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1163    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1164    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1165    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1166
1167    pad1->Draw();
1168    pad2->Draw();
1169    pad3->Draw();
1170    pad4->Draw();
1171   
1172    pad1->cd();
1173    treepainter->DrawHistos("Nph2_Over_Nph_Theta_fluo",drawcut);
1174    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Over_Nph_Theta_fluo");
1175    histo->SetTitle("Shower crash effect on fluo photons nb");
1176    histo->SetXTitle("#theta (deg)");
1177    histo->SetYTitle("Ratio ( first HALF / FULL profiles )");
1178    histo->SetStats(0);
1179   
1180    pad2->cd();
1181    treepainter->DrawHistos("Nph2_Over_Nph_Theta_cer",drawcut);
1182    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Over_Nph_Theta_cer");
1183    histo->SetTitle("Shower crash effect on cerenkov photons nb");
1184    histo->SetXTitle("#theta (deg)");
1185    histo->SetYTitle("Ratio ( first HALF / FULL profiles )");
1186    histo->SetStats(0);
1187   
1188    pad3->cd();
1189    treepainter->DrawHistos("Ne2_Over_Ne_Theta",drawcut);
1190    histo = (TH1*)gPad->GetPrimitive("h_Ne2_Over_Ne_Theta");
1191    histo->SetTitle("Shower crash effect on electrons profile");
1192    histo->SetXTitle("#theta (deg)");
1193    histo->SetYTitle("Ratio ( first HALF / FULL profile )");
1194    histo->GetYaxis()->SetTitleOffset(1.5);
1195    gPad->SetLeftMargin(0.13);
1196    histo->SetStats(0);
1197   
1198    pad4->cd();
1199    treepainter->DrawHistos("Width2_Over_Width_Theta_fluo",drawcut);
1200    histo = (TH1*)gPad->GetPrimitive("h_Width2_Over_Width_Theta_fluo");
1201    histo->SetTitle("Shower crash effect on fluo profile Width");
1202    histo->SetXTitle("#theta (deg)");
1203    histo->SetYTitle("Ratio ( first HALF width / FULL width )");
1204    histo->SetStats(0);
1205   
1206
1207}
1208
1209//__________________________________________________________________________________________________________________________-
1210void Xmax_shift() {
1211   
1212    Double_t min = 0;
1213    Double_t max = 0;
1214    Double_t textpos = 0;
1215    Double_t par[3];
1216    TLatex *t = 0;
1217    Int_t nbin = 0;
1218    TF1* fitfunc = 0;
1219   
1220    TCanvas *cc2 = pad_atmo_photons->FindCanvas("c_treeshower_Xmaxshift");
1221    if(!cc2 ) cc2 = pad_atmo_photons->AddTab("c_treeshower_Xmaxshift","c_treeshower_Xmaxshift");
1222
1223    cc2->Clear();
1224    cc2->cd();
1225
1226    TPad *pad102 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1227    TPad *pad202 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1228    TPad *pad302 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1229    TPad *pad402 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1230
1231    pad102->Draw();
1232    pad202->Draw();
1233    pad302->Draw();
1234    pad402->Draw();
1235
1236    pad102->cd();
1237    treepainter->DrawHistos("Xmaxfluo_Xmax_Theta",drawcut);
1238    histo = (TH1*)gPad->GetPrimitive("h_Xmaxfluo_Xmax_Theta");
1239    histo->SetTitle("Xmax shift between fluo and electrons VS/ #theta");
1240    histo->SetXTitle("#theta (deg)");
1241    histo->SetYTitle("X_{max}Fluo - X_{max}Shower (g/cm^{2})");
1242   
1243    pad202->cd();
1244    treepainter->DrawHistos("Xmaxfluo_Xmax_Energy",drawcut);
1245    histo = (TH1*)gPad->GetPrimitive("h_Xmaxfluo_Xmax_Energy");
1246    histo->SetTitle("Xmax shift between fluo and electrons VS/ Energy");
1247    histo->SetXTitle("log10(Energy) (MeV)");
1248    histo->SetYTitle("Xmax_fluo - Xmax_shower (g/cm^{2})");
1249   
1250    pad302->cd();
1251    treepainter->DrawHistos("Xmaxfluo_Xmax_X1",drawcut);
1252    histo = (TH1*)gPad->GetPrimitive("h_Xmaxfluo_Xmax_X1");
1253    histo->SetTitle("Xmax shift between fluo and electrons VS/ X1");
1254    histo->SetXTitle("X1 (g/cm^{2})");
1255    histo->SetYTitle("Xmax_fluo - Xmax_shower (g/cm^{2})");
1256   
1257    pad402->cd();
1258    treepainter->DrawHistos("Xmaxfluo_Xmax_Hmax",drawcut);
1259    histo = (TH1*)gPad->GetPrimitive("h_Xmaxfluo_Xmax_Hmax");
1260    histo->SetTitle("Xmax shift between fluo and electrons VS/ Hmax");
1261    histo->SetXTitle("Hmax  (km)");
1262    histo->SetYTitle("X_{max}Fluo - X_{max}Shower (g/cm^{2})");
1263
1264
1265    TCanvas * c = pad_atmo_photons->FindCanvas("c_treeshower_Xmaxshift2");
1266    if(!c )   c = pad_atmo_photons->AddTab("c_treeshower_Xmaxshift2","c_treeshower_Xmaxshift2");
1267
1268    c->Clear();
1269    c->cd();
1270
1271    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1272    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1273    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1274    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1275
1276    pad1->Draw();
1277    pad2->Draw();
1278    pad3->Draw();
1279    pad4->Draw();
1280
1281    pad1->cd();
1282    treepainter->DrawHistos("Hmaxfluo_Hmax_Theta",drawcut);
1283    histo = (TH1*)gPad->GetPrimitive("h_Hmaxfluo_Hmax_Theta");
1284    histo->SetTitle("Hmax shift between fluo and electrons VS #theta");
1285    histo->SetXTitle("#theta (deg)");
1286    histo->SetYTitle("Hmax_fluo - Hmax_shower (km)");
1287   
1288    pad2->cd();
1289    treepainter->DrawHistos("Lmaxfluo_Lmax_Theta",drawcut);
1290    histo = (TH1*)gPad->GetPrimitive("h_Lmaxfluo_Lmax_Theta");
1291    histo->SetTitle("Lmax shift between fluo and electrons VS #theta");
1292    histo->SetXTitle("#theta (deg)");
1293    histo->SetYTitle("Lmax_fluo - Lmax_shower (km)");
1294   
1295    pad3->cd();
1296    treepainter->DrawHistos("xmax_ymax_shift",drawcut);
1297    histo = (TH1*)gPad->GetPrimitive("h_xmax_ymax_shift");
1298    histo->SetTitle("Position shift between fluo and electrons");
1299    histo->SetXTitle("Shift in x-coordinate");
1300    histo->SetYTitle("Shift in y-coordinate");
1301   
1302    pad4->cd();
1303    treepainter->DrawHistos("Rmaxfluo_Rmax_Theta",drawcut);
1304    histo = (TH1*)gPad->GetPrimitive("h_Rmaxfluo_Rmax_Theta");
1305    histo->SetTitle("Rmax shift between fluo and electrons VS #theta");
1306    histo->SetXTitle("#theta (deg)");
1307    histo->SetYTitle("Rmax_fluo - Rmax_shower (km)");
1308   
1309    TCanvas * c1 = pad_atmo_photons->FindCanvas("c_treeshower_Xmaxshift3");
1310    if(!c1 )   c1 = pad_atmo_photons->AddTab("c_treeshower_Xmaxshift3","c_treeshower_Xmaxshift3");
1311
1312    c1->Clear();
1313    c1->cd();
1314   
1315    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1316    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1317    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1318    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1319
1320    pad11->Draw();
1321    pad21->Draw();
1322    pad31->Draw();
1323    pad41->Draw();
1324
1325    pad11->cd();
1326    treepainter->DrawHistos("Hmaxcer_Hmax_Theta",drawcut);
1327    histo = (TH1*)gPad->GetPrimitive("h_Hmaxcer_Hmax_Theta");
1328    histo->SetTitle("Hmax shift between ckov and electrons VS #theta");
1329    histo->SetXTitle("#theta (deg)");
1330    histo->SetYTitle("Hmax_ckov - Hmax_shower (km)");
1331   
1332    pad21->cd();
1333    treepainter->DrawHistos("Yield_mean_Hmax",drawcut);
1334    histo = (TH1*)gPad->GetPrimitive("h_Yield_mean_Hmax");
1335    histo->SetTitle("Mean Fluo yield");
1336    histo->SetXTitle("Hmax_show  (km)");
1337    histo->SetYTitle("Yield  (m_{1})");
1338   
1339    pad31->cd();
1340    treepainter->DrawHistos("Yield_max_Hmax",drawcut);
1341    histo = (TH1*)gPad->GetPrimitive("h_Yield_max_Hmax");
1342    histo->SetTitle("Fluo yield at shower maximum");
1343    histo->SetXTitle("Hmax_show  (km)");
1344    histo->SetYTitle("Yield  (m_{1})");
1345   
1346    pad41->cd();
1347    treepainter->DrawHistos("Yield_max_cer_Hmax",drawcut);
1348    histo = (TH1*)gPad->GetPrimitive("h_Yield_max_cer_Hmax");
1349    histo->SetTitle("Ckov yield at shower maximum");
1350    histo->SetXTitle("Hmax_show  (km)");
1351    histo->SetYTitle("Yield  (m_{1})");
1352
1353}
1354
1355//__________________________________________________________________________________________________________________________-
1356void treeatmohistos() {
1357    Bool_t flag = true;
1358    if(!chain&&!etree&&!stree) flag = checktree();     
1359    if(!flag) {
1360        Warning("treeatmohistos()","No Tree found");
1361        return;
1362    }
1363   
1364    Double_t min = 0;
1365    Double_t max = 0;
1366    Double_t textpos = 0;
1367    Double_t par[3];
1368    TLatex *t = 0;
1369   
1370    if(!pad_atmo_photons) pad_atmo_photons = new EGViewer();
1371    TCanvas *c = pad_atmo_photons->FindCanvas("c_treeatmo_histo");
1372    if(!c ) c = pad_atmo_photons->AddTab("c_treeatmo_histo","c_treeatmo_histo");
1373    TCanvas *c2 = pad_atmo_photons->FindCanvas("c_treeatmo_histo2");
1374    if(!c2 ) c2 = pad_atmo_photons->AddTab("c_treeatmo_histo2","c_treeatmo_histo2");
1375    //TCanvas *c3 = pad_atmo_photons->FindCanvas("c_treeatmo_histo3");
1376    //if(!c3 ) c3 = pad_atmo_photons->AddTab("c_treeatmo_histo3","c_treeatmo_histo3");
1377
1378    c->Clear();
1379    c->cd();
1380
1381    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1382    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1383    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1384    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1385
1386    pad1->Draw();
1387    pad2->Draw();
1388    pad3->Draw();
1389    pad4->Draw();
1390
1391    pad1->cd();
1392    pad1->SetLogy(1);
1393    treepainter->DrawHistos("Nph_Fluo_L",drawcut);
1394    histo = (TH1*)gPad->GetPrimitive("h_Nph_Fluo_L");
1395    histo->SetTitle("Fluorescence photons at creation");
1396    histo->SetXTitle("Nph_fluo");
1397   
1398    pad2->cd();
1399    pad2->SetLogy(1);
1400    treepainter->DrawHistos("Nph_Cer_L",drawcut);
1401    histo = (TH1*)gPad->GetPrimitive("h_Nph_Cer_L");
1402    histo->SetTitle("Cerenkov photons at creation");
1403    histo->SetXTitle("Nph_cer");
1404   
1405    pad3->cd();
1406    pad3->SetLogy(1);
1407    treepainter->DrawHistos("Nmax_Fluo",drawcut);
1408    histo = (TH1*)gPad->GetPrimitive("h_Nmax_Fluo");
1409    histo->SetTitle("Fluo photons maximum");
1410    histo->SetXTitle("Nphmax_fluo (cm^{2}/g -- per depth unit)");
1411   
1412    /*
1413    pad4->cd();
1414    pad4->SetLogy(1);
1415    treepainter->DrawHistos("Nmax_Cer",drawcut);
1416    histo = (TH1*)gPad->GetPrimitive("h_Nmax_Cer");
1417    histo->SetTitle("Cerenkov photons maximum");
1418    histo->SetXTitle("Nphmax_cer (cm^{2}/g -- per depth unit)");
1419    */
1420   
1421   
1422////////////////////////////////
1423    c2->Clear();
1424    c2->cd();
1425
1426    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1427    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1428    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1429    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1430
1431    pad11->Draw();
1432    pad21->Draw();
1433    pad31->Draw();
1434    pad41->Draw();
1435
1436    pad11->cd();
1437    pad11->SetLogy(1);
1438    treepainter->DrawHistos("Xmaxph_Fluo",drawcut);
1439    histo = (TH1*)gPad->GetPrimitive("h_Xmaxph_Fluo");
1440    histo->SetTitle("Depth of Fluorescence maximum");
1441    histo->SetXTitle("Xmaxph_fluo (g/cm^{2})");
1442    min = histo->GetBinCenter(histo->GetMaximumBin()) + 10;
1443    max = histo->GetBinLowEdge(histo->GetXaxis()->GetLast()) + histo->GetBinWidth(histo->GetXaxis()->GetLast());
1444    mypol = new TF1("mypol_Xmax_fluo","expo",min,max);   // fit the Xmax_fluo distribution tail (10 g/cm2 after maximum)
1445    mypol->SetLineWidth(1);
1446    mypol->SetLineColor(2);
1447    textpos = histo->GetBinLowEdge(1) + 0.05*(histo->GetBinLowEdge(histo->GetXaxis()->GetLast()) - histo->GetBinLowEdge(1));
1448    histo->Fit("mypol_Xmax_fluo","WOQR");
1449    mypol->GetParameters(par);
1450    t = new TLatex();
1451    char txt[200];
1452    t->SetTextFont(62); 
1453    t->SetTextSize(0.06);
1454    t->SetTextColor(2);
1455    t->SetTextColor(1);
1456    sprintf(txt,"%.2g  exp(%.2g Xfluo_{max})",exp(par[0]),par[1]);
1457    t->DrawLatex(textpos,pow(10,log10(histo->GetMaximum())*0.15),txt);
1458
1459   /* pad21->cd();
1460    pad21->SetLogy(1);
1461    treepainter->DrawHistos("Xmaxph_Cer",drawcut);
1462    histo = (TH1*)gPad->GetPrimitive("h_Xmaxph_Cer");
1463    histo->SetTitle("Depth of Cerenkov maximum");
1464    histo->SetXTitle("Xmaxph_cer (g/cm^{2})");
1465    */
1466    pad31->cd();
1467    treepainter->DrawHistos("Hmax_Fluo",drawcut);
1468    histo = (TH1*)gPad->GetPrimitive("h_Hmax_Fluo");
1469    histo->SetTitle("Height of Fluorescence maximum");
1470    histo->SetXTitle("Hmax_fluo (km)");
1471    pad41->cd();
1472    treepainter->DrawHistos("Hmax_Cer",drawcut);
1473    histo = (TH1*)gPad->GetPrimitive("h_Hmax_Cer");
1474    histo->SetTitle("Height of Cerenkov maximum");
1475    histo->SetXTitle("Hmax_cer (km)");
1476   
1477   
1478/////////////////////////
1479/* c3 disabled too
1480    c3->Clear();
1481    c3->cd();
1482    TPad *pad12 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1483    TPad *pad22 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1484    TPad *pad32 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1485    TPad *pad42 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1486
1487    pad12->Draw();
1488    pad22->Draw();
1489    pad32->Draw();
1490    pad42->Draw();
1491
1492    pad12->cd();
1493    treepainter->DrawHistos("Wlmean_f",drawcut);
1494    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_f");
1495    histo->SetTitle("Fluorescence mean wavelength at creation");
1496    histo->SetXTitle("Wlmean_f (nm)");
1497    pad22->cd();
1498    treepainter->DrawHistos("Wlmean_c",drawcut);
1499    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_c");
1500    histo->SetTitle("Cerenkov mean wavelength at creation");
1501    histo->SetXTitle("Wlmean_c (nm)");
1502    pad32->cd();
1503    //treepainter->DrawHistos("Wlmean_tot",drawcut);
1504    //histo = (TH1*)gPad->GetPrimitive("h_Wlmean_tot");
1505    //histo->SetTitle("Photons mean wavelength at creation");
1506    //histo->SetXTitle("Wlmean_tot (nm)");
1507*/
1508   
1509}
1510
1511//__________________________________________________________________________________________________________________________-
1512void atmo_theta_correl() {
1513    Double_t min = 0;
1514    Double_t max = 0;
1515    Double_t textpos = 0;
1516    Double_t par[3];
1517    TLatex *t = 0;
1518    Int_t nbin = 0;
1519    TF1* fitfunc = 0;
1520    TF1* mypol = 0;   
1521
1522    Bool_t flag = true;
1523    if(!chain&&!etree&&!stree) flag = checktree();     
1524    if(!flag) {
1525        Warning("treeatmoEnergycorrel()","No Tree found");
1526        return;
1527    }
1528   
1529    TCanvas *cc = pad_atmo_photons->FindCanvas("c_treeatmo_correltheta");
1530    if(!cc ) cc = pad_atmo_photons->AddTab("c_treeatmo_correltheta","c_treeatmo_correltheta");
1531    TCanvas *cc2 = pad_atmo_photons->FindCanvas("c_treeatmo_correltheta_ckov");
1532    if(!cc2 ) cc2 = pad_atmo_photons->AddTab("c_treeatmo_correltheta_ckov","c_treeatmo_correltheta_ckov");
1533
1534
1535    // *************** PLOTS FOR BEHAVIOR ********************
1536
1537    cc->Clear();
1538    cc->cd();
1539
1540    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1541    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1542    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1543    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1544
1545    pad11->Draw();
1546    pad21->Draw();
1547    pad31->Draw();
1548    pad41->Draw();
1549
1550    pad11->cd();
1551    treepainter->DrawHistos("Nph2_Theta_fluo",drawcut);
1552    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Theta_fluo");
1553    histo->SetTitle("HALF Fluo photons vs/ theta -- energy normalized");
1554    histo->SetYTitle("HALF fluo photons generated");
1555    histo->SetXTitle("#theta (deg)");
1556    histo->SetStats(0);
1557    min = 0;
1558    max = 0;
1559    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
1560        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
1561        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
1562        if(fabs(min*max) > 0) break;
1563    }
1564    /*
1565    mypol = new TF1("mypol_Nph2_Theta_fluo",my_overcos60,min,max,1);
1566    mypol->SetLineWidth(1);
1567    mypol->SetLineColor(2);
1568    mypol->SetParameters(1.,0);
1569    mypol->SetParLimits(0,0.1,100.);
1570    textpos = min + 0.1*(max-min);
1571    histo->Fit("mypol_Nph2_Theta_fluo","WROQ");
1572    mypol->GetParameters(par);
1573    t = new TLatex();
1574    char txt[200];
1575    t->SetTextFont(62); 
1576    t->SetTextSize(0.06);
1577    t->SetTextColor(2);
1578    t->SetTextColor(1);
1579    sprintf(txt,"%.3g MeV^{-1}/ cos(#theta)",par[0]);
1580    t->DrawLatex(textpos,histo->GetMaximum()*0.74,txt);
1581    */
1582       
1583    pad21->cd();
1584    treepainter->DrawHistos("Nmax_Nph2_fluo_Hmax",drawcut);
1585    histo = (TH1*)gPad->GetPrimitive("h_Nmax_Nph2_fluo_Hmax");
1586    histo->SetTitle("ratio (fluo max #times1000 / Fluo HALF) vs/ Hmax_fluo");
1587    histo->SetYTitle("Fluo_{max} / Fluo_{HALF}  (#times 10^{-3})");
1588    histo->SetXTitle("Hmax fluo (km)");
1589    histo->SetStats(0);
1590   
1591    pad31->cd();
1592    treepainter->DrawHistos("Width2_fluo_Theta",drawcut);
1593    prof_widthTheta = (TProfile*)gPad->GetPrimitive("h_Width2_fluo_Theta");
1594    prof_widthTheta->SetTitle("Half width of Fluo profile vs/ #theta");
1595    prof_widthTheta->SetXTitle("#theta (deg)");
1596    prof_widthTheta->SetYTitle("HALF width of fluo profile (km)");
1597    prof_widthTheta->SetStats(0);
1598    min = 0;
1599    max = 0;
1600    for(Int_t m=0; m<prof_widthTheta->GetNbinsX(); m++) {
1601        if(prof_widthTheta->GetBinContent(m+1) && fabs(min) == 0) min = prof_widthTheta->GetBinLowEdge(m+1);
1602        if(prof_widthTheta->GetBinContent(prof_widthTheta->GetNbinsX() - m) && fabs(max) == 0) max = prof_widthTheta->GetBinLowEdge(prof_widthTheta->GetNbinsX() - m);
1603        if(fabs(min*max) > 0) break;
1604    }
1605    /*
1606    mypol = new TF1("mypol_Width2_fluo_Theta",my_overcos60,min,max,1);
1607    mypol->SetLineWidth(1);
1608    mypol->SetLineColor(2);
1609    mypol->SetParameters(0.01.,0);
1610    mypol->SetParLimits(0,1.,100.);
1611    textpos = min + 0.08*(max-min);
1612    prof_widthTheta->Fit("mypol_Width2_fluo_Theta","WROQ");
1613    mypol->GetParameters(par);
1614    t = new TLatex();
1615    char txt[200];
1616    t->SetTextFont(62); 
1617    t->SetTextSize(0.06);
1618    t->SetTextColor(2);
1619    t->SetTextColor(1);
1620    sprintf(txt,"W_{1/2} = %.3g km /cos(#theta)",par[0]);
1621    t->DrawLatex(textpos,prof_widthTheta->GetMaximum()*0.84,txt);
1622    */
1623   
1624   
1625    pad41->cd();
1626    treepainter->DrawHistos("Width2_fluo_Hmax_fluo",drawcut);
1627    prof_widthHmax = (TProfile*)gPad->GetPrimitive("h_Width2_fluo_Hmax_fluo");
1628    prof_widthHmax->SetTitle("Half width of Fluo profile vs/ fluo maximum height");
1629    prof_widthHmax->SetXTitle("Height of fluo max (km)");
1630    prof_widthHmax->SetYTitle("HALF width of fluo profile (km)");
1631    prof_widthHmax->SetStats(0);
1632    min = 0;
1633    max = 0;
1634    for(Int_t m=0; m<prof_widthHmax->GetNbinsX(); m++) {
1635        if(prof_widthHmax->GetBinContent(m+1) && fabs(min) == 0) min = prof_widthHmax->GetBinLowEdge(m+1);
1636        if(prof_widthHmax->GetBinContent(prof_widthHmax->GetNbinsX() - m) && fabs(max) == 0) max = prof_widthHmax->GetBinLowEdge(prof_widthHmax->GetNbinsX() - m);
1637        if(fabs(min*max) > 0) break;
1638    }
1639    /*
1640    mypol = new TF1("mypol_Width2_fluo_Hmax_fluo",my_expon15,2,max,2);
1641    mypol->SetLineWidth(1);
1642    mypol->SetLineColor(2);
1643    mypol->SetParameters(1.,0.01);
1644    mypol->SetParLimits(0,0.1,10.);
1645    mypol->SetParLimits(1,0.01,10.);
1646    textpos = min + 0.08*(max-min);
1647    prof_widthHmax->Fit("mypol_Width2_fluo_Hmax_fluo","ROQ");
1648    mypol->GetParameters(par);
1649    t = new TLatex();
1650    char txt[200];
1651    t->SetTextFont(62); 
1652    t->SetTextSize(0.06);
1653    t->SetTextColor(2);
1654    t->SetTextColor(1);
1655    sprintf(txt,"W_{1/2} = %.3g km * exp(H_{max}/%.3g)",exp(par[0]),1/par[1]);
1656    t->DrawLatex(textpos,prof_widthHmax->GetMaximum()*0.84,txt);
1657    */
1658
1659
1660////////////////////////// ckov
1661    cc2->Clear();
1662    cc2->cd();
1663
1664    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1665    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1666    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1667    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1668
1669    pad1->Draw();
1670    pad2->Draw();
1671    pad3->Draw();
1672    pad4->Draw();
1673
1674    pad1->cd();
1675    treepainter->DrawHistos("Nmax_Hmax_fluo",drawcut);
1676    histo = (TH1*)gPad->GetPrimitive("h_Nmax_Hmax_fluo");
1677    histo->SetTitle("Fluo max vs/ Hmax -- energy normalized");
1678    histo->SetYTitle("Fluo max generated km^{-1}");
1679    histo->SetXTitle("Hmax fluo (km)");
1680    histo->SetStats(0);
1681   
1682    pad2->cd();
1683    treepainter->DrawHistos("Nph2_Theta_cer",drawcut);
1684    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Theta_cer");
1685    histo->SetTitle("HALF Cerenkov photons vs/ #theta -- energy normalized");
1686    histo->SetYTitle("HALF cerenkov photons generated");
1687    histo->SetXTitle("#theta (deg)");
1688    histo->SetStats(0);
1689   
1690    pad3->cd();
1691    treepainter->DrawHistos("Nph2_Hmax_fluo",drawcut);
1692    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Hmax_fluo");
1693    histo->SetTitle("Fluo HALF vs/ Hmax -- energy normalized");
1694    histo->SetYTitle("Fluo HALF generated");
1695    histo->SetXTitle("Hmax (km)");
1696    histo->SetStats(0);
1697   
1698    pad4->cd();
1699    treepainter->DrawHistos("Nph2_Hmax_cer",drawcut);
1700    histo = (TH1*)gPad->GetPrimitive("h_Nph2_Hmax_cer");
1701    histo->SetTitle("HALF Cerenkov photons vs/ Hmax -- energy normalized");
1702    histo->SetYTitle("HALF cerenkov photons generated");
1703    histo->SetXTitle("Hmax (km)");
1704    histo->SetStats(0);
1705
1706}
1707
1708//__________________________________________________________________________________________________________________________-
1709void atmo_energy_correl() {
1710    Double_t min = 0;
1711    Double_t max = 0;
1712    Double_t textpos = 0;
1713    Double_t par[3];
1714    TLatex *t = 0;
1715    Int_t nbin = 0;
1716    TF1* fitfunc = 0;
1717    TF1* mypol = 0;
1718    if(!pad_atmo_photons) pad_atmo_photons = new EGViewer();
1719   
1720
1721
1722    Bool_t flag = true;
1723    if(!chain&&!etree&&!stree) flag = checktree();     
1724    if(!flag) {
1725        Warning("treeatmoEnergycorrel()","No Tree found");
1726        return;
1727    }
1728   
1729    TCanvas *c = pad_atmo_photons->FindCanvas("c_treeatmo_correlenergy");
1730    if(!c ) c = pad_atmo_photons->AddTab("c_treeatmo_correlenergy","c_treeatmo_correlenergy");
1731    TCanvas *c2 = pad_atmo_photons->FindCanvas("c_treeatmo_correlenergy_F");
1732    if(!c2 ) c2 = pad_atmo_photons->AddTab("c_treeatmo_correlenergy_F","c_treeatmo_correlenergy_F");
1733
1734    c->Clear();
1735    c->cd();
1736
1737    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1738    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1739    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1740    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1741
1742    pad1->Draw();
1743    pad2->Draw();
1744    pad3->Draw();
1745    pad4->Draw();
1746
1747
1748    pad1->cd();
1749    treepainter->DrawHistos("Nph2_Energy_fluo",drawcut);
1750    prof_Nphfluo_Energy = (TProfile*)gPad->GetPrimitive("h_Nph2_Energy_fluo");
1751    prof_Nphfluo_Energy->SetTitle("HALF Fluo photons vs/ Energy -- cos#theta normalized");
1752    prof_Nphfluo_Energy->SetYTitle("(log10) HALF fluo photons generated");
1753    prof_Nphfluo_Energy->SetXTitle("log10(Energy) (MeV)");
1754    prof_Nphfluo_Energy->SetStats(0);
1755    min = 0;
1756    max = 0;
1757    for(Int_t m=0; m<prof_Nphfluo_Energy->GetNbinsX(); m++) {
1758        if(prof_Nphfluo_Energy->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nphfluo_Energy->GetBinLowEdge(m+1);
1759        if(prof_Nphfluo_Energy->GetBinContent(prof_Nphfluo_Energy->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nphfluo_Energy->GetBinLowEdge(prof_Nphfluo_Energy->GetNbinsX() - m);
1760        if(fabs(min*max) > 0) break;
1761    }
1762    mypol = new TF1("mypol_Nph2_Energy_fluo",my_linear,min,max,2);
1763    mypol->SetLineWidth(1);
1764    mypol->SetLineColor(2);
1765    mypol->SetParameter(1.,1.);
1766    mypol->SetParLimits(0,-10.,10.);
1767    mypol->SetParLimits(1,1.,1.);
1768    textpos = min + 0.1*(max-min);
1769    prof_Nphfluo_Energy->Fit("mypol_Nph2_Energy_fluo","RCOQ");
1770    mypol->GetParameters(par);
1771    t = new TLatex();
1772    char txt[200];
1773    t->SetTextFont(62); 
1774    t->SetTextSize(0.06);
1775    t->SetTextColor(2);
1776    t->SetTextColor(1);
1777    sprintf(txt,"Nfluo_{1/2} = %.3g  *  (E/MeV)",pow(10,par[0]));
1778    t->DrawLatex(textpos,prof_Nphfluo_Energy->GetMaximum()*0.54,txt);
1779   
1780   
1781    pad2->cd();
1782    treepainter->DrawHistos("Nph2_Energy_cer",drawcut);
1783    prof_NphcerEnergy = (TProfile*)gPad->GetPrimitive("h_Nph2_Energy_cer");
1784    prof_NphcerEnergy->SetTitle("HALF Cerenkov photons vs/ Energy");
1785    prof_NphcerEnergy->SetYTitle("(log10) HALF cerenkov photons generated");
1786    prof_NphcerEnergy->SetXTitle("log10(Energy) (MeV)");
1787    prof_NphcerEnergy->SetStats(0);
1788    min = 0;
1789    max = 0;
1790    for(Int_t m=0; m<prof_NphcerEnergy->GetNbinsX(); m++) {
1791        if(prof_NphcerEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_NphcerEnergy->GetBinLowEdge(m+1);
1792        if(prof_NphcerEnergy->GetBinContent(prof_NphcerEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_NphcerEnergy->GetBinLowEdge(prof_NphcerEnergy->GetNbinsX() - m);
1793        if(fabs(min*max) > 0) break;
1794    }
1795    mypol = new TF1("mypol_Nph2_Energy_cer",my_linear,min,max,2);  // pente fixee a 1
1796    mypol->SetLineWidth(1);
1797    mypol->SetLineColor(2);
1798    mypol->SetParameter(1.,1.);
1799    mypol->SetParLimits(0,-10.,10.);
1800    mypol->SetParLimits(1,1.,1.);
1801    textpos = min + 0.1*(max-min);
1802    prof_NphcerEnergy->Fit("mypol_Nph2_Energy_cer","WRCOQ");
1803    mypol->GetParameters(par);
1804    t = new TLatex();
1805    char txt[200];
1806    t->SetTextFont(62); 
1807    t->SetTextSize(0.06);
1808    t->SetTextColor(2);
1809    t->SetTextColor(1);
1810    sprintf(txt,"Nckov_{1/2} = %.3g  *  (E/MeV)",pow(10,par[0]));
1811    t->DrawLatex(textpos,prof_NphcerEnergy->GetMaximum()*0.54,txt);
1812
1813
1814    pad3->cd();
1815    treepainter->DrawHistos("Nmax_Energy_fluo",drawcut);
1816    prof_NmaxfluoEnergy = (TProfile*)gPad->GetPrimitive("h_Nmax_Energy_fluo");
1817    prof_NmaxfluoEnergy->SetTitle("Fluorescence maximum vs/ Energy -- cos#theta normalized");
1818    prof_NmaxfluoEnergy->SetYTitle("(log10) Fluo max");
1819    prof_NmaxfluoEnergy->SetXTitle("log10(Energy) (MeV)");
1820    prof_NmaxfluoEnergy->SetStats(0);
1821    min = 0;
1822    max = 0;
1823    for(Int_t m=0; m<prof_NmaxfluoEnergy->GetNbinsX(); m++) {
1824        if(prof_NmaxfluoEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_NmaxfluoEnergy->GetBinLowEdge(m+1);
1825        if(prof_NmaxfluoEnergy->GetBinContent(prof_NmaxfluoEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_NmaxfluoEnergy->GetBinLowEdge(prof_NmaxfluoEnergy->GetNbinsX() - m);
1826        if(fabs(min*max) > 0) break;
1827    }
1828    mypol = new TF1("mypol_Nmax_Energy_fluo",my_linear,min,max,2);
1829    mypol->SetLineWidth(1);
1830    mypol->SetLineColor(2);
1831    mypol->SetParameter(1.,1.);
1832    mypol->SetParLimits(0,-10.,10.);
1833    mypol->SetParLimits(1,1.,1.);
1834    textpos = min + 0.1*(max-min);
1835    prof_NmaxfluoEnergy->Fit("mypol_Nmax_Energy_fluo","WRCOQ");
1836    mypol->GetParameters(par);
1837    t = new TLatex();
1838    char txt[200];
1839    t->SetTextFont(62); 
1840    t->SetTextSize(0.06);
1841    t->SetTextColor(2);
1842    t->SetTextColor(1);
1843    sprintf(txt,"Nfluo_{max} = %.3g  *  (E/MeV) cm^{2}/g",pow(10,par[0]));
1844    t->DrawLatex(textpos,prof_NmaxfluoEnergy->GetMaximum()*0.54,txt);
1845
1846    /*
1847    pad4->cd();
1848    treepainter->DrawHistos("Nmax_Energy_cer",drawcut);
1849    prof_NmaxcerEnergy = (TProfile*)gPad->GetPrimitive("h_Nmax_Energy_cer");
1850    prof_NmaxcerEnergy->SetTitle("Cerenkov maximum vs/ Energy");
1851    prof_NmaxcerEnergy->SetYTitle("(log10) Cerenkov max");
1852    prof_NmaxcerEnergy->SetXTitle("log10(Energy) (MeV)");
1853    prof_NmaxcerEnergy->SetStats(0);
1854    min = 0;
1855    max = 0;
1856    for(Int_t m=0; m<prof_NmaxcerEnergy->GetNbinsX(); m++) {
1857        if(prof_NmaxcerEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_NmaxcerEnergy->GetBinLowEdge(m+1);
1858        if(prof_NmaxcerEnergy->GetBinContent(prof_NmaxcerEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_NmaxcerEnergy->GetBinLowEdge(prof_NmaxcerEnergy->GetNbinsX() - m);
1859        if(fabs(min*max) > 0) break;
1860    }
1861    mypol = new TF1("mypol_Nmax_Energy_cer",my_linear,min,max,2);  // pente fixee a 1
1862    mypol->SetLineWidth(1);
1863    mypol->SetLineColor(2);
1864    mypol->SetParameter(1.,1.);
1865    mypol->SetParLimits(0,-10.,10.);
1866    mypol->SetParLimits(1,1.,1.);
1867    textpos = min + 0.1*(max-min);
1868    prof_NmaxcerEnergy->Fit("mypol_Nmax_Energy_cer","WRCOQ");
1869    mypol->GetParameters(par);
1870    t = new TLatex();
1871    char txt[200];
1872    t->SetTextFont(62); 
1873    t->SetTextSize(0.06);
1874    t->SetTextColor(2);
1875    t->SetTextColor(1);
1876    //sprintf(txt,"Nckov_{max} = %.3g  *  (E/MeV)^{%.3g} cm^{2}/g",pow(10,par[0]),par[1]);
1877    sprintf(txt,"Nckov_{max} = %.3g  *  (E/MeV) cm^{2}/g",pow(10,par[0]));
1878    t->DrawLatex(textpos,prof_NmaxcerEnergy->GetMaximum()*0.54,txt);
1879    */
1880
1881   
1882
1883    // *************** PLOTS FOR SIGMA  *********************************************
1884   
1885    c2->Clear();
1886    c2->cd();
1887    TPad *pad12 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
1888    TPad *pad22 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
1889    TPad *pad32 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
1890    TPad *pad42 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
1891
1892    pad12->Draw();
1893    pad22->Draw();
1894    pad32->Draw();
1895    pad42->Draw();
1896    TProfile* prof = 0;
1897
1898    pad12->cd();
1899    treepainter->DrawHistos("Nph2_Energy_fluo_linear",drawcut);
1900    prof = (TProfile*) gDirectory->Get("h_Nph2_Energy_fluo_linear");
1901    nbin = prof->GetNbinsX();
1902    min = prof->GetBinLowEdge(1);
1903    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1904    resolNphfluoEnergy = new TProfile("resolNphfluoEnergy","#splitline{Normalized sigma of Fluo first HALF}{VS/ Energy}",nbin,min,max,1e-6,1e6);
1905    for(Int_t i=0; i<nbin; i++)
1906        if(prof->GetBinContent(i+1)) resolNphfluoEnergy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1907    resolNphfluoEnergy->SetMarkerStyle(27);
1908    resolNphfluoEnergy->SetXTitle("log10(Energy) (MeV)");
1909    resolNphfluoEnergy->Draw();
1910    gPad->SetTopMargin(0.18);
1911   
1912    pad22->cd();
1913    treepainter->DrawHistos("Nph2_Energy_cer_linear",drawcut);
1914    prof = (TProfile*) gDirectory->Get("h_Nph2_Energy_cer_linear");
1915    nbin = prof->GetNbinsX();
1916    min = prof->GetBinLowEdge(1);
1917    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1918    resolNphcerEnergy = new TProfile("resolNphcerEnergy","#splitline{Normalized sigma of Cerenkov first HALF}{VS/ Energy}",nbin,min,max,1e-6,1e6);
1919    for(Int_t i=0; i<nbin; i++) {
1920        if(prof->GetBinContent(i+1)) resolNphcerEnergy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1921    }
1922    resolNphcerEnergy->SetMarkerStyle(27);
1923    resolNphcerEnergy->SetXTitle("log10(Energy) (MeV)");
1924    resolNphcerEnergy->Draw();
1925    gPad->SetTopMargin(0.18);
1926
1927    pad32->cd();
1928    treepainter->DrawHistos("Nmax_Energy_fluo_linear",drawcut);
1929    prof = (TProfile*) gDirectory->Get("h_Nmax_Energy_fluo_linear");
1930    nbin = prof->GetNbinsX();
1931    min = prof->GetBinLowEdge(1);
1932    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1933    resolNmaxfluoEnergy = new TProfile("resolNmaxfluoEnergy","Normalized sigma of Fluo max VS Energy",nbin,min,max,1e-6,1e6);
1934    for(Int_t i=0; i<nbin; i++) {
1935        if(prof->GetBinContent(i+1)) resolNmaxfluoEnergy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1936    }
1937    resolNmaxfluoEnergy->SetMarkerStyle(27);
1938    resolNmaxfluoEnergy->SetXTitle("log10(Energy) (MeV)");
1939    resolNmaxfluoEnergy->Draw();
1940     
1941    /*
1942    pad42->cd();
1943    treepainter->DrawHistos("Nmax_Energy_cer_linear",drawcut);
1944    prof = (TProfile*) gDirectory->Get("h_Nmax_Energy_cer_linear");
1945    nbin = prof->GetNbinsX();
1946    min = prof->GetBinLowEdge(1);
1947    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
1948    resolNmaxcerEnergy = new TProfile("resolNmaxcerEnergy","Normalized sigma of Cerenkov max VS Energy",nbin,min,max,1e-6,1e6);
1949    for(Int_t i=0; i<nbin; i++) {
1950        if(prof->GetBinContent(i+1)) resolNmaxcerEnergy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
1951    }
1952    resolNmaxcerEnergy->SetMarkerStyle(27);
1953    resolNmaxcerEnergy->SetXTitle("log10(Energy) (MeV)");
1954    resolNmaxcerEnergy->Draw();
1955    */
1956}
1957
1958//__________________________________________________________________________________________________________________________-
1959void treepupilhistos() {
1960    Bool_t flag = true;
1961    if(!chain&&!etree&&!stree) flag = checktree();     
1962    if(!flag) {
1963        delete c; c= 0;
1964        Warning("treepupilhistos()","No Tree found");
1965        return;
1966    }
1967    if(!pad_p_on_pupil) pad_p_on_pupil = new EGViewer();
1968   
1969    TCanvas *c = pad_p_on_pupil->FindCanvas("c_treepupil_histo");
1970    if(!c ) c = pad_p_on_pupil->AddTab("c_treepupil_histo","c_treepupil_histo");
1971    TCanvas *c2 = pad_p_on_pupil->FindCanvas("c_treepupil_histo2");
1972    if(!c2 ) c2 = pad_p_on_pupil->AddTab("c_treepupil_histo2","c_treepupil_histo2");
1973
1974    c->Clear();
1975    c->cd();
1976    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.67,0.49,0.99);
1977    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.67,0.99,0.99);
1978    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.34,0.49,0.66);
1979    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.34,0.99,0.66);
1980    TPad *pad5 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.33);
1981    TPad *pad6 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.33);
1982
1983    pad1->Draw();
1984    pad2->Draw();
1985    pad3->Draw();
1986    pad4->Draw();
1987    pad5->Draw();
1988    pad6->Draw();
1989
1990
1991    pad1->cd();
1992    pad1->SetLogy(1);
1993    treepainter->DrawHistos("Nph_CB_P",drawcut);
1994    histo = (TH1*)gPad->GetPrimitive("h_Nph_CB_P");
1995    histo->SetTitle("Backscattered cerenkov on pupil");
1996    histo->SetXTitle("Nph_CB_P");
1997    pad2->cd();
1998    pad2->SetLogy(1);
1999    treepainter->DrawHistos("Nph_CR_P",drawcut);
2000    histo = (TH1*)gPad->GetPrimitive("h_Nph_CR_P");
2001    histo->SetTitle("Reflected cerenkov on pupil");
2002    histo->SetXTitle("Nph_CR_P");
2003    pad3->cd();
2004    pad3->SetLogy(1);
2005    treepainter->DrawHistos("Nph_Fluo_P",drawcut);
2006    histo = (TH1*)gPad->GetPrimitive("h_Nph_Fluo_P");
2007    histo->SetTitle("Fluorescence on pupil");
2008    histo->SetXTitle("Nph_Fluo_P");
2009    pad4->cd();
2010    pad4->SetLogy(1);
2011    treepainter->DrawHistos("Nph_tot_P",drawcut);
2012    histo = (TH1*)gPad->GetPrimitive("h_Nph_tot_P");
2013    histo->SetTitle("Total nb of photons on pupil");
2014    histo->SetXTitle("Nph_tot_P");
2015    pad5->cd();
2016    pad5->SetLogy(1);
2017    treepainter->DrawHistos("Nmax_Fluo_P",drawcut);
2018    histo = (TH1*)gPad->GetPrimitive("h_Nmax_Fluo_P");
2019    histo->SetTitle("Fluorescence maximum on pupil");
2020    histo->SetXTitle("Nmax_Fluo_P");
2021    pad6->cd();
2022    pad6->SetLogy(1);
2023    treepainter->DrawHistos("Sig_Fluo_P",drawcut);
2024    histo = (TH1*)gPad->GetPrimitive("h_Sig_Fluo_P");
2025    histo->SetTitle("Half Width of fluo time profile on pupil");
2026    histo->SetXTitle("Half width (GTU)");
2027
2028
2029//////////////////////////
2030    c2->Clear();
2031    c2->cd();
2032    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2033    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2034    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2035    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2036
2037    pad11->Draw();
2038    pad21->Draw();
2039    pad31->Draw();
2040    pad41->Draw();
2041
2042    pad11->cd();
2043    treepainter->DrawHistos("Wlmean_f_p",drawcut);
2044    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_f_p");
2045    histo->SetTitle("Mean wavelength for fluo on pupil");
2046    histo->SetXTitle("Wlmean_f_p");
2047    pad21->cd();
2048    treepainter->DrawHistos("Wlmean_cb_p",drawcut);
2049    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_cb_p");
2050    histo->SetTitle("Mean wavelength for backscattered cerenkov on pupil");
2051    histo->SetXTitle("Wlmean_cb_p");
2052    pad31->cd();
2053    treepainter->DrawHistos("Wlmean_cr_p",drawcut);
2054    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_cr_p");
2055    histo->SetTitle("Mean wavelength for reflected cerenkov on pupil");
2056    histo->SetXTitle("Wlmean_cr_p");
2057    pad41->cd();
2058    treepainter->DrawHistos("Wlmean_tot_p",drawcut);
2059    histo = (TH1*)gPad->GetPrimitive("h_Wlmean_tot_p");
2060    histo->SetTitle("Mean wavelength of photons on pupil");
2061    histo->SetXTitle("Wlmean_tot_p");
2062}
2063
2064//__________________________________________________________________________________________________________________________-
2065void pupil_theta_correl() {
2066    Double_t min = 0;
2067    Double_t max = 0;
2068    Double_t textpos = 0;
2069    Double_t par[3];
2070    TLatex *t = 0;
2071    Int_t nbin = 0;
2072    TF1* fitfunc = 0;
2073    Bool_t flag = true;
2074    if(!chain&&!etree&&!stree) flag = checktree();     
2075    if(!flag) {
2076        Warning("treepupilEnergycorrel()","No Tree found");
2077        return;
2078    }
2079   
2080    TCanvas *c2 = pad_p_on_pupil->FindCanvas("c_treepupil_correltheta");
2081    if(!c2 ) c2 = pad_p_on_pupil->AddTab("c_treepupil_correltheta","c_treepupil_correltheta");   
2082   
2083    c2->Clear();
2084    c2->cd();
2085
2086    TPad *pad11 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2087    TPad *pad21 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2088    TPad *pad31 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2089    TPad *pad41 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2090
2091    pad11->Draw();
2092    pad21->Draw();
2093    pad31->Draw();
2094    pad41->Draw();
2095
2096    pad11->cd();
2097    treepainter->DrawHistos("Nph_Theta_CB_P",drawcut);
2098    histo = (TH1*)gPad->GetPrimitive("h_Nph_Theta_CB_P");
2099    histo->SetTitle("#splitline{Backscattered Cerenkov on pupil}{VS/ #theta -- energy normalized}");
2100    histo->SetXTitle("#theta (deg)");
2101    histo->SetYTitle("Backscattered Cerenkov on pupil");
2102    histo->SetStats(0);
2103    gPad->SetTopMargin(0.2);
2104   
2105    pad21->cd();
2106    treepainter->DrawHistos("Nph_Theta_CR_P",drawcut);
2107    histo = (TH1*)gPad->GetPrimitive("h_Nph_Theta_CR_P");
2108    histo->SetTitle("#splitline{Reflected Cerenkov on pupil}{VS/ #theta -- energy normalized}");
2109    histo->SetXTitle("#theta (deg)");
2110    histo->SetYTitle("Reflected Cerenkov on pupil");
2111    histo->SetStats(0);
2112    gPad->SetTopMargin(0.2);
2113   
2114    pad31->cd();
2115    treepainter->DrawHistos("Nph_Theta_Fluo2_P",drawcut);
2116    histo = (TH1*)gPad->GetPrimitive("h_Nph_Theta_Fluo2_P");
2117    histo->SetTitle("HALF Fluo on pupil vs/ #theta -- energy normalized");
2118    histo->SetXTitle("#theta (deg)");
2119    histo->SetYTitle("HALF Fluorescence on pupil");
2120    histo->SetStats(0);
2121    min = 0;
2122    max = 0;
2123    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
2124        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
2125        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
2126        if(fabs(min*max) > 0) break;
2127    }
2128    /*
2129    mypol = new TF1("mypol_Nph_Theta_Fluo2_P",my_overcos20,min,max,1);
2130    mypol->SetLineWidth(1);
2131    mypol->SetLineColor(2);
2132    mypol->SetParameters(1.e-11,0);
2133    mypol->SetParLimits(0,0.1e-11,100.e-11);
2134    textpos = min + 0.1*(max-min);
2135    histo->Fit("mypol_Nph_Theta_Fluo2_P","WROQ");
2136    mypol->GetParameters(par);
2137    t = new TLatex();
2138    char txt[200];
2139    t->SetTextFont(62); 
2140    t->SetTextSize(0.06);
2141    t->SetTextColor(2);
2142    t->SetTextColor(1);
2143    sprintf(txt,"%.3g EeV^{-1}/ cos(#theta)",par[0]*1e12);
2144    t->DrawLatex(textpos,histo->GetMaximum()*0.74,txt);
2145    */
2146    gPad->SetTopMargin(0.15);
2147   
2148    pad41->cd();
2149    treepainter->DrawHistos("Nph_airscat_effect",drawcut);
2150    histo = (TH1*)gPad->GetPrimitive("h_Nph_airscat_effect");
2151    histo->SetTitle("Increase of signal due to air scattered ckov");
2152    histo->SetYTitle("total / (fluo+reflected ckov)");
2153    histo->SetXTitle("#theta (deg)");
2154   
2155/////////////////////////////////////   
2156    TCanvas *c = pad_p_on_pupil->FindCanvas("c_treepupil_correltheta2");
2157    if(!c ) c = pad_p_on_pupil->AddTab("c_treepupil_correltheta2","c_treepupil_correltheta2");   
2158   
2159    c->Clear();
2160    c->cd();
2161
2162    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2163    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2164    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2165    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2166
2167    pad1->Draw();
2168    pad2->Draw();
2169    pad3->Draw();
2170    pad4->Draw();
2171   
2172    pad1->cd();
2173    treepainter->DrawHistos("Nph_Theta_Fluo_P_crash",drawcut);
2174    histo = (TH1*)gPad->GetPrimitive("h_Nph_Theta_Fluo_P_crash");
2175    histo->SetTitle("Crash effect on fluo on pupil");
2176    histo->SetYTitle("Half/full FLuo on pupil");
2177    histo->SetXTitle("#theta (deg)");
2178   
2179    pad2->cd();
2180    treepainter->DrawHistos("TimeWidth2_Theta",drawcut);
2181    histo = (TH1*)gPad->GetPrimitive("h_TimeWidth2_Theta");
2182    histo->SetTitle("HALF width of fluo time profile on pupil");
2183    histo->SetYTitle("Half time width FLuo on pupil (nb of GTU)");
2184    histo->SetXTitle("#theta (deg)");
2185    min = 0;
2186    max = 0;
2187    for(Int_t m=0; m<histo->GetNbinsX(); m++) {
2188        if(histo->GetBinContent(m+1) && fabs(min) == 0) min = histo->GetBinLowEdge(m+1);
2189        if(histo->GetBinContent(histo->GetNbinsX() - m) && fabs(max) == 0) max = histo->GetBinLowEdge(histo->GetNbinsX() - m);
2190        if(fabs(min*max) > 0) break;
2191    }
2192    /*
2193    mypol = new TF1("mypol_TimeWidth2_Theta",my_overcos20,min,max,1);
2194    mypol->SetLineWidth(1);
2195    mypol->SetLineColor(2);
2196    mypol->SetParameters(1.,0);
2197    mypol->SetParLimits(0,0.1,100.);
2198    textpos = min + 0.1*(max-min);
2199    histo->Fit("mypol_TimeWidth2_Theta","WROQ");
2200    mypol->GetParameters(par);
2201    t = new TLatex();
2202    char txt[200];
2203    t->SetTextFont(62); 
2204    t->SetTextSize(0.06);
2205    t->SetTextColor(2);
2206    t->SetTextColor(1);
2207    sprintf(txt,"%.3g / cos(#theta)  GTU",par[0]);
2208    t->DrawLatex(textpos,histo->GetMaximum()*0.74,txt);
2209    */
2210    gPad->SetTopMargin(0.15);
2211}
2212
2213//__________________________________________________________________________________________________________________________-
2214void pupil_energy_correl() {
2215    Double_t min = 0;
2216    Double_t max = 0;
2217    Double_t textpos = 0;
2218    Double_t par[3];
2219    TLatex *t = 0;
2220    Int_t nbin = 0;
2221    TF1* fitfunc = 0;
2222    Bool_t flag = true;
2223    if(!chain&&!etree&&!stree) flag = checktree();     
2224    if(!flag) {
2225        Warning("treepupilEnergycorrel()","No Tree found");
2226        return;
2227    }
2228   
2229    TCanvas *c = pad_p_on_pupil->FindCanvas("c_treepupil_correlenergy");
2230    if(!c ) c = pad_p_on_pupil->AddTab("c_treepupil_correlenergy","c_treepupil_correlenergy");
2231    TCanvas *cc1 = pad_p_on_pupil->FindCanvas("c_treepupil_correlenergy_F");
2232    if(!cc1 ) cc1 = pad_p_on_pupil->AddTab("c_treepupil_correlenergy_F","c_treepupil_correlenergy_F");
2233
2234    c->Clear();
2235    c->cd();
2236
2237    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2238    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2239    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2240    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2241
2242    pad1->Draw();
2243    pad2->Draw();
2244    pad3->Draw();
2245    pad4->Draw();
2246       
2247    pad1->cd();
2248    treepainter->DrawHistos("Nph2_Energy_Fluo_P",drawcut);
2249    prof_Nph2fluoEnergy = (TProfile*)gPad->GetPrimitive("h_Nph2_Energy_Fluo_P");
2250    prof_Nph2fluoEnergy->SetTitle("(log10) HALF fluo on pupil vs/ Energy -- cos(#theta) normalization");
2251    prof_Nph2fluoEnergy->SetXTitle("Energy (MeV)");
2252    prof_Nph2fluoEnergy->SetYTitle("HALF fluorescence on pupil");
2253    prof_Nph2fluoEnergy->SetStats(0);
2254    min = 0;
2255    max = 0;
2256    for(Int_t m=0; m<prof_Nph2fluoEnergy->GetNbinsX(); m++) {
2257        if(prof_Nph2fluoEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nph2fluoEnergy->GetBinLowEdge(m+1);
2258        if(prof_Nph2fluoEnergy->GetBinContent(prof_Nph2fluoEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nph2fluoEnergy->GetBinLowEdge(prof_Nph2fluoEnergy->GetNbinsX() - m);
2259        if(fabs(min*max) > 0) break;
2260    }
2261    mypol = new TF1("mypol_Nph2_Energy_Fluo_P",my_linear,min,max,2);  // pente fixee a 1
2262    mypol->SetLineWidth(1);
2263    mypol->SetLineColor(2);
2264    mypol->SetParameter(1.,1.);
2265    mypol->SetParLimits(0,-20.,10.);
2266    mypol->SetParLimits(1,1.,1.);
2267    textpos = min + 0.1*(max-min);
2268    prof_Nph2fluoEnergy->Fit("mypol_Nph2_Energy_Fluo_P","WRCOQ");
2269    mypol->GetParameters(par);
2270    t = new TLatex();
2271    char txt[200];
2272    t->SetTextFont(62); 
2273    t->SetTextSize(0.06);
2274    t->SetTextColor(2);
2275    t->SetTextColor(1);
2276    sprintf(txt,"Nfluo_{1/2} = %.2g  *  (E/EeV)",pow(10,par[0])*1e12);
2277    t->DrawLatex(textpos,prof_Nph2fluoEnergy->GetMaximum()*0.44,txt);
2278   
2279    pad2->cd();
2280    treepainter->DrawHistos("Nph_Energy_CB_P",drawcut);
2281    prof_NphbackEnergy = (TProfile*)gPad->GetPrimitive("h_Nph_Energy_CB_P");
2282    prof_NphbackEnergy->SetTitle("(log10) Backscattered Cerenkov on pupil vs/ Energy -- cos(#theta) normalized");
2283    prof_NphbackEnergy->SetXTitle("Energy (MeV)");
2284    prof_NphbackEnergy->SetYTitle("Backscattered cerenkov on pupil");
2285    prof_NphbackEnergy->SetStats(0);
2286    min = 0;
2287    max = 0;
2288    for(Int_t m=0; m<prof_NphbackEnergy->GetNbinsX(); m++) {
2289        if(prof_NphbackEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_NphbackEnergy->GetBinLowEdge(m+1);
2290        if(prof_NphbackEnergy->GetBinContent(prof_NphbackEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_NphbackEnergy->GetBinLowEdge(prof_NphbackEnergy->GetNbinsX() - m);
2291        if(fabs(min*max) > 0) break;
2292    }
2293    mypol = new TF1("mypol_Nph_Energy_CB_P",my_linear,min,max,2);
2294    mypol->SetLineWidth(1);
2295    mypol->SetLineColor(2);
2296    mypol->SetParameter(1.,1.);
2297    mypol->SetParLimits(0,-20.,10);
2298    mypol->SetParLimits(1,1.,1.);
2299    textpos = min + 0.1*(max-min);
2300    prof_NphbackEnergy->Fit("mypol_Nph_Energy_CB_P","WRCOQ");
2301    mypol->GetParameters(par);
2302    t = new TLatex();
2303    char txt[200];
2304    t->SetTextFont(62); 
2305    t->SetTextSize(0.06);
2306    t->SetTextColor(2);
2307    t->SetTextColor(1);
2308    sprintf(txt,"Nckov_back = %.2g  *  (E/EeV)",pow(10,par[0])*1e12);
2309    t->DrawLatex(textpos,prof_NphbackEnergy->GetMaximum()*0.44,txt);
2310
2311    pad3->cd();
2312    treepainter->DrawHistos("Nmax_Energy_Fluo_P",drawcut);
2313    prof_Nmaxfluo_pupil_Energy = (TProfile*)gPad->GetPrimitive("h_Nmax_Energy_Fluo_P");
2314    prof_Nmaxfluo_pupil_Energy->SetTitle("(log10) Fluo max on pupil vs/ Energy -- cos(#theta) normalized");
2315    prof_Nmaxfluo_pupil_Energy->SetXTitle("Energy (MeV)");
2316    prof_Nmaxfluo_pupil_Energy->SetYTitle("Fluo max on pupil");
2317    prof_Nmaxfluo_pupil_Energy->SetStats(0);
2318    min = 0;
2319    max = 0;
2320    for(Int_t m=0; m<prof_Nmaxfluo_pupil_Energy->GetNbinsX(); m++) {
2321        if(prof_Nmaxfluo_pupil_Energy->GetBinContent(m+1) && fabs(min) == 0) min = prof_Nmaxfluo_pupil_Energy->GetBinLowEdge(m+1);
2322        if(prof_Nmaxfluo_pupil_Energy->GetBinContent(prof_Nmaxfluo_pupil_Energy->GetNbinsX() - m) && fabs(max) == 0) max = prof_Nmaxfluo_pupil_Energy->GetBinLowEdge(prof_Nmaxfluo_pupil_Energy->GetNbinsX() - m);
2323        if(fabs(min*max) > 0) break;
2324    }
2325    mypol = new TF1("mypol_Nmax_Energy_Fluo_P",my_linear,min,max,2);
2326    mypol->SetLineWidth(1);
2327    mypol->SetLineColor(2);
2328    mypol->SetParameter(1.,1.);
2329    mypol->SetParLimits(0,-20.,10.);
2330    mypol->SetParLimits(1,1.,1.);
2331    textpos = min + 0.1*(max-min);
2332    prof_Nmaxfluo_pupil_Energy->Fit("mypol_Nmax_Energy_Fluo_P","WRCOQ");
2333    mypol->GetParameters(par);
2334    t = new TLatex();
2335    char txt[200];
2336    t->SetTextFont(62); 
2337    t->SetTextSize(0.06);
2338    t->SetTextColor(2);
2339    t->SetTextColor(1);
2340    sprintf(txt,"Nfluo_{max} = %.2g  *  (E/EeV) GTU^{-1}",pow(10,par[0])*1e12);
2341    t->DrawLatex(textpos,prof_Nmaxfluo_pupil_Energy->GetMaximum()*0.30,txt);
2342   
2343    pad4->cd();
2344    treepainter->DrawHistos("Nph_Energy_CR_P",drawcut);
2345    prof_NphrefEnergy = (TProfile*)gPad->GetPrimitive("h_Nph_Energy_CR_P");
2346    prof_NphrefEnergy->SetYTitle("(log10) Reflected cerenkov on pupil");
2347    prof_NphrefEnergy->SetXTitle("Energy (MeV)");
2348    prof_NphrefEnergy->SetTitle("Reflected Cerenkov on pupil vs/ Energy");
2349    prof_NphrefEnergy->SetStats(0);
2350    min = 0;
2351    max = 0;
2352    for(Int_t m=0; m<prof_NphrefEnergy->GetNbinsX(); m++) {
2353        if(prof_NphrefEnergy->GetBinContent(m+1) && fabs(min) == 0) min = prof_NphrefEnergy->GetBinLowEdge(m+1);
2354        if(prof_NphrefEnergy->GetBinContent(prof_NphrefEnergy->GetNbinsX() - m) && fabs(max) == 0) max = prof_NphrefEnergy->GetBinLowEdge(prof_NphrefEnergy->GetNbinsX() - m);
2355        if(fabs(min*max) > 0) break;
2356    }
2357    mypol = new TF1("mypol_Nph_Energy_CR_P",my_linear,min,max,2);  // pente fixee a 1
2358    mypol->SetLineWidth(1);
2359    mypol->SetLineColor(2);
2360    mypol->SetParameter(1.,1.);
2361    mypol->SetParLimits(0,-20.,10.);
2362    mypol->SetParLimits(1,1.,1.);
2363
2364    textpos = min + 0.07*(max-min);
2365    prof_NphrefEnergy->Fit("mypol_Nph_Energy_CR_P","WRCOQ");
2366    mypol->GetParameters(par);
2367    t = new TLatex();
2368    char txt[200];
2369    t->SetTextFont(62); 
2370    t->SetTextSize(0.06);
2371    t->SetTextColor(2);
2372    t->SetTextColor(1);
2373    sprintf(txt,"Nckov_reflected = %.2g  *  (E/EeV)",pow(10,par[0])*1e12);
2374    t->DrawLatex(textpos,prof_NphrefEnergy->GetMaximum()*0.34,txt);
2375   
2376   
2377    // *************** PLOTS FOR SIGMA  ********************************************************
2378   
2379    cc1->Clear();
2380    cc1->cd();
2381    TPad *pad111 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2382    TPad *pad112 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2383    TPad *pad113 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2384    TPad *pad114 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2385
2386    pad111->Draw();
2387    pad112->Draw();
2388    pad113->Draw();
2389    pad114->Draw();
2390    TProfile* prof = 0;
2391   
2392
2393    pad111->cd();
2394    treepainter->DrawHistos("Nph2_Energy_Fluo_P_linear",drawcut);
2395    prof = (TProfile*)gDirectory->Get("h_Nph2_Energy_Fluo_P_linear");
2396    nbin = prof->GetNbinsX();
2397    min = prof->GetBinLowEdge(1);
2398    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
2399    resolNph2_pupil_fluoEnergy = new TProfile("resolNph2_pupil_fluoEnergy","#splitline{Normalized #sigma of Fluo first HALF}{VS/ Energy  (on pupil) -- cos(#theta) normalized}",nbin,min,max,1e-6,1e6);
2400    for(Int_t i=0; i<nbin; i++)
2401        if(prof->GetBinContent(i+1)) resolNph2_pupil_fluoEnergy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
2402    resolNph2_pupil_fluoEnergy->SetMarkerStyle(27);
2403    resolNph2_pupil_fluoEnergy->SetXTitle("log10(Energy) (MeV)");
2404    resolNph2_pupil_fluoEnergy->Draw();
2405    gPad->SetTopMargin(0.18);
2406   
2407    pad112->cd();
2408    treepainter->DrawHistos("Nph_Energy_CB_P_linear",drawcut);
2409    prof = (TProfile*)gDirectory->Get("h_Nph_Energy_CB_P_linear");
2410    nbin = prof->GetNbinsX();
2411    min = prof->GetBinLowEdge(1);
2412    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
2413    resol_Nphcerback_Energy = new TProfile("resol_Nphcerback_Energy","#splitline{Normalized #sigma of backscattered ckov}{VS/ Energy   (on pupil) -- cos(#theta) normalized}",nbin,min,max,1e-6,1e6);
2414    for(Int_t i=0; i<nbin; i++) {
2415        if(prof->GetBinContent(i+1)) resol_Nphcerback_Energy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
2416    }
2417    resol_Nphcerback_Energy->SetMarkerStyle(27);
2418    resol_Nphcerback_Energy->SetXTitle("log10(Energy) (MeV)");
2419    resol_Nphcerback_Energy->Draw();
2420    gPad->SetTopMargin(0.18);
2421
2422    pad113->cd();
2423    treepainter->DrawHistos("Nmax_Energy_Fluo_P_linear",drawcut);
2424    prof = (TProfile*)gDirectory->Get("h_Nmax_Energy_Fluo_P_linear");
2425    nbin = prof->GetNbinsX();
2426    min = prof->GetBinLowEdge(1);
2427    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
2428    resolNmaxfluo_pupil_Energy = new TProfile("resolNmaxfluo_pupil_Energy","Normalized #sigma of Fluo max on pupil VS Energy",nbin,min,max,1e-6,1e6);
2429    for(Int_t i=0; i<nbin; i++) {
2430        if(prof->GetBinContent(i+1)) resolNmaxfluo_pupil_Energy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
2431    }
2432    resolNmaxfluo_pupil_Energy->SetMarkerStyle(27);
2433    resolNmaxfluo_pupil_Energy->SetXTitle("log10(Energy) (MeV)");
2434    resolNmaxfluo_pupil_Energy->Draw();
2435     
2436    pad114->cd();
2437    treepainter->DrawHistos("Nph_Energy_CR_P_linear",drawcut);
2438    prof = (TProfile*)gDirectory->Get("h_Nph_Energy_CR_P_linear");
2439    nbin = prof->GetNbinsX();
2440    min = prof->GetBinLowEdge(1);
2441    max =  prof->GetBinLowEdge(nbin) + prof->GetBinWidth(nbin);
2442    resol_Nohcerreflec_Energy = new TProfile("resol_Nohcerreflec_Energy","Normalized #sigma of reflected Ckov VS Energy  (on pupil)",nbin,min,max,1e-6,1e6);
2443    for(Int_t i=0; i<nbin; i++) {
2444        if(prof->GetBinContent(i+1)) resol_Nohcerreflec_Energy->Fill(prof->GetBinCenter(i+1),prof->GetBinError(i+1)/prof->GetBinContent(i+1));
2445    }
2446    resol_Nohcerreflec_Energy->SetMarkerStyle(27);
2447    resol_Nohcerreflec_Energy->SetXTitle("log10(Energy) (MeV)");
2448    resol_Nohcerreflec_Energy->Draw();
2449}
2450
2451//__________________________________________________________________________________________________________________________-
2452void transmission() {
2453    Double_t min = 0;
2454    Double_t max = 0;
2455    Double_t textpos = 0;
2456    Double_t par[3];
2457    TLatex *t = 0;
2458    Int_t nbin = 0;
2459    TF1* fitfunc = 0;
2460    Bool_t flag = true;
2461    if(!chain&&!etree&&!stree) flag = checktree();     
2462    if(!flag) {
2463        Warning("treepupilEnergycorrel()","No Tree found");
2464        return;
2465    }
2466   
2467    TCanvas *c3 = pad_p_on_pupil->FindCanvas("c_treepupil_correl_trans");
2468    if(!c3 ) c3 = pad_p_on_pupil->AddTab("c_treepupil_correl_trans","c_treepupil_correl_trans");
2469
2470    c3->Clear();
2471    c3->cd();
2472    TPad *pad12 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2473    TPad *pad22 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2474    TPad *pad32 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2475    TPad *pad42 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2476
2477    pad12->Draw();
2478    pad22->Draw();
2479    pad32->Draw();
2480    pad42->Draw();
2481   
2482    pad12->cd();
2483    treepainter->DrawHistos("Trans_max_mean",drawcut);
2484    histo = (TH1*)gPad->GetPrimitive("h_Trans_max_mean");
2485    histo->SetTitle("Relation between T_{mean} and T_{fluomax} for FLUO");
2486    histo->SetXTitle("#theta (deg)");
2487    histo->SetYTitle("(T_{fluomax} - T_{mean})/T_{mean}");
2488    histo->SetStats(0);
2489   
2490    pad22->cd();
2491    treepainter->DrawHistos("Trans_GTUmax_max",drawcut);
2492    histo = (TH1*)gPad->GetPrimitive("h_Trans_GTUmax_max");
2493    histo->SetTitle("Relative difference between T_{maxfluo} and T_{GTUmax} for FLUO");
2494    histo->SetXTitle("(T_{GTUmax} - T_{fluomax})/T_{fluomax}");
2495   
2496    pad32->cd();
2497    treepainter->DrawHistos("Trans_ckov_Theta",drawcut);
2498    histo = (TH1*)gPad->GetPrimitive("h_Trans_ckov_Theta");
2499    histo->SetTitle("Cerenkov transmission");
2500    histo->SetXTitle("#theta (deg)");
2501    histo->SetYTitle("Ratio   on_pupil / created");
2502    histo->SetStats(0);
2503   
2504    pad42->cd();
2505    treepainter->DrawHistos("Trans_fluo_Theta",drawcut);
2506    histo = (TH1*)gPad->GetPrimitive("h_Trans_fluo_Theta");
2507    histo->SetTitle("Fluo transmission");
2508    histo->SetXTitle("#theta (deg)");
2509    histo->SetYTitle("Ratio   on_pupil / created");
2510    histo->SetStats(0);
2511}
2512
2513//__________________________________________________________________________________________________________________________-
2514void Miscellaneous() {
2515    Double_t min = 0;
2516    Double_t max = 0;
2517    Double_t textpos = 0;
2518    Double_t par[3];
2519    TLatex *t = 0;
2520    Int_t nbin = 0;
2521    TF1* fitfunc = 0;
2522    Bool_t flag = true;
2523    if(!chain&&!etree&&!stree) flag = checktree();     
2524    if(!flag) {
2525        Warning("treepupilEnergycorrel()","No Tree found");
2526        return;
2527    }
2528
2529    if(!pad_p_on_pupil) pad_p_on_pupil = new EGViewer();
2530   
2531    TCanvas *c = pad_p_on_pupil->FindCanvas("c_treepupil_miscel");
2532    if(!c ) c = pad_p_on_pupil->AddTab("c_treepupil_miscel","c_treepupil_miscel");
2533
2534
2535/////////////////////////
2536    c->Clear();
2537    c->cd();
2538
2539    TPad *pad1 = new TPad("Pad1","Pad1",0.01,0.51,0.49,0.99);
2540    TPad *pad2 = new TPad("Pad2","Pad2",0.51,0.51,0.99,0.99);
2541    TPad *pad3 = new TPad("Pad3","Pad3",0.01,0.01,0.49,0.49);
2542    TPad *pad4 = new TPad("Pad4","Pad4",0.51,0.01,0.99,0.49);
2543
2544    pad1->Draw();
2545    pad2->Draw();
2546    pad3->Draw();
2547    pad4->Draw();
2548   
2549    pad1->cd();
2550    treepainter->DrawHistos("TimeWidth2_Hmax",drawcut);
2551    histo = (TH1*)gPad->GetPrimitive("h_TimeWidth2_Hmax");
2552    histo->SetTitle("Time HALF width on pupil vs/ height of shower maximum");
2553    histo->SetXTitle("Height of fluo max (km)");
2554    histo->SetYTitle("Fluo time half width on pupil (GTU)");
2555    histo->SetStats(0);
2556   
2557    pad2->cd();
2558    treepainter->DrawHistos("Yield_max_mean",drawcut);
2559    histo = (TH1*)gPad->GetPrimitive("h_Yield_max_mean");
2560    histo->SetTitle("Difference between FLUO Yield_{mean} and Yield_{atmax}");
2561    histo->SetXTitle("#theta (deg)");
2562    histo->SetYTitle("(Y_{fluomax} - Y_{mean}) / Y_{mean}");
2563    histo->SetStats(0);
2564   
2565    pad3->cd();
2566    treepainter->DrawHistos("Yield_GTUmax_max",drawcut);
2567    histo = (TH1*)gPad->GetPrimitive("h_Yield_GTUmax_max");
2568    histo->SetTitle("Difference between FLUO Yield_{GTUmax} and Yield_{atmax}");
2569    histo->SetXTitle("#theta (deg)");
2570    histo->SetYTitle("(Y_{GTUmax} - Y_{atmax}) / Y_{atmax}");
2571    histo->SetStats(0);
2572   
2573    pad4->cd();
2574    treepainter->DrawHistos("Xmax_shift",drawcut);
2575    histo = (TH1*)gPad->GetPrimitive("h_Xmax_shift");
2576    histo->SetTitle("Xmax shift due to scattered cerenkov");
2577    histo->SetXTitle("#theta (deg)");
2578    histo->SetYTitle("(Xscat_{GTUmax} - X_{GTUmax}) / X_{GTUmax}");
2579}
2580
2581//__________________________________________________________________________________________________________________________-
2582void simpletreeviewer() {   
2583    stree=(TTree*)gDirectory->Get("stree"); 
2584    if(!stree) {
2585        Warning("simpletreeviewer()","No stree found");
2586        Printf("I do not find simple results tree.");
2587        return;
2588    }
2589    stree->StartViewer();
2590}
2591
2592//__________________________________________________________________________________________________________________________-
2593void savesimpletree() {
2594    if(!stree) {
2595      Printf("Cannot save stree, it has not been generated");
2596      return;
2597    }
2598
2599    TFile *f = 0;
2600    delete f;
2601    Bool_t wantosave = true;
2602   
2603    TGFileInfo* info = new TGFileInfo();
2604    while ( info->fFilename == 0 || !f || f->IsZombie() ) {
2605        TGFileDialog* filedialog = new TGFileDialog(gClient->GetRoot(), 0, kFDSave, info);
2606        if( info->fFilename && !SAVED ) f = new TFile( info->fFilename,"create");
2607        else {wantosave = false; break;}
2608    }
2609    if(wantosave && !SAVED) stree->Write();
2610    else if(SAVED) Info("savesimpletree()","stree is already saved");
2611    SAVED = true;
2612   
2613    stree->Write();
2614    f->Close();
2615}
2616
2617//__________________________________________________________________________________________________________________________-
2618Double_t my_linear(Double_t *x, Double_t *par) {
2619    return par[0] + par[1]*x[0];
2620}
2621
2622//__________________________________________________________________________________________________________________________-
2623Double_t my_expon12(Double_t *x, Double_t *par) {
2624    if(x[0] > 12. || x[0]<0.2) {
2625       TF1::RejectPoint();
2626    }
2627    return exp(par[0] + par[1]*x[0]);
2628}
2629
2630//__________________________________________________________________________________________________________________________-
2631Double_t my_expon15(Double_t *x, Double_t *par) {
2632    if(x[0] > 15. || x[0]<0.2) {
2633       TF1::RejectPoint();
2634    }
2635    return exp(par[0] + par[1]*x[0]);
2636}
2637
2638//__________________________________________________________________________________________________________________________-
2639Double_t my_logcos(Double_t *x, Double_t *par) {
2640    if(cos(x[0]*TMath::Pi()/180.) <= 0) {
2641       TF1::RejectPoint();
2642       return 0;
2643    }
2644    if(x[0] > 60 || x[0]<5) {
2645       TF1::RejectPoint();
2646    }
2647    return par[0] + par[1]*log(cos(x[0]*TMath::Pi()/180.));
2648}
2649
2650//__________________________________________________________________________________________________________________________-
2651Double_t my_overcos60(Double_t *x, Double_t *par) {
2652    if(cos(x[0]*TMath::Pi()/180.) <= 0.) {
2653       TF1::RejectPoint();
2654       return 0;
2655    }
2656    if(x[0] > 60) {
2657       TF1::RejectPoint();
2658    }
2659    return par[0] / cos(x[0]*TMath::Pi()/180.);
2660}
2661
2662//__________________________________________________________________________________________________________________________-
2663Double_t my_overcos20(Double_t *x, Double_t *par) {
2664    if(cos(x[0]*TMath::Pi()/180.) <= 0.) {
2665       TF1::RejectPoint();
2666       return 0;
2667    }
2668    if(x[0] > 20) {
2669       TF1::RejectPoint();
2670    }
2671    return par[0] / cos(x[0]*TMath::Pi()/180.);
2672}
2673
Note: See TracBrowser for help on using the repository browser.