Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ETreePainter.cc @ 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: 82.0 KB

Line
1	// $Id: ETreePainter.cc 2923 2011-06-12 19:30:46Z mabl $
2	// Author: Anne Stutz 2005/03/03
3
4	// $Id: ETreePainter.cc 2923 2011-06-12 19:30:46Z mabl $
5	// Author: Anne Stutz 2005/03/03
6
7	/*****************************************************************************
8	* ESAF: Euso Simulation and Analysis Framework *
9	* *
10	* Id: ETreePainter *
11	* Package: <packagename> *
12	* Coordinator: <coordinator> *
13	* *
14	*****************************************************************************/
15
16	//_____________________________________________________________________________
17	//
18	// ETreePainter
19	//
20	// <extensive class description>
21	//
22	// Config file parameters
23	// ======================
24	//
25	// <parameter name>: <parameter description>
26	// -Valid options: <available options>
27	//
28	//
29	// WARNING : the relations between showerstep ID and corresponding bunchofphotons ID are :
30	// - array[i] for showerstep
31	// - array[2i] for fluo bunch -> bunchID = 2i+1
32	// - array[2i+1] for ckov bunch -> bunchID = 2i+2
33	//
34	// //TOFIX : GTUmax treatment MUST be changed if fluo backscattering
35	//
36	//
37
38	#include <vector>
39	#include "ETreePainter.hh"
40	#include "ESystemOfUnits.hh"
41	#include "EAtmosphere.hh"
42	#include "ETruth.hh"
43	#include "EShower.hh"
44	#include "EShowerStep.hh"
45	#include "EBunchPhotons.hh"
46	#include "EEvent.hh"
47
48	#include <TVirtualGeoTrack.h>
49	#include <TGeoVolume.h>
50	#include <TGeoManager.h>
51	#include <TGeoTube.h>
52	#include <TTree.h>
53	#include <TChain.h>
54	#include <TH1F.h>
55	#include <TF1.h>
56	#include <TNtuple.h>
57	#include <TCut.h>
58	#include <iostream>
59	#include <TRefArray.h>
60	#include <TMath.h>
61	#include <ConfigFileParser.hh>
62	#include "Config.hh"
63
64	ClassImp(ETreePainter)
65
66	extern Double_t Zv(const TVector3&);
67
68	using namespace sou;
69	using namespace std;
70
71	// TOOL FUNCTION
72	//_____________________________________________________________________________
73	Double_t ETreePainter_my_gaussian_fit(Double_t x, Double_t par) {
74	//
75	// par[3] allows to fit a sub-range only
76	//
77
78	if(x[0] > par[3]) {
79	TF1::RejectPoint();
80	return 0;
81	}
82
83	return par[0]exp(-0.5pow((x[0] - par[1]) / par[2],2));
84	}
85
86	//_____________________________________________________________________________
87	ETreePainter::ETreePainter( TTree* etree ) {
88	//
89	// Constructor
90	//
91
92	// if stree already implemented, only Draw() method can be used
93	if(!etree) Fatal("ETreePainter ctor","TTree given in argument is NULL");
94	if(strncmp(etree->GetName(),"stree",5) == 0) {
95	fStree = etree;
96	Info("ETreePainter()","Better to CHECK VALUES used for [altEUSO, pupil radius, GTU length] at Stree building");
97	}
98	else {
99	ConfigFileParser* pConf = Config::Get()->GetCF("Electronics","MacroCell");
100	fGTU = pConf->GetNum("MacroCell.fGtuTimeLength");
101	fGTU/=1000.;
102
103	ConfigFileParser* pConf2 = Config::Get()->GetCF("Electronics","EusoDetector");
104	fRadius = pConf2->GetNum("EusoDetector.fMaxRadius");
105
106	ConfigFileParser* pConf3 = Config::Get()->GetCF("General","Euso");
107	Double_t fHeight = pConf3->GetNum("Euso.fAltitude");
108
109	Info("ETreePainter()","Building Stree");
110	fEUSO = TVector3(0,0,fHeight*km);
111	fRadius = fRadius*mm;
112	fGTU = fGTU*microsecond;
113	fTree = etree;
114	fChain = NULL;
115	fEv = NULL;
116	fStree = NULL;
117	Build();
118	Info("ETreePainter()","BUILT FROM A TREE\nUSE pupil radius = %f mm, euso altitude = %f km, GTU length = %f ms",fRadius/mm,fEUSO.Z()/km,fGTU/microsecond);
119	}
120	}
121
122	//_____________________________________________________________________________
123	ETreePainter::ETreePainter( TChain* chain ) {
124	//
125	// Constructor
126	//
127
128	ConfigFileParser* pConf = Config::Get()->GetCF("Electronics","MacroCell");
129	fGTU = pConf->GetNum("MacroCell.fGtuTimeLength");
130	fGTU/=1000.;
131
132	ConfigFileParser* pConf2 = Config::Get()->GetCF("Electronics","EusoDetector");
133	fRadius = pConf2->GetNum("EusoDetector.fMaxRadius");
134
135	ConfigFileParser* pConf3 = Config::Get()->GetCF("General","Euso");
136	Double_t fHeight = pConf3->GetNum("Euso.fAltitude");
137
138
139
140	Info("ETreePainter()","Building Stree");
141	fEUSO = TVector3(0,0,fHeight*km);
142	fRadius = fRadius*mm;
143	fGTU = fGTU*microsecond;
144	fChain = chain;
145	fTree = chain;
146	fEv = NULL;
147	fStree = NULL;
148	Build();
149	Info("ETreePainter()","BUILT FROM A CHAIN\nUSE pupil radius = %f mm, euso altitude = %f km, GTU length = %f ms",fRadius/mm,fEUSO.Z()/km,fGTU/microsecond);
150	}
151
152	//_____________________________________________________________________________
153	ETreePainter::~ETreePainter() {
154	//
155	// Destructor
156	//
157	}
158
159	//_____________________________________________________________________________
160	void ETreePainter::Build() {
161	//
162	// Build
163	//
164
165	if ( !fTree ) {
166	Info("Build()","ETree object is NULL. Painter made zombie.");
167	MakeZombie();
168	return;
169	}
170
171	fEv = new EEvent();
172
173	if ( fChain ) {
174	SetBranchesStatus(fChain); // MUST be before setbranches()
175	fEv->SetBranches(fChain);
176	}
177	else {
178	SetBranchesStatus(fTree);
179	fEv->SetBranches(fTree);
180	}
181
182	fEvTot = (Int_t)fTree->GetEntries();
183
184	BuildTree();
185
186
187	for (Int_t i=0; i< fEvTot; i++) {
188	Printf( "Entry %d: %d bytes read", i, fTree->GetEntry(i) );
189	fAtmo = fEv->GetAtmosphere();
190	fTruth = fEv->GetTruth();
191	fShower = fEv->GetShower();
192
193	if(fAtmo && fTruth && fShower) FillTrack();
194	else Warning("Build()","EAtmosphere or ETruth or EShower object is NULL");
195	FillTree();
196	}
197	}
198
199	//_____________________________________________________________________________
200	void ETreePainter::SetBranchesStatus(TTree* t) {
201	//
202	// for fChain or fTree
203	//
204	t->SetBranchStatus("*",0);
205	t->SetBranchStatus("fTrue*",1);
206
207	// shower
208	t->SetBranchStatus("fDir*",1);
209	t->SetBranchStatus("fInit*",1);
210	t->SetBranchStatus("fHitGround",1);
211	t->SetBranchStatus("fNumSteps",1);
212	t->SetBranchStatus("fStep*",1);
213
214	// atmosphere
215	t->SetBranchStatus("fNumSingles",1);
216	t->SetBranchStatus("fNumBunches",1);
217	t->SetBranchStatus("fBunch*",1);
218	t->SetBranchStatus("fSingle*",1);
219	}
220
221
222	//_____________________________________________________________________________
223	void ETreePainter::BuildTree() {
224	//
225	// build a simple tree
226	//
227
228	if ( !fStree ) {
229	fStree = new TTree("stree","a Tree with simple datas");
230
231	// flags for checks
232	fStree->Branch("Good_fit", &fTrack.fGood_fit, "fGood_fit/I");
233	fStree->Branch("H_max_last_bin", &fTrack.fH_max_last_bin, "fH_max_last_bin/I");
234	fStree->Branch("X_max_last_bin", &fTrack.fX_max_last_bin, "fX_max_last_bin/I");
235	fStree->Branch("L_max_last_bin", &fTrack.fL_max_last_bin, "fL_max_last_bin/I");
236	fStree->Branch("W_max_last_bin", &fTrack.fW_max_last_bin, "fW_max_last_bin/I");
237	fStree->Branch("H_fluomax_last_bin", &fTrack.fH_fluomax_last_bin, "fH_fluomax_last_bin/I");
238	fStree->Branch("X_fluomax_last_bin", &fTrack.fX_fluomax_last_bin, "fX_fluomax_last_bin/I");
239	fStree->Branch("L_fluomax_last_bin", &fTrack.fL_fluomax_last_bin, "fL_fluomax_last_bin/I");
240	fStree->Branch("W_fluomax_last_bin", &fTrack.fW_fluomax_last_bin, "fW_fluomax_last_bin/I");
241	fStree->Branch("H_ckovmax_last_bin", &fTrack.fH_ckovmax_last_bin, "fH_ckovmax_last_bin/I");
242	fStree->Branch("X_ckovmax_last_bin", &fTrack.fX_ckovmax_last_bin, "fX_ckovmax_last_bin/I");
243
244	// truth
245	fStree->Branch("Energy", &fTrack.fEnergy, "fEnergy/F");
246	fStree->Branch("Theta", &fTrack.fTheta, "fTheta/F");
247	fStree->Branch("Phi", &fTrack.fPhi, "fPhi/F");
248
249	fStree->Branch("H1", &fTrack.fH1, "fH1/F");
250	fStree->Branch("InitPosX", &fTrack.fInitPosX, "fInitPosX/F");
251	fStree->Branch("InitPosY", &fTrack.fInitPosY, "fInitPosY/F");
252	fStree->Branch("InitPosZ", &fTrack.fInitPosZ, "fInitPosZ/F");
253	fStree->Branch("X1", &fTrack.fX1, "fX1/F");
254
255	fStree->Branch("TrueImpactX", &fTrack.fTrueImpactX, "fTrueImpactX/F");
256	fStree->Branch("TrueImpactY", &fTrack.fTrueImpactY, "fTrueImpactY/F");
257	fStree->Branch("TrueImpactZ", &fTrack.fTrueImpactZ, "fTrueImpactZ/F");
258	fStree->Branch("TrueMaxPosX", &fTrack.fTrueMaxPosX, "fTrueMaxPosX/F");
259	fStree->Branch("TrueMaxPosY", &fTrack.fTrueMaxPosY, "fTrueMaxPosY/F");
260	fStree->Branch("TrueMaxPosZ", &fTrack.fTrueMaxPosZ, "fTrueMaxPosZ/F");
261	fStree->Branch("TrueHmax", &fTrack.fTrueHmax, "fTrueHmax/F");
262	fStree->Branch("TrueXmax", &fTrack.fTrueXmax, "fTrueXmax/F");
263	fStree->Branch("TrueDmax", &fTrack.fTrueDmax, "fTrueDmax/F");
264
265	fStree->Branch("Latitude", &fTrack.fLatitude, "fLatitude/F");
266	fStree->Branch("Longitude", &fTrack.fLongitude, "fLongitude/F");
267	fStree->Branch("Date", &fTrack.fDate, "fDate/F");
268	fStree->Branch("Hclouds", &fTrack.fHclouds, "fHclouds/F");
269	fStree->Branch("Clouds_thick", &fTrack.fCloudsThick, "fCloudsThick/F");
270	fStree->Branch("Clouds_OD", &fTrack.fCloudsOD, "fCloudsOD/F");
271
272	// shower
273	fStree->Branch("Ne", &fTrack.fNe, "fNe/F");
274	fStree->Branch("Ne2", &fTrack.fNe2, "fNe2/F");
275	fStree->Branch("Nemax", &fTrack.fNemax, "fNemax/F");
276	fStree->Branch("Width_show", &fTrack.fWidthshow, "fWidthshow/F");
277	fStree->Branch("Width2_show", &fTrack.fWidth2show, "fWidth2show/F");
278	fStree->Branch("Xmax_show", &fTrack.fXmaxshow, "fXmaxshow/F");
279	fStree->Branch("Hmax_show", &fTrack.fHmaxshow, "fHmaxshow/F");
280	fStree->Branch("Lmax_show", &fTrack.fLmaxshow, "fLmaxshow/F");
281
282	// photons at creation
283	fStree->Branch("Nph_fluo_L", &fTrack.fNph_f_L, "fNph_f_L/F");
284	fStree->Branch("Nph2_fluo_L", &fTrack.fNph2_f_L, "fNph2_f_L/F");
285	fStree->Branch("Nmax_fluo_L", &fTrack.fNmax_f_L, "fNmax_f_L/F");
286	fStree->Branch("Width_fluo", &fTrack.fWidth_f, "fWidth_f/F");
287	fStree->Branch("Width2_fluo", &fTrack.fWidth2_f, "fWidth2_f/F");
288	fStree->Branch("Hmax_fluo", &fTrack.fHmax_f, "fHmax_f/F");
289	fStree->Branch("Lmax_fluo", &fTrack.fLmax_f, "fLmax_f/F");
290	fStree->Branch("xmax_fluo", &fTrack.fxmax_f, "fxmax_f/F");
291	fStree->Branch("ymax_fluo", &fTrack.fymax_f, "fymax_f/F");
292	fStree->Branch("Xmaxph_fluo", &fTrack.fXmaxph_f, "fXmaxph_f/F");
293	fStree->Branch("Yield_f", &fTrack.fYieldmean, "fYieldmean/F");
294	fStree->Branch("Yieldmax_fluo", &fTrack.fYieldHmax, "fYieldHmax/F");
295
296	fStree->Branch("Nph_cer_L", &fTrack.fNph_c_L, "fNph_c_L/F");
297	fStree->Branch("Nph2_cer_L", &fTrack.fNph2_c_L, "fNph2_c_L/F");
298	fStree->Branch("Hmax_cer", &fTrack.fHmax_c, "fHmax_c/F");
299	fStree->Branch("Yieldmax_ckov", &fTrack.fYieldHmax_c, "fYieldHmax_c/F");
300
301	// old and/or so far useless data
302	//fStree->Branch("Nph_fluo", &fTrack.fNph_f, "fNph_f/F");
303	//fStree->Branch("Nph2_fluo", &fTrack.fNph2_f, "fNph2_f/F");
304	//fStree->Branch("Nmax_fluo", &fTrack.fNmax_f, "fNmax_f/F");
305
306	//fStree->Branch("Nph_cer", &fTrack.fNph_c, "fNph_c/F");
307	//fStree->Branch("Nph2_cer", &fTrack.fNph2_c, "fNph2_c/F");
308	//fStree->Branch("Nmax_cer", &fTrack.fNmax_c, "fNmax_c/F"); // not used
309	//fStree->Branch("Xmaxph_cer", &fTrack.fXmaxph_c, "fXmaxph_c/F"); // not used
310
311	//fStree->Branch("Wlmean_f", &fTrack.fWlmean_f, "fWlmean_f/F"); // not used
312	//fStree->Branch("Wlmean_c", &fTrack.fWlmean_c, "fWlmean_c/F"); // not used
313	//fStree->Branch("Wlmean_tot", &fTrack.fWlmean_tot, "fWlmean_tot/F"); // not used
314
315	// on pupil
316	fStree->Branch("FluoInOmega", &fTrack.fFluoInOmega, "fFluoInOmega/F");
317	fStree->Branch("CkovInOmega", &fTrack.fCkovInOmega, "fCkovInOmega/F");
318	fStree->Branch("Nph_fluo_pupil", &fTrack.fNph_f_p, "fNph_f_p/F");
319	fStree->Branch("Nph_cer_dir_pupil", &fTrack.fNph_cdirect_p, "fNph_cdirect_p/F");
320	fStree->Branch("Nph2_fluo_pupil", &fTrack.fNph2_f_p, "fNph2_f_p/F");
321	fStree->Branch("Nph_cer_air_pupil", &fTrack.fNph_airscat_p, "fNph_airscat_p/F");
322	fStree->Branch("Nph_cer_cloud_pupil",&fTrack.fNph_cloudscat_p,"fNph_cloudscat_p/F");
323	fStree->Branch("Nph_cer_refl_pupil", &fTrack.fNph_cr_p, "fNph_cr_p/F");
324	fStree->Branch("Nmax_fluo2_pupil", &fTrack.fNmax_f2_p, "fNmax_f2_p/F");
325	fStree->Branch("Nmax_fluo_pupil", &fTrack.fNmax_f_p, "fNmax_f_p/F");
326	fStree->Branch("Nmax_fluo_pupil_raw",&fTrack.fNmax_f_p_raw, "fNmax_f_p_raw/F");
327	fStree->Branch("Nmax_tot_pupil", &fTrack.fNmax_tot_p, "fNmax_tot_p/F");
328	fStree->Branch("Nmax_tot_pupil_raw", &fTrack.fNmax_tot_p_raw, "fNmax_tot_p_raw/F");
329	fStree->Branch("Nmax_cr_pupil", &fTrack.fNmax_cr_p, "fNmax_cr_p/F");
330	fStree->Branch("Nmax_crtot_pupil", &fTrack.fNmax_crtot_p, "fNmax_crtot_p/F");
331
332	fStree->Branch("Width_fluo2_pupil", &fTrack.fWidth_f2_p, "fWidth_f2_p/F");
333	fStree->Branch("Width_fluo_pupil", &fTrack.fWidth_f_p, "fWidth_f_p/F");
334	fStree->Branch("Width_tot_pupil", &fTrack.fWidth_tot_p, "fWidth_tot_p/F");
335	fStree->Branch("Yield_fluoGTUmax", &fTrack.fYieldGTUmax, "fYieldGTUmax/F");
336	fStree->Branch("NbBunchInGTUmax", &fTrack.fNbBunches, "fNbBunches/I");
337	fStree->Branch("X_GTUmax", &fTrack.fXGTUmax, "fXGTUmax/F");
338	fStree->Branch("X_GTUmax_tot", &fTrack.fXGTUmax_tot, "fXGTUmax_tot/F");
339	fStree->Branch("Trans_fluo_atMax", &fTrack.fTransMax, "fTransMax/F");
340	fStree->Branch("Trans_fluo_GTUmax", &fTrack.fTransGTUmax, "fTransGTUmax/F");
341	fStree->Branch("Wlmean_f_p", &fTrack.fWlmean_f_p, "fWlmean_f_p/F");
342	fStree->Branch("Wlmean_cb_p", &fTrack.fWlmean_cb_p, "fWlmean_cb_p/F");
343	fStree->Branch("Wlmean_cr_p", &fTrack.fWlmean_cr_p, "fWlmean_cr_p/F");
344	fStree->Branch("Wlmean_tot_p", &fTrack.fWlmean_tot_p, "fWlmean_tot_p/F");
345	}
346	}
347
348	//_____________________________________________________________________________
349	void ETreePainter::FillTrack() {
350	//
351	// fill fTrack
352	//
353
354	// fill truth data
355	if(fTruth) {
356	fTrack.fEnergy = fTruth->GetTrueEnergy(); // en MeV
357	fTrack.fTheta = fTruth->GetTrueTheta()*TMath::RadToDeg();
358	fTrack.fPhi = fTruth->GetTruePhi()*TMath::RadToDeg();
359
360	fTrack.fH1 = Zv(fTruth->GetTrueInitPos())/ km;
361	fTrack.fInitPosX = fTruth->GetTrueInitPos().X()/km;
362	fTrack.fInitPosY = fTruth->GetTrueInitPos().Y()/km;
363	fTrack.fInitPosZ = fTruth->GetTrueInitPos().Z()/km;
364	fTrack.fX1 = fTruth->GetTrueX1()/g*cm2;
365
366	fTrack.fTrueImpactX = fTruth->GetTrueEarthImpact().X()/km;
367	fTrack.fTrueImpactY = fTruth->GetTrueEarthImpact().Y()/km;
368	fTrack.fTrueImpactZ = fTruth->GetTrueEarthImpact().Z()/km;
369	fTrack.fTrueMaxPosX = fTruth->GetTrueShowerMaxPos().X()/km;
370	fTrack.fTrueMaxPosY = fTruth->GetTrueShowerMaxPos().Y()/km;
371	fTrack.fTrueMaxPosZ = fTruth->GetTrueShowerMaxPos().Z()/km;
372	fTrack.fTrueHmax = Zv(fTruth->GetTrueShowerMaxPos())/km;
373	fTrack.fTrueXmax = fTruth->GetTrueShowerXMax()/g*cm2;
374	fTrack.fTrueDmax = (fEUSO - fTruth->GetTrueShowerMaxPos()).Mag() / km;
375
376	fTrack.fLatitude = fTruth->GetTrueLatitude();
377	fTrack.fLongitude = fTruth->GetTrueLongitude();
378	fTrack.fDate = fTruth->GetTrueDate();
379	fTrack.fHclouds = fTruth->GetTrueHclouds() /km;
380	fTrack.fCloudsThick = fTruth->GetTrueCloudsthick() /km;
381	fTrack.fCloudsOD = fTruth->GetTrueCloudsOD();
382
383	}
384	else Info("FillTrack","ETruth object EMPTY");
385
386	if(fShower && fAtmo) {
387
388	TH1F* th(0);
389	TF1* mypol2 = 0;
390	TF1* mygaus = 0;
391	Double_t par[10];
392	Double_t TimeMin = 0.;
393	Double_t TimeMax = 0.;
394	Int_t maxbin = 0;
395	Double_t maxpos = 0;
396	Int_t leftbin = 1;
397	Int_t rightbin = 1;
398	Int_t NbinsAroundMax = 0; // to avoid step effects when searching Xmax, Hmax etc.. (fit around maximum)
399	Bool_t DoFit = false; // if true -> curves are smoothed before determination of plot features
400	fTrack.fGood_fit = 1;
401
402
403	////////////////////////////////////////////////////////////
404	////////////////////// SHOWER DATA /////////////////////////
405	////////////////////////////////////////////////////////////
406
407	// init
408	Int_t numsteps = fShower->GetNumSteps();
409	TVector3 TrackDir = fShower->GetDir(); // track direction in MES
410	EShowerStep* shstep = 0;
411	Int_t indexofmax = 0; //TOFIX : also used for fTransMax
412	Float_t NORM_weight = 0.;
413	Bool_t flag = false;
414	fTrack.fYieldmean = 0.;
415	//Float_t integrate = 0.;
416	//Float_t integrate2 = 0.;
417
418
419
420	// poly2 fit of shower profile (longitudinal - in g/cm2) around its maxumum -> Xmax_shower Nemax Ne Ne2
421	th = BuildLongitudinalHisto("Shower_Longit_prof_gram");
422	leftbin = 1;
423	NbinsAroundMax = Int_t(30. / th->GetBinWidth(th->GetMaximumBin())) + 1; // 30 g/cm2 seems ok
424	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
425	rightbin = th->GetNbinsX();
426	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
427	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
428	th->Fit("mypol2","RQON");
429	mypol2->GetParameters(par);
430	fTrack.fNe = th->Integral("width");
431	maxpos = -0.5*par[1] / par[2];
432	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
433	fTrack.fNemax = par[0] - 0.25par[1]par[1]/par[2];
434	fTrack.fXmaxshow = -0.5*par[1] / par[2];
435	maxbin = th->FindBin(fTrack.fXmaxshow);
436	fTrack.fX_max_last_bin = th->GetNbinsX() - maxbin;
437	fTrack.fNe2 = th->Integral(1,maxbin - 1,"width");
438	fTrack.fNe2 += th->GetBinContent(maxbin) * (fTrack.fXmaxshow - th->GetBinLowEdge(maxbin));
439	}
440	else {
441	fTrack.fGood_fit = 0;
442	maxbin = th->GetMaximumBin();
443	fTrack.fXmaxshow = th->GetBinCenter(maxbin);
444	fTrack.fX_max_last_bin = th->GetNbinsX() - maxbin;
445	if(DoFit) Warning("FillTrack()","Fit PB <fXmaxshow>, par[2] is zero -> max_last_bin = %d",fTrack.fX_max_last_bin);
446	fTrack.fNemax = th->GetMaximum();
447	fTrack.fNe2 = th->Integral(1,maxbin - 1,"width");
448	fTrack.fNe2 += th->GetBinContent(maxbin) * 0.5*th->GetBinWidth(maxbin);
449	}
450	delete th; th = 0;
451	delete mypol2; mypol2 = 0;
452
453
454
455
456
457	// shower Lmax + width
458	th = BuildLongitudinalHisto("Shower_Longit_prof_km");
459	leftbin = 1;
460	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
461	rightbin = th->GetNbinsX();
462	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
463	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
464	th->Fit("mypol2","RQON");
465	mypol2->GetParameters(par);
466	maxpos = -0.5*par[1] / par[2];
467	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
468	fTrack.fLmaxshow = -0.5*par[1] / par[2];
469	maxbin = th->FindBin(fTrack.fLmaxshow);
470	fTrack.fL_max_last_bin = th->GetNbinsX() - maxbin;
471	}
472	else {
473	fTrack.fGood_fit = 0;
474	maxbin = th->GetMaximumBin();
475	fTrack.fLmaxshow = th->GetBinCenter(maxbin);
476	fTrack.fL_max_last_bin = th->GetNbinsX() - maxbin;
477	if(DoFit) Warning("FillTrack()","Fit PB <fLmaxshow>, par[2] is zero -> max_last_bin = %d",fTrack.fL_max_last_bin);
478	}
479
480	// shower width
481	leftbin = 1;
482	rightbin = th->GetNbinsX();
483	Bool_t stopleft(false), stopright(false);
484	maxpos = -0.5*par[1] / par[2];
485	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
486	maxpos = par[0] - 0.25par[1]par[1]/par[2];
487	for(Int_t i=0; i<th->GetNbinsX(); i++) {
488	if((th->GetBinContent(i+1) < 0.5*maxpos) && (!stopleft)) leftbin = i+1;
489	else stopleft = true; // bin found
490	if((th->GetBinContent(th->GetNbinsX() - i) < 0.5*maxpos) && (!stopright)) rightbin = th->GetNbinsX() - i;
491	else stopright = true;
492	}
493	// for linear interpolation within a step
494	/* does not work so far : not fit used instead
495	leftx = th->GetBinLowEdge(1);
496	rightx = th->GetBinLowEdge(th->GetNbinsX());
497	if(rightbin > 1) rightx = th->GetBinLowEdge(rightbin) + (0.5maxpos - th->GetBinContent(rightbin-1))/(th->GetBinContent(rightbin-1) - th->GetBinContent(rightbin))th->GetBinWidth(rightbin);
498	if(leftbin < th->GetNbinsX()) leftx = th->GetBinLowEdge(leftbin+1) + (0.5maxpos - th->GetBinContent(leftbin+1))/(th->GetBinContent(leftbin) - th->GetBinContent(leftbin+1))th->GetBinWidth(leftbin);
499	fTrack.fWidthshow = rightx - leftx;
500	fTrack.fW_max_last_bin = th->GetNbinsX() - th->FindBin(rightx);
501	if(fTrack.fWidthshow < 0) {
502	Warning("FillTrack()","Shower WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
503	fTrack.fWidthshow = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
504	fTrack.fW_max_last_bin = -10; // -10 <=> there is a pb
505	}
506	fTrack.fWidth2show = fTrack.fLmaxshow - leftx;
507	*/
508	fTrack.fWidthshow = th->GetBinLowEdge(rightbin) + th->GetBinWidth(rightbin) - th->GetBinLowEdge(leftbin);
509	fTrack.fW_max_last_bin = th->GetNbinsX() - rightbin; // no linear interpolation
510	if(fTrack.fWidthshow < 0) {
511	Warning("FillTrack()","Shower WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
512	fTrack.fWidthshow = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
513	fTrack.fW_max_last_bin = -10; // -10 <=> there is a pb
514	}
515	fTrack.fWidth2show = th->GetBinLowEdge(maxbin) + 0.5*th->GetBinWidth(maxbin) - th->GetBinLowEdge(leftbin);
516	}
517	else {
518	maxpos = th->GetBinContent(maxbin);
519	for(Int_t i=0; i<th->GetNbinsX(); i++) {
520	if((th->GetBinContent(i+1) < 0.5*maxpos) && (!stopleft)) leftbin = i+1;
521	else stopleft = true; // bin found
522	if((th->GetBinContent(th->GetNbinsX() - i) < 0.5*maxpos) && (!stopright)) rightbin = th->GetNbinsX() - i;
523	else stopright = true;
524	}
525	fTrack.fWidthshow = th->GetBinLowEdge(rightbin) + th->GetBinWidth(rightbin) - th->GetBinLowEdge(leftbin);
526	fTrack.fW_max_last_bin = th->GetNbinsX() - rightbin; // no linear interpolation
527	if(fTrack.fWidthshow < 0) {
528	Warning("FillTrack()","Shower WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
529	fTrack.fWidthshow = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
530	fTrack.fW_max_last_bin = -10; // -10 <=> there is a pb
531	}
532	fTrack.fWidth2show = th->GetBinLowEdge(maxbin) + 0.5*th->GetBinWidth(maxbin) - th->GetBinLowEdge(leftbin);
533	}
534	delete th; th = 0;
535	delete mypol2; mypol2 = 0;
536
537
538
539
540	// gaussian fit of shower profile (vs altitude - in km) around its maxumum -> Hmax_shower
541	th = BuildAltHisto("alt_Shower");
542	leftbin = 1;
543	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
544	rightbin = th->GetNbinsX();
545	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
546	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
547	th->Fit("mypol2","RQON");
548	mypol2->GetParameters(par);
549	maxpos = -0.5*par[1] / par[2];
550	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
551	fTrack.fHmaxshow = -0.5*par[1] / par[2];
552	maxbin = th->FindBin(fTrack.fHmaxshow);
553	fTrack.fH_max_last_bin = maxbin - 1;
554	}
555	else {
556	fTrack.fGood_fit = 0;
557	maxbin = th->GetMaximumBin();
558	fTrack.fHmaxshow = th->GetBinCenter(maxbin);
559	fTrack.fH_max_last_bin = maxbin - 1;
560	if(DoFit) Warning("FillTrack()","Fit PB <fHmaxshow>, par[2] is zero -> max_last_bin = %d",fTrack.fH_max_last_bin);
561	}
562	delete th; th = 0;
563	delete mypol2; mypol2 = 0;
564
565
566
567	// NOW : kept only for yield calculations
568	// (OLD : integrate shower profile - find the showerstep of Ne maximum - calculate mean fluo yield (fYieldmean))
569	for(Int_t i=0; i<numsteps; i++) {
570	shstep = fShower->GetStep(i);
571	if(!flag) {
572	if((shstep->GetXf()/g*cm2) > fTrack.fXmaxshow) {
573	indexofmax = i;
574	flag = true;
575	}
576	}
577	fTrack.fYieldmean += fAtmo->GetBunch(2i)->GetYield()m * shstep->GetNumElectrons();
578	NORM_weight += shstep->GetNumElectrons();
579	}
580	if(NORM_weight) fTrack.fYieldmean /= NORM_weight;
581	fTrack.fYieldHmax = fAtmo->GetBunch(2indexofmax)->GetYield()m; // get fluo yield associated to shower max
582	fTrack.fYieldHmax_c = fAtmo->GetBunch(2indexofmax+1)->GetYield()m; // get fluo yield associated to shower max
583
584	const Int_t indexofmax_cst = indexofmax;
585
586
587
588
589
590	////////////////////////////////////////////////////////////
591	////////////////////// ATMOSPHERE DATA /////////////////////
592	////////////////////////////////////////////////////////////
593
594	///////////////////// AT PHOTONS CREATION ///////////////////
595
596	// from data at creation VS atlitude
597	// gaussian fit profiles around their maxumum -> Hmax
598	th = BuildAltHisto("alt_Bunch_Nph_F");
599	leftbin = 1;
600	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
601	rightbin = th->GetNbinsX();
602	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
603	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
604	th->Fit("mypol2","RQON");
605	mypol2->GetParameters(par);
606	maxpos = -0.5*par[1] / par[2];
607	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
608	fTrack.fHmax_f = -0.5*par[1] / par[2];
609	maxbin = th->FindBin(fTrack.fHmax_f);
610	fTrack.fH_fluomax_last_bin = maxbin - 1;
611	}
612	else {
613	fTrack.fGood_fit = 0;
614	maxbin = th->GetMaximumBin();
615	fTrack.fHmax_f = th->GetBinCenter(maxbin);
616	fTrack.fH_fluomax_last_bin = maxbin - 1;
617	if(DoFit) Warning("FillTrack()","Fit PB <fHmax_fluo>, par[2] is zero -> max_last_bin = %d",fTrack.fH_fluomax_last_bin);
618	}
619	delete th; th = 0;
620	delete mypol2; mypol2 = 0;
621
622
623	th = BuildAltHisto("alt_Bunch_Nph_C");
624	leftbin = 1;
625	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
626	rightbin = th->GetNbinsX();
627	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
628	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
629	th->Fit("mypol2","RQON");
630	mypol2->GetParameters(par);
631	fTrack.fH_ckovmax_last_bin = th->FindBin(par[1]) - 1;
632	maxpos = -0.5*par[1] / par[2];
633	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
634	fTrack.fHmax_c = -0.5*par[1] / par[2];
635	maxbin = th->FindBin(fTrack.fHmax_c);
636	fTrack.fH_ckovmax_last_bin = maxbin - 1;
637	}
638	else {
639	fTrack.fGood_fit = 0;
640	maxbin = th->GetMaximumBin();
641	fTrack.fHmax_c = th->GetBinCenter(maxbin);
642	fTrack.fH_ckovmax_last_bin = maxbin - 1;
643	if(DoFit) Warning("FillTrack()","Fit PB <fHmax_ckov>, par[2] is zero -> max_last_bin = %d",fTrack.fH_ckovmax_last_bin);
644	}
645	delete th; th = 0;
646	delete mypol2; mypol2 = 0;
647
648
649
650	// from longitudinal data in grammage
651	// - nb - HALF nb - nb at max - Xmax - of photons at creation
652	th = BuildLongitudinalHisto("Fluo_Longit_prof_gram");
653	leftbin = 1;
654	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
655	rightbin = th->GetNbinsX();
656	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
657	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
658	th->Fit("mypol2","RQON");
659	mypol2->GetParameters(par);
660	maxpos = -0.5*par[1] / par[2];
661	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
662	fTrack.fXmaxph_f = -0.5*par[1] / par[2];
663	maxbin = th->FindBin(fTrack.fXmaxph_f);
664	fTrack.fX_fluomax_last_bin = th->GetNbinsX() - maxbin;
665	// NOT USED ANYMORE : to be adapted (with par[2] and Dofit) if used in the future
666	//fTrack.fNmax_f = par[0];
667	//fTrack.fNmax_f = th->GetMaximum(); //nofit
668	//fTrack.fNph_f = th->Integral("width");
669	//fTrack.fNph2_f = th->Integral(1,maxbin - 1,"width");
670	//fTrack.fNph2_f += th->GetBinContent(maxbin) * (fTrack.fXmaxph_f - th->GetBinLowEdge(maxbin));
671	//fTrack.fNph2_f += th->GetBinContent(th->GetMaximumBin()) * 0.5*th->GetBinWidth(th->GetMaximumBin()); //nofit
672	}
673	else {
674	fTrack.fGood_fit = 0;
675	maxbin = th->GetMaximumBin();
676	fTrack.fXmaxph_f = th->GetBinCenter(maxbin);
677	fTrack.fX_fluomax_last_bin = th->GetNbinsX() - maxbin;
678	if(DoFit) Warning("FillTrack()","Fit PB <fXmax_fluo>, par[2] is zero -> max_last_bin = %d",fTrack.fX_fluomax_last_bin);
679	}
680	delete th; th = 0;
681	delete mypol2; mypol2 = 0;
682
683
684
685	// idem for ckov
686	/* not used so far : to be adapted for fit (with par[2] and Dofit) if used in the future
687	th = BuildLongitudinalHisto("Ckov_Longit_prof_gram");
688	leftbin = 1;
689	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
690	rightbin = th->GetNbinsX();
691	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
692	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
693	th->Fit("mypol2","RQON");
694	mypol2->GetParameters(par);
695	//fTrack.fXmaxph_c = par[1]; //DELETE
696	//fTrack.fXmaxph_c = th->GetBinCenter(th->GetMaximumBin());
697	//maxbin = th->FindBin(fTrack.fXmaxph_c); //DELETE
698	maxbin = th->GetMaximumBin();
699	//fTrack.fNmax_c = par[0]; //DELETE
700	//fTrack.fNmax_c = th->GetMaximum();
701	//fTrack.fNph_c = th->Integral("width");
702	//fTrack.fNph2_c = th->Integral(1,maxbin - 1,"width");
703	//fTrack.fNph2_c += th->GetBinContent(maxbin) * (fTrack.fXmaxph_c - th->GetBinLowEdge(maxbin));
704	//fTrack.fNph2_c += th->GetBinContent(maxbin) * 0.5*th->GetBinWidth(th->GetMaximumBin()); //DELETE
705	fTrack.fX_ckovmax_last_bin = th->GetNbinsX() - th->GetMaximumBin();//TEMP
706	delete th; th = 0;
707	delete mypol2; mypol2 = 0;
708	*/
709
710
711
712	// Nmax_fluo calculated per km - Lmax, xmax, ymax
713	th = BuildLongitudinalHisto("Fluo_Longit_prof_km");
714	leftbin = 1;
715	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
716	rightbin = th->GetNbinsX();
717	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
718	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
719	th->Fit("mypol2","RQON");
720	mypol2->GetParameters(par);
721	fTrack.fNph_f_L = th->Integral("width");
722	maxpos = -0.5*par[1] / par[2];
723	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
724	fTrack.fNmax_f_L = par[0] - 0.25 * par[1] * par[1] / par[2];
725	fTrack.fLmax_f = -0.5 * par[1] / par[2];
726	maxbin = th->FindBin(fTrack.fLmax_f);
727	fTrack.fL_fluomax_last_bin = th->GetNbinsX() - maxbin;
728	fTrack.fNph2_f_L = th->Integral(1,maxbin - 1,"width");
729	fTrack.fNph2_f_L += th->GetBinContent(maxbin) * (fTrack.fLmax_f - th->GetBinLowEdge(maxbin));
730	}
731	else {
732	fTrack.fGood_fit = 0;
733	maxbin = th->GetMaximumBin();
734	fTrack.fLmax_f = th->GetBinCenter(maxbin);
735	fTrack.fL_fluomax_last_bin = th->GetNbinsX() - maxbin;
736	if(DoFit) Warning("FillTrack()","Fit PB <fLmax_fluo>, par[2] is zero -> max_last_bin = %d",fTrack.fL_fluomax_last_bin);
737	fTrack.fNmax_f_L = th->GetMaximum(); //nofit
738	fTrack.fNph2_f_L = th->Integral(1,maxbin - 1,"width");
739	fTrack.fNph2_f_L += th->GetBinContent(maxbin) * 0.5*th->GetBinWidth(maxbin);
740	}
741	fTrack.fxmax_f = fTrack.fInitPosX + fTrack.fLmax_f*TrackDir.X();
742	fTrack.fymax_f = fTrack.fInitPosY + fTrack.fLmax_f*TrackDir.Y();
743
744
745	// width - HALF width - of fluo profile (km)
746	//maxbin = th->GetMaximumBin(); //nofit
747	//maxpos = th->GetBinContent(maxbin); //nofit
748	leftbin = 1;
749	rightbin = th->GetNbinsX();
750	stopleft = false;
751	stopright = false;
752	maxpos = -0.5*par[1] / par[2];
753	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
754	maxpos = par[0] - 0.25par[1]par[1]/par[2];
755	for(Int_t i=0; i<th->GetNbinsX(); i++) {
756	if((th->GetBinContent(i+1) < 0.5*maxpos) && (!stopleft)) leftbin = i+1;
757	else stopleft = true; // bin found
758	if((th->GetBinContent(th->GetNbinsX() - i) < 0.5*maxpos) && (!stopright)) rightbin = th->GetNbinsX() - i;
759	else stopright = true;
760	}
761	// for linear interpolation within a step
762	/* does not work so far : no fits used instead
763	leftx = th->GetBinLowEdge(1);
764	rightx = th->GetBinLowEdge(th->GetNbinsX());
765	if(rightbin > 1) rightx = th->GetBinLowEdge(rightbin) + (0.5maxpos - th->GetBinContent(rightbin-1))/(th->GetBinContent(rightbin-1) - th->GetBinContent(rightbin))th->GetBinWidth(rightbin);
766	if(leftbin < th->GetNbinsX()) leftx = th->GetBinLowEdge(leftbin+1) + (0.5maxpos - th->GetBinContent(leftbin+1))/(th->GetBinContent(leftbin) - th->GetBinContent(leftbin+1))th->GetBinWidth(leftbin);
767	fTrack.fWidth_f = rightx - leftx;
768	fTrack.fW_fluomax_last_bin = th->GetNbinsX() - th->FindBin(rightx);
769	if(fTrack.fWidth_f < 0) {
770	Warning("FillTrack()","FLuo WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
771	fTrack.fWidth_f = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
772	fTrack.fW_fluomax_last_bin = -10; // -10 <=> there is a pb
773	}
774	fTrack.fWidth2_f = fTrack.fLmax_f - leftx;
775	*/
776	fTrack.fWidth_f = th->GetBinLowEdge(rightbin) + th->GetBinWidth(rightbin) - th->GetBinLowEdge(leftbin);
777	fTrack.fW_fluomax_last_bin = th->GetNbinsX() - rightbin; // no linear interpolation
778	if(fTrack.fWidth_f < 0) {
779	Warning("FillTrack()","Fluo WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
780	fTrack.fWidth_f = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
781	fTrack.fW_fluomax_last_bin = -10; // -10 <=> there is a pb
782	}
783	fTrack.fWidth2_f = th->GetBinLowEdge(maxbin) + 0.5*th->GetBinWidth(maxbin) - th->GetBinLowEdge(leftbin);
784	}
785	else {
786	maxpos = th->GetBinContent(maxbin);
787	for(Int_t i=0; i<th->GetNbinsX(); i++) {
788	if((th->GetBinContent(i+1) < 0.5*maxpos) && (!stopleft)) leftbin = i+1;
789	else stopleft = true; // bin found
790	if((th->GetBinContent(th->GetNbinsX() - i) < 0.5*maxpos) && (!stopright)) rightbin = th->GetNbinsX() - i;
791	else stopright = true;
792	}
793	fTrack.fWidth_f = th->GetBinLowEdge(rightbin) + th->GetBinWidth(rightbin) - th->GetBinLowEdge(leftbin);
794	fTrack.fW_fluomax_last_bin = th->GetNbinsX() - rightbin; // no linear interpolation
795	if(fTrack.fWidth_f < 0) {
796	Warning("FillTrack()","Fluo WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
797	fTrack.fWidth_f = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
798	fTrack.fW_fluomax_last_bin = -10; // -10 <=> there is a pb
799	}
800	fTrack.fWidth2_f = th->GetBinLowEdge(maxbin) + 0.5*th->GetBinWidth(maxbin) - th->GetBinLowEdge(leftbin);
801	}
802	delete th; th = 0;
803	delete mypol2; mypol2 = 0;
804
805
806
807
808	// Ckov Nph_L and Nph2_L calculated per km
809	th = BuildLongitudinalHisto("Ckov_Longit_prof_km");
810	leftbin = 1;
811	if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
812	rightbin = th->GetNbinsX();
813	if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
814	mypol2 = new TF1("mypol2","pol2",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
815	th->Fit("mypol2","RQON");
816	mypol2->GetParameters(par);
817	fTrack.fNph_c_L = th->Integral("width");
818	maxpos = -0.5*par[1] / par[2];
819	if(par[2] && DoFit && (maxpos <th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin)) && (maxpos > th->GetBinLowEdge(leftbin))) {
820	Double_t localmax = -0.5 * par[1] / par[2];
821	maxbin = th->FindBin(localmax);
822	fTrack.fNph2_c_L = th->Integral(1,maxbin - 1,"width");
823	fTrack.fNph2_c_L += th->GetBinContent(maxbin) * (localmax - th->GetBinLowEdge(maxbin));
824	}
825	else {
826	fTrack.fGood_fit = 0;
827	maxbin = th->GetMaximumBin();
828	if(DoFit) Warning("FillTrack()","Fit PB <fLmax_ckov>, par[2] is zero -> max_last_bin = %d", th->GetNbinsX() - maxbin);
829	fTrack.fNph2_c_L = th->Integral(1,maxbin - 1,"width");
830	fTrack.fNph2_c_L += th->GetBinContent(maxbin) * 0.5*th->GetBinWidth(maxbin);
831	}
832	delete th; th = 0;
833	delete mypol2; mypol2 = 0;
834
835
836
837
838	////////////////////// ON PUPIL /////////////////////////////
839
840	// from data on pupil VS time
841	// fluo
842	Double_t tmax_pupil = 0.;
843	th = BuildTofHisto("Single_Nph_F_P_t");
844	fTrack.fNph_f_p = th->Integral();
845	fTrack.fNmax_f_p_raw = th->GetMaximum();
846	TimeMin = th->GetBinLowEdge(1);
847	TimeMax = th->GetBinLowEdge(th->GetXaxis()->GetLast())+th->GetBinWidth(th->GetXaxis()->GetLast());
848	mygaus = new TF1("mygaus","gaus",TimeMin,TimeMax);
849	th->Fit("mygaus","RQON");
850	mygaus->GetParameters(par);
851	fTrack.fNph2_f_p = th->Integral(1,th->FindBin(par[1])); // first HALF of profile
852	fTrack.fNmax_f_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
853	tmax_pupil = par[1]*microsecond; // records absolute Tmax in ms for below search of fXGTUmax - fYieldGTUmax etc..
854	fTrack.fWidth_f_p = par[2]/(fGTU/microsecond);
855	delete th; th = 0;
856	delete mygaus; mygaus = 0;
857
858	// direct ckov
859	th = BuildTofHisto("Single_Nph_Cdir_P_t");
860	fTrack.fNph_cdirect_p = th->Integral();
861	delete th; th = 0;
862
863	// ground reflected ckov
864	th = BuildTofHisto("Single_Nph_CR_P_t");
865	Int_t localbin = th->GetMaximumBin();
866	Double_t time_cr_max = th->GetBinCenter(localbin); // used here below for fTrack.fNmax_crtot_p
867	fTrack.fNph_cr_p = th->Integral();
868	fTrack.fNmax_cr_p = th->GetBinContent(localbin)/(fGTU/microsecond);
869	delete th; th = 0;
870
871	// air scattered ckov
872	th = BuildTofHisto("Single_Nph_CB_P_t");
873	fTrack.fNph_airscat_p = th->Integral();
874	fTrack.fNmax_crtot_p = fTrack.fNmax_cr_p + th->GetBinContent(th->FindBin(time_cr_max))/(fGTU/microsecond);
875	delete th; th = 0;
876
877	// clouds scattered ckov
878	th = BuildTofHisto("Single_Nph_CB_cloud_P_t");
879	fTrack.fNph_cloudscat_p = th->Integral();
880	delete th; th = 0;
881
882	// direct fluo + scattered ckov ->
883	// 1. fit total
884	// 2. fit only the first profile half (to substract a part of ckov component)
885	Double_t tmax_tot_pupil = 0.;
886	th = BuildTofHisto("Single_Nph_F-CB_P_t");
887	fTrack.fNmax_tot_p_raw = th->GetMaximum();
888	TimeMin = th->GetBinLowEdge(1);
889	TimeMax = th->GetBinLowEdge(th->GetXaxis()->GetLast())+th->GetBinWidth(th->GetXaxis()->GetLast());
890	mygaus = new TF1("mygaus","gaus",TimeMin,TimeMax);
891	th->Fit("mygaus","RQON");
892	mygaus->GetParameters(par);
893	fTrack.fNmax_tot_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
894	fTrack.fWidth_tot_p = par[2]/(fGTU/microsecond);
895	tmax_tot_pupil = par[1]*microsecond; // records absolute Tmax in ms for below search of fXGTUmax_tot
896	delete mygaus; mygaus = 0;
897	mygaus = new TF1("mygaus",ETreePainter_my_gaussian_fit,TimeMin,TimeMax,4);
898	mygaus->SetParameters(par[0],par[1],par[2],par[1]);
899	mygaus->SetParLimits(0,0,2*par[0]);
900	mygaus->SetParLimits(1,TimeMin,TimeMax);
901	mygaus->SetParLimits(2,0.,(TimeMax - TimeMin));
902	mygaus->SetParLimits(3,par[1],par[1]); // max_tot bound -> for fitting the FIRST half of total profile
903	th->Fit("mygaus","RQON");
904	mygaus->GetParameters(par);
905	fTrack.fNmax_f2_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
906	fTrack.fWidth_f2_p = par[2]/(fGTU/microsecond);
907	delete th; th = 0;
908	delete mygaus; mygaus = 0;
909
910
911
912
913	///////////////// GTUmax RELATED + singles related (mean trans)/////////////////////////
914
915	// Calculate fYieldGTUmax - fXGTUmax - fXGTUmax_tot
916	// fFluoInOmega and fCkovInOmega : get nb of photons BEFORE last TRANSMISSION
917	vector<Int_t> list_bunchID; // ID list of bunches in GTUmax
918	vector<Int_t> list_bunchID_new; // ID list of bunches in GTUmax, one ID appear only once
919	list_bunchID.clear();
920	list_bunchID_new.clear();
921	vector<Int_t> list_bunchID_tot; // same object as here above, but for GTUmax_tot
922	vector<Int_t> list_bunchID_tot_new;
923	list_bunchID_tot.clear();
924	list_bunchID_tot_new.clear();
925	Int_t ns = fAtmo->GetNumSingles();
926	fTrack.fFluoInOmega = 0.;
927	fTrack.fCkovInOmega = 0.;
928	// only FLUO (scattered ckov used only to get shift in time on pupil (tmax_tot_pupil)
929	for (Int_t i=0; i<ns; i++) {
930	ESinglePhoton* s = fAtmo->GetSingle(i);
931	if(s->GetType() == 0) {
932	fTrack.fFluoInOmega++;
933	if ( !s->IsAbsorbed() ) {
934	// get Id of bunches which have generated pupil photons with arrival time within : tmax_pupil +/- (GTU/2)
935	if ( fabs((s->GetTof() + s->GetDate()) - tmax_pupil) <= 0.5*fGTU) {
936	list_bunchID.push_back(Int_t(s->GetBunchId()));
937	}
938	// idem but within : tmax_pupil_tot +/- (GTU/2)
939	if ( fabs((s->GetTof() + s->GetDate()) - tmax_tot_pupil) <= 0.5*fGTU) {
940	list_bunchID_tot.push_back(Int_t(s->GetBunchId()));
941	}
942	}
943	}
944	else if(s->GetType() == 1) fTrack.fCkovInOmega++;
945	}
946
947
948	// to get each ID only once in the ID-lists
949	Bool_t flagg = false; // true if value already exists in new list
950	for(size_t i=0; i<list_bunchID.size(); i++) {
951	flagg = false;
952	for(size_t j=0; j<list_bunchID_new.size(); j++) {
953	if(list_bunchID_new[j] == list_bunchID[i]) flagg = true;
954	}
955	if(!flagg) list_bunchID_new.push_back(list_bunchID[i]);
956	}
957	// idem for ID-lists_tot
958	for(size_t i=0; i<list_bunchID_tot.size(); i++) {
959	flagg = false;
960	for(size_t j=0; j<list_bunchID_tot_new.size(); j++) {
961	if(list_bunchID_tot_new[j] == list_bunchID_tot[i]) flagg = true;
962	}
963	if(!flagg) list_bunchID_tot_new.push_back(list_bunchID_tot[i]);
964	}
965	fTrack.fNbBunches = Int_t(list_bunchID_new.size());
966
967
968	// to get mean transmission (nb of bunches considered == fNbBunches, thus given by GTUmax <-> X translation)
969	// - at Xmax
970	// - at GTUmax
971	Float_t ph_atmax = 0.;
972	Float_t ph_trans_atmax = 0.;
973	Float_t ph_atGTUmax = 0.;
974	Float_t ph_trans_atGTUmax = 0.;
975	Int_t idty = 0;
976
977	// GTUmax photons
978	EBunchPhotons* b = 0;
979	ESinglePhoton* s = 0;
980	// Bunch ID points here to FLUO bunches (look at above lines)
981	for(size_t j=0; j<list_bunchID_new.size(); j++) {
982	b = fAtmo->GetBunch(list_bunchID_new[j]-1); // j-1 because ID j <-> array[j-1]
983	if(b->GetType() != 0) Warning("FillTrack","1. SHOULD be FLUO bunch !!!");
984	ph_atGTUmax += b->GetNumDirect("fluo");
985	TRefArrayIter iter(b->GetSinglesArray());
986	while((s = (ESinglePhoton*)iter.Next()))
987	if(!s->IsAbsorbed()) ph_trans_atGTUmax++; //TOFIX : to be changed if fluo backscattering
988	}
989
990	// (shower-)Xmax photons (nb of bunches considered ('fNbBunches') comes from GTUmax studies, to get ~same statistics in both case)
991	// if idty is around bunch(Xmax) ID <-> if idty is between maxID-2(nb/2) and maxID+2(nb/2 +1 if nb is odd)
992	// "2*" because one fluo bunch step further <-> two bunches ID step further
993	// "nb" is the nb of bunches within GTUmax"
994	// "2*indexofmax_cst" is the ID of showermax bunch
995	idty = 2*indexofmax_cst +1 - fTrack.fNbBunches - fTrack.fNbBunches%2;
996	if(idty < 0) {
997	Warning("FillTrack()","TRONCATED shower ? -> indexofmax = %d, nbBunches = %d",indexofmax_cst,fTrack.fNbBunches);
998	while(idty <= 0) idty ++;
999	Warning("FillTrack()","so idty set to = %d",idty);
1000	}
1001	for(Int_t w=0; w<fTrack.fNbBunches; w++) {
1002	idty += 2*w;
1003	if((idty-1) > fAtmo->GetNumBunches()) {
1004	Warning("FillTrack()","idty is too high -> idty = %d",idty);
1005	fTrack.fNbBunches = w;
1006	Warning("FillTrack()","so fTrack.fNbBunches set to = %d",fTrack.fNbBunches);
1007	break;
1008	}
1009	b = fAtmo->GetBunch(idty-1); // idty-1 because ID idty <-> array[idty-1]
1010	if(!b) { // occurs when end of list is reached
1011	Warning("FillTrack","fAtmo->GetBunch(%d) returned NULL pointer",idty-1);
1012	if(ph_atmax*ph_trans_atmax == 0) { ph_atmax=-1; ph_trans_atmax=-1; }
1013	continue;
1014	}
1015	if(b->GetType() != 0) Warning("FillTrack","2. SHOULD be FLUO bunch !!!");
1016	ph_atmax += b->GetNumDirect("fluo");
1017	TRefArrayIter iter(b->GetSinglesArray());
1018	while((s = (ESinglePhoton*)iter.Next()))
1019	if(!s->IsAbsorbed()) ph_trans_atmax++; //TOFIX : to be changed if fluo backscattering
1020	}
1021
1022	fTrack.fTransMax = ph_trans_atmax/ph_atmax;
1023	fTrack.fTransGTUmax = ph_trans_atGTUmax/ph_atGTUmax;
1024
1025	fTrack.fYieldGTUmax = 0.;
1026	fTrack.fXGTUmax = 0.;
1027	for(size_t j=0; j<list_bunchID_new.size(); j++) {
1028	// mean yield
1029	fTrack.fYieldGTUmax += fAtmo->GetBunch(list_bunchID_new[j] - 1)->GetYield()*m;
1030	// get grammage position from corresponding shower step (bunchID is necessarily of FLUO type)
1031	// XGTUmax is taken at the middle between first and last bunch involved
1032	fTrack.fXGTUmax += 0.5(fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXi())cm2/g;
1033	}
1034	if(list_bunchID_new.size()) {
1035	fTrack.fYieldGTUmax /= Double_t(list_bunchID_new.size());
1036	fTrack.fXGTUmax /= Double_t(list_bunchID_new.size());
1037	}
1038
1039	// idem for fXGTUmax_tot
1040	fTrack.fXGTUmax_tot = 0.;
1041	for(size_t j=0; j<list_bunchID_tot_new.size(); j++) {
1042	// get grammage position from corresponding shower step (bunchID is necessarily of FLUO type)
1043	// XGTUmax is taken at the middle between first and last bunch involved
1044	fTrack.fXGTUmax_tot += 0.5(fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXi())cm2/g;
1045	}
1046	if(list_bunchID_tot_new.size()) fTrack.fXGTUmax_tot /= Double_t(list_bunchID_tot_new.size());
1047
1048
1049
1050	/////////////////// WAVELENGTH //////////////////////////
1051
1052	// from wavelength spectra BEFORE and AFTER transmission
1053	// before
1054	/* mpt really used so far
1055	th = BuildWlHisto("Bunch_Wl_F");
1056	fTrack.fWlmean_f = th->GetMean();
1057	delete th; th = 0;
1058	//th = BuildWlHisto("Bunch_Wl_C");
1059	fTrack.fWlmean_c = th->GetMean();
1060	delete th; th = 0;
1061	th = BuildWlHisto("Bunch_Wl_tot");
1062	fTrack.fWlmean_tot = th->GetMean();
1063	delete th; th = 0;
1064	*/
1065	// after
1066	th = BuildWlHisto("Single_Wl_F_P");
1067	fTrack.fWlmean_f_p = th->GetMean();
1068	delete th; th = 0;
1069	th = BuildWlHisto("Single_Wl_CB_P");
1070	fTrack.fWlmean_cb_p = th->GetMean();
1071	delete th; th = 0;
1072	th = BuildWlHisto("Single_Wl_CR_P");
1073	fTrack.fWlmean_cr_p = th->GetMean();
1074	delete th; th = 0;
1075	th = BuildWlHisto("Single_Wl_tot_P");
1076	fTrack.fWlmean_tot_p = th->GetMean();
1077	delete th; th = 0;
1078	delete mygaus; mygaus = 0;
1079	}
1080
1081	else Info("FillTrack","EAtmosphere OR EShower object EMPTY");
1082
1083	}
1084
1085	//_____________________________________________________________________________
1086	void ETreePainter::FillTree() {
1087	//
1088	// fill stree
1089	//
1090
1091	BuildTree(); // should be useless
1092	fStree->Fill();
1093	}
1094
1095	//_____________________________________________________________________________
1096	void ETreePainter::PrepareHistos(Option_t *opt , Float_t min, Float_t max, UInt_t Nbins) {
1097	//
1098	// define histos features : range, binwidth..
1099	//
1100
1101	if(!fStree) {
1102	cout <<"NO STREE FOR PLOTS !!\n";
1103	return;
1104	}
1105	TString name("h_");
1106	TString option(opt);
1107	name += option.Data();
1108
1109	if(Nbins > 0 && fabs(min*max) > 0) new TH1F(name.Data(),name.Data(),Nbins,min,max);
1110	//else cout <<"\nHistogram '"<<option.Data()<<"' has not been configured\n";
1111	}
1112
1113	//_____________________________________________________________________________
1114	void ETreePainter::DrawHistos( Option_t opt , const TCut& drawcutt, Option_t drawopt, Float_t min, Float_t max, UInt_t Nbins) {
1115	//
1116	//
1117	//
1118
1119	PrepareHistos(opt,min,max,Nbins);
1120	TString option(opt);
1121	TString drawopt_str(drawopt);
1122	TCut drawcut(drawcutt);
1123	drawcut += drawopt_str.Data();
1124
1125	// ************************** RAW HISTOS ****************************
1126
1127	// truth histo
1128	if(option == "Energy") fStree->Draw("TMath::Log10(Energy)>>h_Energy",drawcut,"");
1129	else if(option == "Theta") fStree->Draw("Theta>>h_Theta",drawcut,"");
1130	else if(option == "Phi") fStree->Draw("Phi>>h_Phi",drawcut,"");
1131	else if(option == "H1") fStree->Draw("H1>>h_H1",drawcut,"");
1132	else if(option == "Hmax") fStree->Draw("TrueHmax>>h_Hmax",drawcut,"");
1133	else if(option == "X1") fStree->Draw("X1>>h_X1",drawcut,"");
1134	else if(option == "Xmax") fStree->Draw("TrueXmax>>h_Xmax",drawcut,"");
1135	else if(option == "Xmax_X1") fStree->Draw("TrueXmax - X1>>h_Xmax_X1",drawcut,"");
1136	// shower histo
1137	else if(option == "Nemax") fStree->Draw("Nemax>>h_Nemax",drawcut,"");
1138
1139	// atmosphere histos
1140	else if(option == "Nph_Fluo_L") fStree->Draw("Nph_fluo_L>>h_Nph_Fluo_L",drawcut,"");
1141	else if(option == "Nph_Cer_L") fStree->Draw("Nph_cer_L>>h_Nph_Cer_L",drawcut,"");
1142	else if(option == "Nmax_Fluo") fStree->Draw("Nmax_fluo_L>>h_Nmax_Fluo",drawcut,"");
1143	else if(option == "Nmax_Cer") fStree->Draw("Nmax_cer>>h_Nmax_Cer",drawcut,"");
1144	else if(option == "Hmax_Fluo") fStree->Draw("Hmax_fluo>>h_Hmax_Fluo",drawcut,"");
1145	else if(option == "Hmax_Cer") fStree->Draw("Hmax_cer>>h_Hmax_Cer",drawcut,"");
1146	else if(option == "Xmaxph_Fluo") fStree->Draw("Xmaxph_fluo>>h_Xmaxph_Fluo",drawcut,"");
1147	else if(option == "Xmaxph_Cer") fStree->Draw("Xmaxph_cer>>h_Xmaxph_Cer",drawcut,"");
1148	else if(option == "Wlmean_f") fStree->Draw("Wlmean_f>>h_Wlmean_f",drawcut,"");
1149	else if(option == "Wlmean_c") fStree->Draw("Wlmean_c>>h_Wlmean_c",drawcut,"");
1150	else if(option == "Wlmean_tot") fStree->Draw("Wlmean_tot>>h_Wlmean_tot",drawcut,"");
1151
1152	// on pupil histos
1153	else if(option == "Nph_Fluo_P") fStree->Draw("Nph_fluo_pupil>>h_Nph_Fluo_P",drawcut,"");
1154	else if(option == "Nmax_Fluo_P") fStree->Draw("Nmax_fluo_pupil>>h_Nmax_Fluo_P",drawcut,"");
1155	else if(option == "Sig_Fluo_P") fStree->Draw("Width_fluo_pupil>>h_Sig_Fluo_P",drawcut,"");
1156	//else if(option == "Nph_CB_P") fStree->Draw("Nph_cer_back_pupil>>h_Nph_CB_P",drawcut,"");
1157	else if(option == "Nph_CB_P") fStree->Draw("Nph_cer_air_pupil>>h_Nph_CB_P",drawcut,"");
1158	else if(option == "Nph_CR_P") fStree->Draw("Nph_cer_refl_pupil>>h_Nph_CR_P",drawcut,"");
1159	//else if(option == "Nph_tot_P") fStree->Draw("Nph_fluo_pupil+Nph_cer_back_pupil+Nph_cer_refl_pupil>>h_Nph_tot_P",drawcut,"");
1160	else if(option == "Nph_tot_P") fStree->Draw("Nph_fluo_pupil+Nph_cer_air_pupil+Nph_cer_refl_pupil>>h_Nph_tot_P",drawcut,"");
1161	else if(option == "Wlmean_f_p") fStree->Draw("Wlmean_f_p>>h_Wlmean_f_p",drawcut,"");
1162	else if(option == "Wlmean_cb_p") fStree->Draw("Wlmean_cb_p>>h_Wlmean_cb_p",drawcut,"");
1163	else if(option == "Wlmean_cr_p") fStree->Draw("Wlmean_cr_p>>h_Wlmean_cr_p",drawcut,"");
1164	else if(option == "Wlmean_tot_p") fStree->Draw("Wlmean_tot_p>>h_Wlmean_tot_p",drawcut,"");
1165
1166
1167
1168
1169	// ************************** CORRELATIONS ****************************
1170
1171	// show EAS crash on ground -> it is better to use half integrals
1172	else if(option == "Nph2_Over_Nph_Hmax_fluo") fStree->Draw("Nph2_fluo_L/Nph_fluo_L:Hmax_show>>h_Nph2_Over_Nph_Hmax_fluo",drawcut,"");
1173	else if(option == "Nph2_Over_Nph_Theta_fluo") fStree->Draw("Nph2_fluo_L/Nph_fluo_L:Theta>>h_Nph2_Over_Nph_Theta_fluo",drawcut,"");
1174	else if(option == "Nph2_Over_Nph_Hmax_cer") fStree->Draw("Nph2_cer_L/Nph_cer_L:Hmax_show>>h_Nph2_Over_Nph_Hmax_cer",drawcut,"");
1175	else if(option == "Nph2_Over_Nph_Theta_cer") fStree->Draw("Nph2_cer_L/Nph_cer_L:Theta>>h_Nph2_Over_Nph_Theta_cer",drawcut,"");
1176	else if(option == "Ne2_Over_Ne_Hmax") fStree->Draw("Ne2/Ne:Hmax_show>>h_Ne2_Over_Ne_Hmax",drawcut,"");
1177	else if(option == "Ne2_Over_Ne_Theta") fStree->Draw("Ne2/Ne:Theta>>h_Ne2_Over_Ne_Theta",drawcut,"");
1178	else if(option == "Width2_Over_Width_Hmax_fluo") fStree->Draw("Width2_fluo/Width_fluo:Hmax_show>>h_Width2_Over_Width_Hmax_fluo",drawcut,"");
1179	else if(option == "Width2_Over_Width_Hmax_shower") fStree->Draw("Width2_show/Width_show:Hmax_show>>h_Width2_Over_Width_Hmax_shower",drawcut,"");
1180	else if(option == "Width2_Over_Width_Theta_fluo") fStree->Draw("Width2_fluo/Width_fluo:Theta>>h_Width2_Over_Width_Theta_fluo",drawcut,"");
1181
1182	// shower correlations
1183	// hmax -- Xmax
1184	else if(option == "Hmax_Theta") fStree->Draw("Hmax_fluo:Theta>>h_Hmax_Theta",drawcut,"");
1185	else if(option == "Xmax_Energy") fStree->Draw("TrueXmax:log10(Energy)>>h_Xmax_Energy",drawcut,"");
1186	// energy
1187
1188	else if(option == "Nemax_Energy") fStree->Draw("log10(Nemax):log10(Energy)>>h_Nemax_Energy",drawcut,"profs");
1189	else if(option == "Nemax_Energy_linear") fStree->Draw("Nemax:log10(Energy)>>h_Nemax_Energy_linear",drawcut,"goffprofs");
1190	else if(option == "Ne2_Energy") fStree->Draw("log10(Ne2):log10(Energy)>>h_Ne2_Energy",drawcut,"profs");
1191	else if(option == "Ne2_Energy_linear") fStree->Draw("Ne2:log10(Energy)>>h_Ne2_Energy_linear",drawcut,"goffprofs");
1192	// theta
1193	else if(option == "Ne_Theta") fStree->Draw("Ne/Energy:Theta>>h_Ne_Theta",drawcut,"");
1194	else if(option == "Ne2_Theta") fStree->Draw("Ne2/Energy:Theta>>h_Ne2_Theta",drawcut,"");
1195	else if(option == "Nemax_Theta") fStree->Draw("Nemax/Energy:Theta>>h_Nemax_Theta",drawcut,"");
1196	else if(option == "Widthshow_Hmax") fStree->Draw("Width2_show/Width_show:Hmax_show>>h_Widthshow_Hmax",drawcut,"");
1197
1198	// photons in atmosphere and shower correlations
1199	else if(option == "Nphmaxfluo_Nemax") fStree->Draw("log10(Nmax_fluo_L):log10(Nemax)>>h_Nphmaxfluo_Nemax",drawcut,"profs");
1200	else if(option == "Nphmaxfluo_Nemax_linear") fStree->Draw("Nmax_fluo_Lcos(ThetaTMath::Pi()/180):log10(Nemax)>>h_Nphmaxfluo_Nemax_linear",drawcut,"profs");
1201	else if(option == "Nph2fluo_Ne2") fStree->Draw("log10(Nph2_fluo_L):log10(Ne2)>>h_Nph2fluo_Ne2",drawcut,"profs");
1202	else if(option == "Nph2fluo_Ne2_linear") fStree->Draw("Nph2_fluo_Lcos(ThetaTMath::Pi()/180):log10(Ne2)>>h_Nph2fluo_Ne2_linear",drawcut,"profs");
1203	else if(option == "Nph2ckov_Ne2") fStree->Draw("log10(Nph2_cer_L):log10(Ne2)>>h_Nph2ckov_Ne2",drawcut,"profs");
1204	else if(option == "Nph2ckov_Ne2_linear") fStree->Draw("Nph2_cer_L:log10(Ne2)>>h_Nph2ckov_Ne2_linear",drawcut,"profs");
1205	else if(option == "Nph2fluo_Nemax") fStree->Draw("log10(Nph2_fluo_L):log10(Nemax)>>h_Nph2fluo_Nemax",drawcut,"profs");
1206	else if(option == "Nph2fluo_Nemax_linear") fStree->Draw("Nph2_fluo_Lcos(ThetaTMath::Pi()/180):log10(Nemax)>>h_Nph2fluo_Nemax_linear",drawcut,"profs");
1207
1208	else if(option == "Width2_fluo_Hmax_fluo") fStree->Draw("Width2_fluo:Hmax_fluo>>h_Width2_fluo_Hmax_fluo",drawcut,"");
1209	else if(option == "Width2_fluo_Theta") fStree->Draw("Width2_fluo:Theta>>h_Width2_fluo_Theta",drawcut,"");
1210
1211	else if(option == "Xmaxfluo_Xmax_Theta") fStree->Draw("Xmaxph_fluo - Xmax_show:Theta>>h_Xmaxfluo_Xmax_Theta",drawcut,"");
1212	else if(option == "Xmaxfluo_Xmax_Hmax") fStree->Draw("Xmaxph_fluo - Xmax_show:Hmax_show>>h_Xmaxfluo_Xmax_Hmax",drawcut,"");
1213	else if(option == "Xmaxfluo_Xmax_Energy") fStree->Draw("Xmaxph_fluo - Xmax_show:log10(Energy)>>h_Xmaxfluo_Xmax_Energy",drawcut,"");
1214	else if(option == "Hmaxfluo_Hmax_Theta") fStree->Draw("Hmax_fluo - Hmax_show:Theta>>h_Hmaxfluo_Hmax_Theta",drawcut,"");
1215	else if(option == "Hmaxcer_Hmax_Theta") fStree->Draw("Hmax_cer - Hmax_show:Theta>>h_Hmaxcer_Hmax_Theta",drawcut,"");
1216	else if(option == "Xmaxfluo_Xmax_X1") fStree->Draw("Xmaxph_fluo - Xmax_show:X1>>h_Xmaxfluo_Xmax_X1",drawcut,"");
1217	else if(option == "Lmaxfluo_Lmax_Theta") fStree->Draw("Lmax_fluo - Lmax_show:Theta>>h_Lmaxfluo_Lmax_Theta",drawcut,"");
1218	else if(option == "Rmaxfluo_Rmax_Theta") fStree->Draw("sqrt((ymax_fluo - TrueMaxPosY)(ymax_fluo - TrueMaxPosY) +(xmax_fluo - TrueMaxPosX)(xmax_fluo - TrueMaxPosX)):Theta>>h_Rmaxfluo_Rmax_Theta",drawcut,"");
1219	else if(option == "xmax_ymax_shift") fStree->Draw("Hmax_show:ymax_fluo - TrueMaxPosY:xmax_fluo - TrueMaxPosX>>h_xmax_ymax_shift",drawcut,"");
1220
1221	// photons in atmosphere correlations VS Energy
1222	else if(option == "Nmax_Energy_fluo") fStree->Draw("log10(Nmax_fluo_Lcos(ThetaTMath::Pi()/180)):log10(Energy)>>h_Nmax_Energy_fluo",drawcut,"profs");
1223	else if(option == "Nmax_Energy_cer") fStree->Draw("log10(Nmax_cer):log10(Energy)>>h_Nmax_Energy_cer",drawcut,"profs");
1224	else if(option == "Nph2_Energy_fluo") fStree->Draw("log10(Nph2_fluo_Lcos(ThetaTMath::Pi()/180)):log10(Energy)>>h_Nph2_Energy_fluo",drawcut,"profs");
1225	else if(option == "Nph2_Energy_cer") fStree->Draw("log10(Nph2_cer_L):log10(Energy)>>h_Nph2_Energy_cer",drawcut,"profs");
1226	else if(option == "Nmax_Energy_fluo_linear") fStree->Draw("Nmax_fluo_Lcos(ThetaTMath::Pi()/180):log10(Energy)>>h_Nmax_Energy_fluo_linear",drawcut,"goffprofs");
1227	else if(option == "Nmax_Energy_cer_linear") fStree->Draw("Nmax_cer:log10(Energy)>>h_Nmax_Energy_cer_linear",drawcut,"goffprofs");
1228	else if(option == "Nph2_Energy_fluo_linear") fStree->Draw("Nph2_fluo_Lcos(ThetaTMath::Pi()/180):log10(Energy)>>h_Nph2_Energy_fluo_linear",drawcut,"goffprofs");
1229	else if(option == "Nph2_Energy_cer_linear") fStree->Draw("Nph2_cer_L:log10(Energy)>>h_Nph2_Energy_cer_linear",drawcut,"goffprofs");
1230
1231	// photons in atmosphere correlations VS theta
1232	else if(option == "Nph2_Theta_fluo") fStree->Draw("Nph2_fluo_L/Energy:Theta>>h_Nph2_Theta_fluo",drawcut,"");
1233	else if(option == "Nph2_Theta_cer") fStree->Draw("Nph2_cer_L/Energy:Theta>>h_Nph2_Theta_cer",drawcut,"");
1234	else if(option == "Nmax_Hmax_fluo") fStree->Draw("Nmax_fluo_L/Energy:Hmax_fluo>>h_Nmax_Hmax_fluo",drawcut,"");
1235	else if(option == "Nmax_Nph2_fluo_Hmax") fStree->Draw("Nmax_fluo_L*1000/Nph2_fluo_L:Hmax_fluo>>h_Nmax_Nph2_fluo_Hmax",drawcut,"");
1236	else if(option == "Nmax_Theta_cer") fStree->Draw("Nmax_cer/Energy:Theta>>h_Nmax_Theta_cer",drawcut,"");
1237	else if(option == "Nph2_Hmax_cer") fStree->Draw("Nph2_cer_L/Energy:Hmax_show>>h_Nph2_Hmax_cer",drawcut,"");
1238	else if(option == "Nph2_Hmax_fluo") fStree->Draw("Nph2_fluo_L/Energy:Hmax_show>>h_Nph2_Hmax_fluo",drawcut,"");
1239	else if(option == "Yield_max_Hmax") fStree->Draw("Yieldmax_fluo:Hmax_show>>h_Yield_max_Hmax",drawcut,"");
1240	else if(option == "Yield_max_cer_Hmax") fStree->Draw("Yieldmax_ckov:Hmax_show>>h_Yield_max_cer_Hmax",drawcut,"");
1241	else if(option == "Yield_mean_Hmax") fStree->Draw("Yield_f:Hmax_show>>h_Yield_mean_Hmax",drawcut,"");
1242
1243	// photons on pupil correlations VS Energy
1244	else if(option == "Nph_Energy_CB_P") fStree->Draw("log10(Nph_cer_air_pupilcos(ThetaTMath::Pi()/180)):log10(Energy)>>h_Nph_Energy_CB_P",drawcut,"profs");
1245	else if(option == "Nph_Energy_CR_P") fStree->Draw("log10(Nph_cer_refl_pupil):log10(Energy)>>h_Nph_Energy_CR_P",drawcut,"profs");
1246	else if(option == "Nph2_Energy_Fluo_P") fStree->Draw("log10(Nph2_fluo_pupilcos(ThetaTMath::Pi()/180)):log10(Energy)>>h_Nph2_Energy_Fluo_P",drawcut,"profs");
1247	else if(option == "Nmax_Energy_Fluo_P") fStree->Draw("log10(Nmax_fluo_pupilcos(ThetaTMath::Pi()/180)):log10(Energy)>>h_Nmax_Energy_Fluo_P",drawcut,"profs");
1248	else if(option == "Nph_Energy_CB_P_linear") fStree->Draw("Nph_cer_air_pupilcos(ThetaTMath::Pi()/180):log10(Energy)>>h_Nph_Energy_CB_P_linear",drawcut,"goffprofs");
1249	else if(option == "Nph_Energy_CR_P_linear") fStree->Draw("Nph_cer_refl_pupil:log10(Energy)>>h_Nph_Energy_CR_P_linear",drawcut,"goffprofs");
1250	else if(option == "Nph2_Energy_Fluo_P_linear") fStree->Draw("Nph2_fluo_pupilcos(ThetaTMath::Pi()/180):log10(Energy)>>h_Nph2_Energy_Fluo_P_linear",drawcut,"goffprofs");
1251	else if(option == "Nmax_Energy_Fluo_P_linear") fStree->Draw("Nmax_fluo_pupilcos(ThetaTMath::Pi()/180):log10(Energy)>>h_Nmax_Energy_Fluo_P_linear",drawcut,"goffprofs");
1252
1253	// photons on pupil correlations VS theta
1254	//else if(option == "Nph_Theta_CB_P") fStree->Draw("Nph_cer_back_pupil/Energy:Theta>>h_Nph_Theta_CB_P",drawcut,"");
1255	else if(option == "Nph_Theta_CB_P") fStree->Draw("Nph_cer_air_pupil/Energy:Theta>>h_Nph_Theta_CB_P",drawcut,"");
1256	else if(option == "Nph_Theta_CR_P") fStree->Draw("Nph_cer_refl_pupil/Energy:Theta>>h_Nph_Theta_CR_P",drawcut,"");
1257	else if(option == "Nph_Theta_Fluo2_P") fStree->Draw("Nph2_fluo_pupil/Energy:Theta>>h_Nph_Theta_Fluo2_P",drawcut,"");
1258	else if(option == "Nmax_fluo_Theta_P") fStree->Draw("Nmax_fluo_pupil/Energy:Theta>>h_Nmax_fluo_Theta_P",drawcut,"");
1259	else if(option == "Nph_Theta_Fluo_P_crash") fStree->Draw("Nph2_fluo_pupil/Nph_fluo_pupil:Theta>>h_Nph_Theta_Fluo_P_crash",drawcut,"");
1260	//else if(option == "Nph_airscat_effect") fStree->Draw("(Nph_fluo_pupil+Nph_cer_back_pupil+Nph_cer_refl_pupil)/(Nph_fluo_pupil+Nph_cer_refl_pupil):Theta>>h_Nph_airscat_effect",drawcut,"");
1261	else if(option == "Nph_airscat_effect") fStree->Draw("(Nph_fluo_pupil+Nph_cer_air_pupil+Nph_cer_refl_pupil)/(Nph_fluo_pupil+Nph_cer_refl_pupil):Theta>>h_Nph_airscat_effect",drawcut,"");
1262
1263	else if(option == "TimeWidth2_Theta") fStree->Draw("Width_fluo_pupil:Theta>>h_TimeWidth2_Theta",drawcut,"");
1264
1265	// Transmissions fluo (+ckov) : ratios (inomega / onpupil)
1266	else if(option == "Trans_fluo_Theta") fStree->Draw("Nph_fluo_pupil/FluoInOmega:Theta>>h_Trans_fluo_Theta",drawcut,"");
1267	else if(option == "Trans_max_mean") fStree->Draw("(Trans_fluo_atMax - (FluoInOmega/Nph_fluo_pupil))/(FluoInOmega/Nph_fluo_pupil):Theta>>h_Trans_max_mean",drawcut,"");
1268	else if(option == "Trans_GTUmax_max") fStree->Draw("(Trans_fluo_GTUmax - Trans_fluo_atMax)/Trans_fluo_atMax>>h_Trans_GTUmax_max",drawcut,"");
1269	//else if(option == "Trans_ckov_Theta") fStree->Draw("(Nph_cer_back_pupil+Nph_cer_refl_pupil)/CkovInOmega:Theta>>h_Trans_ckov_Theta",drawcut,"");
1270	else if(option == "Trans_ckov_Theta") fStree->Draw("(Nph_cer_air_pupil+Nph_cer_refl_pupil)/CkovInOmega:Theta>>h_Trans_ckov_Theta",drawcut,"");
1271
1272	// MISCELLANEOUS//////////////////////////
1273
1274	// fluo yield properties
1275	else if(option == "Yield_GTUmax_max") fStree->Draw("(Yield_fluoGTUmax - Yieldmax_fluo)/Yieldmax_fluo:Theta>>h_Yield_GTUmax_max",drawcut,"");
1276
1277	// XGTUmax shift due to scattered ckov
1278	else if(option == "Xmax_shift") fStree->Draw("(X_GTUmax_tot - X_GTUmax)/X_GTUmax:Theta>>h_Xmax_shift",drawcut,"");
1279
1280	// correlations between width of time profile on pupil VS hmax
1281	else if(option == "TimeWidth2_Hmax") fStree->Draw("Width_fluo_pupil:Hmax_show>>h_TimeWidth2_Hmax",drawcut,"");
1282
1283
1284	else Printf("<ETreePainter::DrawHistos()> : %s OPTION not allowed", option.Data());
1285	}
1286
1287	//_____________________________________________________________________________
1288	TH1F* ETreePainter::BuildAltHisto( Option_t *opt ) {
1289	//
1290	// generate histos as function of altitude (Xaxis in km)
1291	//
1292
1293	TString option(opt);
1294
1295	// if no data available
1296	if ( !fAtmo \|\| !fShower) {
1297	Info("BuildAltHistos()","EAtmosphere (or EShower) object is NULL. Painter made zombie.");
1298	MakeZombie();
1299	return 0;
1300	}
1301
1302	Int_t n = fAtmo->GetNumBunches();
1303
1304	// if photons data are empty
1305	if ( n <= 0 ) {
1306	Info("BuildAltHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
1307	MakeZombie();
1308	return 0;
1309	}
1310
1311	// determine bins for development w.r.t altitude
1312	Double_t temp[n+1];
1313	Int_t Nbins = 0;
1314	for (Int_t i=0; i<n; i++) {
1315	if (fAtmo->GetBunch(i)->GetType() == 0 ) {
1316	temp[Nbins] = Zv(fAtmo->GetBunch(i)->GetShowerPosi())/km;
1317	temp[Nbins+1] = Zv(fAtmo->GetBunch(i)->GetShowerPosf())/km;
1318	Nbins++;
1319	}
1320	}
1321
1322	// fill X-axis array for bunches altitude
1323	Double_t AltBins[Nbins+1];
1324	Int_t incrm = 0;
1325	for (Int_t i=0; i<Nbins+1; i++) {
1326	AltBins[Nbins - incrm] = temp[i];
1327	incrm++;
1328	}
1329
1330
1331
1332	TH1F* h = 0;
1333
1334	// fill histograms from bunches
1335	EShowerStep* shstep = 0;
1336	Int_t bin(-1);
1337	for (Int_t i=0; i<n; i++) {
1338	if(i%2 == 0) shstep = fShower->GetStep(i/2);
1339	Double_t alt = 0.5*( Zv(fAtmo->GetBunch(i)->GetShowerPosi())
1340	+ Zv(fAtmo->GetBunch(i)->GetShowerPosf()) )/km;
1341	Double_t Nph = fAtmo->GetBunch(i)->GetWeight();
1342	Double_t Nelec = shstep->GetNumElectrons();
1343
1344	if (fAtmo->GetBunch(i)->GetType()==0) {
1345	if(option == "alt_Bunch_Nph_F") {
1346	bin = h->FindBin(alt);
1347	h->SetBinContent(bin,Nph/h->GetBinWidth(bin));
1348	}
1349	if(option == "alt_Shower") {
1350	bin = h->FindBin(alt);
1351	h->SetBinContent(bin,Nelec);
1352	}
1353
1354	}
1355	if (fAtmo->GetBunch(i)->GetType()==1) {
1356	if(option =="alt_Bunch_Nph_C" ) {
1357	bin = h->FindBin(alt);
1358	h->SetBinContent(bin,Nph/h->GetBinWidth(bin));
1359	}
1360	}
1361	}
1362
1363	return h;
1364	}
1365
1366	//_____________________________________________________________________________
1367	TH1F* ETreePainter::BuildLongitudinalHisto( Option_t* opt ) {
1368	//
1369	// generate histos along the track (Longitudinal Profiles)
1370	// (Xaxis in km or in g/cm2)
1371	//
1372
1373	TString option(opt);
1374
1375	// if no data available
1376	if ( !fAtmo \|\| !fShower) {
1377	Info("BuildLongitudinalHistos()","EAtmosphere or EShower (needed for this method) object is NULL. Painter made zombie.");
1378	MakeZombie();
1379	return 0;
1380	}
1381
1382	Int_t n = fAtmo->GetNumBunches();
1383	Int_t nsh = fShower->GetNumSteps();
1384	if(n != 2*nsh) Warning("BuildLongitudinalHistos()","Nb of bunches and Nb of showersteps SHOULD BE THE SAME");
1385
1386	// if photons data are empty
1387	if ( n <= 0 ) {
1388	Info("BuildLongitudinalHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
1389	MakeZombie();
1390	return 0;
1391	}
1392
1393	// determine bins for longitudinal development (both in km and g/cm2)
1394	UInt_t Nbins = UInt_t(nsh);
1395	Double_t LongBins[Nbins+1];
1396	Double_t GramBins[Nbins+1];
1397	TVector3 initpos = fShower->GetStep(0)->GetPosi();
1398	EShowerStep* shstep = 0;
1399	for (Int_t i=0; i<nsh; i++) {
1400	shstep = fShower->GetStep(i);
1401	LongBins[i] = (shstep->GetPosi() - initpos).Mag()/km;
1402	LongBins[i+1] = (shstep->GetPosf() - initpos).Mag()/km;
1403	GramBins[i] = shstep->GetXi()*cm2/g;
1404	GramBins[i+1] = shstep->GetXf()*cm2/g;
1405	}
1406
1407	TH1F* h = 0;
1408	if( option == "Fluo_Longit_prof_km") h = new TH1F("Fluo_Longit_prof_km","Fluo_Longitudinal_profile_km",Nbins,LongBins);
1409	else if(option == "Fluo_Longit_prof_gram") h = new TH1F("Fluo_Longit_prof_gram","Fluo_Longit_prof_gram",Nbins,GramBins);
1410	else if(option == "Ckov_Longit_prof_km") h = new TH1F("Ckov_Longit_prof_km","Ckov_Longit_prof_km",Nbins,LongBins);
1411	else if(option == "Ckov_Longit_prof_gram" ) h = new TH1F("Ckov_Longit_prof_gram","Ckov_Longit_prof_gram",Nbins,GramBins);
1412	else if(option == "Shower_Longit_prof_gram" ) h = new TH1F("Shower_Longit_prof_gram","Shower_Longit_prof_gram",Nbins,GramBins);
1413	else if(option == "Shower_Longit_prof_km" ) h = new TH1F("Shower_Longit_prof_km","Shower_Longit_prof_km",Nbins,LongBins);
1414
1415
1416	// fill histograms from fluo bunches
1417	// photons density is plotted (per km --or-- per g/cm2)
1418	Int_t bin(-1);
1419	for (Int_t i=0; i<n; i++) {
1420	if(i%2 == 0) shstep = fShower->GetStep(i/2);
1421	Double_t dist = ((0.5*(shstep->GetPosf() + shstep->GetPosi())) - initpos).Mag()/km;
1422	Double_t grammage = (0.5(shstep->GetXf() + shstep->GetXi()))cm2/g;
1423	Double_t nph = fAtmo->GetBunch(i)->GetWeight();
1424	Double_t nelec = shstep->GetNumElectrons();
1425
1426	if (fAtmo->GetBunch(i)->GetType()==0) {
1427	if(option == "Fluo_Longit_prof_km") {
1428	bin = h->FindBin(dist);
1429	h->SetBinContent(bin,nph/h->GetBinWidth(bin));
1430	}
1431	if(option == "Fluo_Longit_prof_gram") {
1432	bin = h->FindBin(grammage);
1433	h->SetBinContent(bin,nph/h->GetBinWidth(bin));
1434	}
1435	if(option == "Shower_Longit_prof_gram") {
1436	bin = h->FindBin(grammage);
1437	h->SetBinContent(bin,nelec);
1438	}
1439	if(option == "Shower_Longit_prof_km") {
1440	bin = h->FindBin(dist);
1441	h->SetBinContent(bin,nelec);
1442	}
1443	}
1444	if (fAtmo->GetBunch(i)->GetType()==1) {
1445	if(option == "Ckov_Longit_prof_km") {
1446	bin = h->FindBin(dist);
1447	h->SetBinContent(bin,nph/h->GetBinWidth(bin));
1448	}
1449	if(option =="Ckov_Longit_prof_gram" ) {
1450	bin = h->FindBin(grammage);
1451	h->SetBinContent(bin,nph/h->GetBinWidth(bin));
1452	}
1453	}
1454	}
1455
1456	return h;
1457	}
1458
1459	//_____________________________________________________________________________
1460	TH1F* ETreePainter::BuildTofHisto( Option_t* opt ) {
1461	//
1462	// generate histos of photons on pupil as function of time since PRIMARY cosmic ray first interaction
1463	// (Xaxis in microsecond -> binwidth in GTU !!)
1464	//
1465
1466	TString option(opt);
1467
1468	// if no data available
1469	if ( !fAtmo ) {
1470	Info("BuildTofHistos()","EAtmosphere object is NULL. Painter made zombie.");
1471	MakeZombie();
1472	return 0;
1473	}
1474
1475	Int_t ns = fAtmo->GetNumSingles();
1476
1477	// if photons data are empty
1478	if ( ns <= 0 ) {
1479	//Info("BuildTofHisto()","No SinglePhoton generated in Atmosphere");
1480	}
1481
1482	// determine bins for time on pupil
1483	Double_t TimeMin = 3000.*microsecond;
1484	Double_t TimeMax = 0;
1485	for (Int_t i=0; i<ns; i++) {
1486	if ( !(fAtmo->GetSingle(i)->IsAbsorbed()) ) {
1487	if ( TimeMin > (fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate()) )
1488	TimeMin = fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate();
1489	if ( TimeMax < (fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate()) )
1490	TimeMax = fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate();
1491	}
1492	}
1493
1494	// build histogram vs time on pupil
1495	Float_t min = (Int_t)(TimeMin/microsecond)-1;
1496	Int_t nbins = (Int_t)((TimeMax-TimeMin)/(fGTU)) + 1;
1497	TH1F* h = 0;
1498	if(option == "Single_Nph_F_P_t") h = new TH1F("Single_Nph_F_P_t","Nph_Fluo_on_Pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1499	else if(option == "Single_Nph_Cdir_P_t") h = new TH1F("Single_Nph_CR_P_t","Nph_Cerenkov direct on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1500	else if(option == "Single_Nph_CB_P_t") h = new TH1F("Single_Nph_CB_P_t","Nph_Cerenkov back on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1501	else if(option == "Single_Nph_CR_P_t") h = new TH1F("Single_Nph_CR_P_t","Nph_Cerenkov refl on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1502	else if(option == "Single_Nph_CB_cloud_P_t") h = new TH1F("Single_Nph_CB_cloud_P_t","Nph_Cerenkov clouds scat on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1503	else if(option == "Single_Nph_F-CB_P_t") h = new TH1F("Single_Nph_F-CB_P_t","Nph_fluo_ckovBack on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
1504
1505	// fill histograms from singles
1506	for (Int_t i=0; i<ns; i++) {
1507	Double_t tof = (fAtmo->GetSingle(i)->GetTof() + fAtmo->GetSingle(i)->GetDate()) / microsecond;
1508	if ( !(fAtmo->GetSingle(i)->IsAbsorbed()) ) {
1509	if ( fAtmo->GetSingle(i)->GetType() == 0 ) {
1510	if(option == "Single_Nph_F_P_t") h->Fill(tof,1);
1511	if(option == "Single_Nph_F-CB_P_t") h->Fill(tof,1);
1512	}
1513	if ( fAtmo->GetSingle(i)->GetType() == 1 ) {
1514	if ( fAtmo->GetSingle(i)->GetHistory()==0) {
1515	if(option == "Single_Nph_Cdir_P_t") h->Fill(tof,1);
1516	}
1517	if ( fAtmo->GetSingle(i)->GetHistory()==2) {
1518	if(option == "Single_Nph_CB_P_t") h->Fill(tof,1);
1519	if(option == "Single_Nph_F-CB_P_t") h->Fill(tof,1);
1520	}
1521	if ( fAtmo->GetSingle(i)->GetHistory()==3) {
1522	if(option == "Single_Nph_CB_cloud_P_t") h->Fill(tof,1);
1523	}
1524	if ( fAtmo->GetSingle(i)->GetHistory()==1 )
1525	if(option == "Single_Nph_CR_P_t") h->Fill(tof,1);
1526	}
1527	}
1528	}
1529
1530	return h;
1531	}
1532
1533	//_____________________________________________________________________________
1534	TH1F* ETreePainter::BuildWlHisto( Option_t* opt ) {
1535	//
1536	// generate wavelength spectra of photons on pupil BEFORE and AFTER last transmission
1537	// generate also wavelength spectra of photons in atmosphere at creation
1538	// (Xaxis in nm)
1539	//
1540
1541	TString option(opt);
1542
1543	// if no data available
1544	if ( !fAtmo ) {
1545	Info("BuildWlHistos()","EAtmosphere object is NULL. Painter made zombie.");
1546	MakeZombie();
1547	return 0;
1548	}
1549
1550	Int_t n = fAtmo->GetNumBunches();
1551	Int_t ns = fAtmo->GetNumSingles();
1552
1553	// if photons data are empty
1554	if ( n <= 0 ) {
1555	Info("BuildWlHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
1556	MakeZombie();
1557	return 0;
1558	}
1559
1560	TH1F* h = 0;
1561
1562	// bunches at creation spectra
1563	//else if(opt == "Bunch_Wl_tot") h = new TH1F("Bunch_Wl_tot", "Nph_tot vs Lambda",150,300,450);
1564
1565	// photons on pupil spectra
1566	if(option == "Single_Wl_F_P") h = new TH1F("Single_Wl_F_P", "Nph_Fluo_on_pupil vs Lambda",150,300,450);
1567	else if(option == "Single_Wl_CB_P") h = new TH1F("Single_Wl_CB_P", "Nph_Cerenkov_back_on_pupil vs Lambda",150,300,450);
1568	else if(option == "Single_Wl_CR_P") h = new TH1F("Single_Wl_CR_P","Nph_Cerenkov refl on pupil vs lambda",150,300,450);
1569	else if(option == "Single_Wl_tot_P") h = new TH1F("Single_Wl_tot_P","Nph_tot on pupil vs lambda",150,300,450);
1570	//TH1F *th16_2 = GetWlHisto("Single_Wl_ckov_CloudBackscat_P", "Nph_Cerenkov_Cloudback_on_pupil vs Lambda");
1571	//TH1F *th16_tot = GetWlHisto("Single_Wl_CB_tot_P", "Nph_Cerenkov_Total_back_on_pupil vs Lambda"); // not used so far
1572
1573
1574	// fill histograms from bunches
1575	for (Int_t i=0; i<n; i++) {
1576	Double_t Nph = fAtmo->GetBunch(i)->GetWeight();
1577	Int_t NumWl = fAtmo->GetBunch(i)->GetNumWavelengths();
1578	const Float_t* lambda = 0;
1579	const Float_t* wlweight = 0;
1580	wlweight = fAtmo->GetBunch(i)->GetTable("weight");
1581	lambda = fAtmo->GetBunch(i)->GetTable("lambda"); // bin centers..
1582	TArrayF Xarray(NumWl+1); // bin low edge
1583	Float_t binw = (lambda[NumWl-1] - lambda[0]) / Float_t(NumWl);
1584	Xarray[0] = (lambda[0] - 0.5*binw) / nm;
1585	for(Int_t m=0; m<NumWl; m++) Xarray[m+1] = (lambda[m] + 0.5*binw) / nm;
1586
1587	if (fAtmo->GetBunch(i)->GetType()==0 && option == "Bunch_Wl_F") {
1588	// fluorescence bunches
1589	if(!h && NumWl>2) h = new TH1F("Bunch_Wl_F","Nph_F vs Lambda",Xarray.GetSize() - 1,Xarray.GetArray()); // NumWl to avoid fake bunches (see ShowerLS.cc)
1590	if(h) for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph);
1591	}
1592
1593	else if (fAtmo->GetBunch(i)->GetType()==1 && option == "Bunch_Wl_C") {
1594	// Cerenkov bunches
1595	if(!h) h = new TH1F("Bunch_Wl_C","Nph_C vs Lambda",Xarray.GetSize() - 1,Xarray.GetArray());
1596	for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph);
1597	}
1598
1599	// problematic for Xaxis definition + not really usefull so far
1600	//if(opt == "Bunch_Wl_tot") {
1601	// for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph );
1602	//}
1603
1604	}
1605
1606	// fill histograms from singles
1607	for (Int_t i=0; i<ns; i++) {
1608	Double_t wl = fAtmo->GetSingle(i)->GetWl()/nm;
1609	if (!(fAtmo->GetSingle(i)->IsAbsorbed())) {
1610
1611	// all photons
1612	if(option == "Single_Wl_tot_P") h->Fill(wl,1);
1613
1614	// for fluorescence photons
1615	if ( fAtmo->GetSingle(i)->GetType() == 0 && option == "Single_Wl_F_P") {
1616	h->Fill(wl,1);
1617	}
1618
1619	// for Cerenkov photons
1620	if ( fAtmo->GetSingle(i)->GetType() == 1 ) {
1621	// for backscattered ones
1622	if ( fAtmo->GetSingle(i)->GetHistory()==2 && option == "Single_Wl_CB_P") {
1623	h->Fill(wl,1);
1624	}
1625	// for reflected ones
1626	if ( fAtmo->GetSingle(i)->GetHistory()==1 && option == "Single_Wl_CR_P") {
1627	h->Fill(wl,1);
1628	}
1629	}
1630	}
1631	}
1632
1633	return h;
1634	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: JEM-EUSO/esaf_cc_at_lal/packages/common/eventviewer/src/ETreePainter.cc @ 114

Download in other formats: