Context Navigation

TrackDirection2Module.cc @ 181

Last change on this file since 181 was 181, checked in by moretto, 11 years ago
changes from Alejandro r:3036
File size: 96.7 KB

Line
1	// $Id: TrackDirection2Module.cc 3036 2013-07-08 09:25:47Z guzman $
2	// Author: elena Nov, 19 2004
3
4	/*****************************************************************************
5	* ESAF: Euso Simulation and Analysis Framework *
6	* *
7	* Id: TrackDirection2Module *
8	* Package: Fitting *
9	* Coordinator: <coordinator> *
10	* *
11	*****************************************************************************/
12
13	//_____________________________________________________________________________
14	//
15	// TrackDirection2Module
16	//
17	// This module is devoted to the reconstruction of the shower direction.
18	// It implements some different algorithms (descibed in detail in ....)
19	// Essentially the module keeps the points found by pattern recogniton
20	// (clustering or Hough transform) and first find the plane that contains the
21	// track and the detector (TDP).
22	// The reserarch of TDP can be done in the following ways:
23	// - further selection of points with hough transform
24	// - further selection of points with shape selection method
25	// - no further selection
26	// In each case the TDP is founded by a fit (least squares, median or hough)
27	// of the x-t, y-t projections of points on the plane z=0.
28	//
29	// Then the shower direction is reconstructed by one of the following methods:
30	// - analytical approximated 1 AA1()
31	// - analytical approximated 2 AA2()
32	// - numerical exact 1 NE1()
33	// - numerical exact 2 NE2()
34	// - analytical exact 1 AE1()
35	//
36	// The AA1() method is in each case used in order to initialize the fit for all
37	// other methods.
38	//
39	// Config file parameters
40	// ======================
41	//
42	// fNumHitsMinimum : minimum number of hits per pixel in a GTU
43	// if >0 the threshold is absolute
44	// if <0 this threshold is relative to the background mean B
45	// true threshold = B + \|fNumHitsMinimum\|*sqrt(B)
46	// fNumPointsMinimum : minimum number of points to proceed with reconstruction
47	// fSignalToNoise : take pixels with signal/sqrt(noise) >= fSignalToNoise number
48	// [valid for ProcessOnlySignal option enabled in RootInputModule]
49	// fUseHough [bool] : use hough transform to find the TDP
50	// fUseHoughwithselection [bool] : use the points selected with Hough transform
51	// in the reconstruction of direction
52	// fUseShape [bool] : use shape selection method to find the TDP
53	// fUseShapeInModules [bool] : use the points selected with shape selection
54	// in the reconstruction of direction
55	// fDebugInfo [bool] : compute and display some debug informations
56	//
57	// fMulti1 : multiplier for shape selection method
58	// fMulti2 : multiplier for shape selection method
59	//
60	// fAA1FitMethod : fit method for AA1 algorithm
61	// - Valid options : linear (least squares fit)
62	// median (median fit)
63	// hough (hough fit)
64	//
65	// fMethod : direction reconstruction method
66	// - Valid options : AA1 \|\| AA2 \|\| NE1 \|\| NE2 \|\| AE1 (single algorithm)
67	// all (executes all algorithms)
68	//
69	// fFixTmaxNumeric [bool] : fix shower maximum parameters in numerical fits
70	// fErrAngle : angular error [deg]
71	//
72	// fStat1 : angular error [deg] value 1 for statistics
73	// fStat2 : angular error [deg] value 2 for statistics
74	// - Statistics are made between 0 < err < fStat1 and fStat1 < err < fStat2
75	//
76	// fDoGraphUseShape [bool] : save some debug graph in rootfile
77	//
78	// fMinuitOutputLevel : set the MINUIT output display level
79	//
80
81	#include "EusoCluster.hh"
82	#include "TrackDirection2Module.hh"
83	#include "RecoEvent.hh"
84	#include "RecoGlobalData.hh"
85	#include "RecoPixelData.hh"
86	#include "LeastSquaresFit.hh"
87	#include "MedianFit.hh"
88	#include "HoughFit.hh"
89	#include "RecoRootEvent.hh"
90	#include "EConst.hh"
91	#include "Config.hh"
92	#include "ERunParameters.hh"
93
94	#include <TGraph.h>
95	#include <TCanvas.h>
96	#include <TMinuit.h>
97	#include <TVector3.h>
98	#include <TH2F.h>
99	#include <TLegend.h>
100	#include <TLinearFitter.h>
101	#include <fstream>
102	#include <TDirectory.h>
103
104	using namespace TMath;
105	using namespace sou;
106	using namespace EConst;
107
108	//FCN function called by numeric method NE1
109	void ChiSquareNE1(Int_t &npar, Double_t gin, Double_t &chisqnorm, Double_t par, Int_t iflag);
110
111	//FCN function called by numeric method NE2
112	void ChiSquareNE2(Int_t &npar, Double_t gin, Double_t &chisqnorm, Double_t par, Int_t iflag);
113
114	ClassImp(TrackDirection2Module)
115	ClassImp(ContainerData)
116
117	//_____________________________________________________________________________
118	TrackDirection2Module::TrackDirection2Module() : RecoModule("TrackDirection2") {
119	//
120	// ctor
121	//
122	}
123
124	//_____________________________________________________________________________
125	TrackDirection2Module::~TrackDirection2Module() {
126	//
127	// dtor
128	//
129	}
130
131	//_____________________________________________________________________________
132	Bool_t TrackDirection2Module::Init() {
133	//
134	// Initialization of variables
135	//
136
137	Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
138
139	useLFaa1 = kFALSE; useMFaa1 = kFALSE; useHFaa1 = kFALSE;
140	useLFplane = kFALSE; useMFplane = kFALSE; useHFplane = kFALSE;
141
142	fSignalToNoise = (Double_t)Conf()->GetNum("TrackDirection2Module.fSignalToNoise");
143
144	fStat1 = Conf()->GetNum("TrackDirection2Module.fStat1");
145	fStat2 = Conf()->GetNum("TrackDirection2Module.fStat2");
146	fDoGraphUseShape = Conf()->GetBool("TrackDirection2Module.fDoGraphUseShape");
147	fNumPointsMin = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumPointsMin");
148	fNumPointsMax = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumPointsMax");
149	fNumHitsMinimum = (Int_t)Conf()->GetNum("TrackDirection2Module.fNumHitsMinimum");
150	fErrAngle = Conf()->GetNum("TrackDirection2Module.fErrAngle")/RadToDeg();
151
152	if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="linear") useLFplane = kTRUE;
153	else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="median") useMFplane = kTRUE;
154	else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="hough") useHFplane = kTRUE;
155	else if (Conf()->GetStr("TrackDirection2Module.fTDPFitMethod")=="rootrobust") useRFplane = kTRUE;
156	else Msg(EsafMsg::Panic) << "Wrong config fTDPFitMethod value in TrackDirection2Module.cfg" << MsgDispatch;
157
158	if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="linear") useLFaa1 = kTRUE;
159	else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="median") useMFaa1 = kTRUE;
160	else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="hough") useHFaa1 = kTRUE;
161	else if (Conf()->GetStr("TrackDirection2Module.fAA1FitMethod")=="rootrobust") useRFaa1 = kTRUE;
162	else Msg(EsafMsg::Panic) << "Wrong config fAA1FitMethod value in TrackDirection2Module.cfg" << MsgDispatch;
163
164	if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA1") fMethodIdentifier = kAA1;
165	else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE1") fMethodIdentifier = kNE1;
166	else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AA2") fMethodIdentifier = kAA2;
167	else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="NE2") fMethodIdentifier = kNE2;
168	else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="AE1") fMethodIdentifier = kAE1;
169	else if (Conf()->GetStr("TrackDirection2Module.fMethod")=="all") fMethodIdentifier = kALL;
170	else Msg(EsafMsg::Panic) << "Wrong config fMethod value in TrackDirection2Module.cfg" << MsgDispatch;
171
172	fDoHough = Conf()->GetBool("TrackDirection2Module.fUseHough");
173	fOptionSelectionHough = Conf()->GetBool("TrackDirection2Module.fUseHoughwithselection");
174	fDoShapeSelection = Conf()->GetBool("TrackDirection2Module.fUseShape");
175	if ( fDoShapeSelection) {
176	fMulti1 = Conf()->GetNum("TrackDirection2Module.fMulti1");
177	fMulti2 = Conf()->GetNum("TrackDirection2Module.fMulti2");
178	}
179
180	fUseShapeSelectioninModules = Conf()->GetBool("TrackDirection2Module.fUseShapeInModules");
181	fFixTmaxNumeric = Conf()->GetBool("TrackDirection2Module.fFixTmaxNumeric");
182	fDebugInfo = Conf()->GetBool("TrackDirection2Module.fDebugInfo");
183	fMinuitOutputLevel = (Int_t)Conf()->GetNum("TrackDirection2Module.fMinuitOutputLevel");
184	if ( fMinuitOutputLevel < -1 \|\| fMinuitOutputLevel > 3)
185	Msg(EsafMsg::Panic) << "Wrong fMinuitOutputLevel value in TrackDirection2Module.cfg" << MsgDispatch;
186
187	ConfigFileParser *pConfigGeneralEuso = Config::Get()->GetCF("General","Euso");
188	fHISS = pConfigGeneralEuso->GetNum("Euso.fAltitude");
189
190	ConfigFileParser* pConf = Config::Get()->GetCF("Electronics","MacroCell");
191	fGTU = pConf->GetNum("MacroCell.fGtuTimeLength");
192	fGTU/=1000.;
193
194
195
196	return kTRUE;
197	}
198
199	//_____________________________________________________________________________
200	Bool_t TrackDirection2Module::PreProcess() {
201	//
202	// Pre-process of the reco event
203	//
204
205	fData.Clear();
206	fAA1done = kFALSE;
207	fAA1nan = kFALSE;
208	fNumPoints = 0;
209	fNumHits = 0;
210	fNumPointsSel = 0;
211	fNumHitsSel = 0;
212	fQuality = -1;
213
214	fCentroid.SetXYZ(0,0,0);
215	fNorm.SetXYZ(0,0,0);
216	fW.SetXYZ(0,0,0);
217	fU.SetXYZ(0,0,0);
218	fTrueDir.SetXYZ(0,0,0);
219	fTrueNorm.SetXYZ(0,0,0);
220	fTrueMax.SetXYZ(0,0,0);
221	fEASDir.SetXYZ(0,0,0);
222
223	fAngularSpeed = 0;
224	fBeta = 0;
225	fBetaInit = 0;
226	fHmax = 0;
227	fRmax = 0;
228	fTrueTheta = 0;
229	fTruePhi = 0;
230	fTHETAloc = 0;
231	fTHETAreco = 0;
232	fPHIreco = 0;
233	fTmaxFit = 0;
234	fDeltaTheta = 0;
235	fDeltaPhi = 0;
236	fDeltaEASDir = 0;
237	fDeltaTDP = 0;
238
239	return kTRUE;
240	}
241
242	//_____________________________________________________________________________
243	Bool_t TrackDirection2Module::Process(RecoEvent *ev) {
244
245	fEv = ev;
246	RecoEventHeader header = fEv->GetHeader();
247	ERunParameters *runpars = header.GetRunPars();
248	// gtu length in microseconds
249	fGtuLength = fEv->GetHeader().GetGtuLength()/microsecond;
250
251	if(fDoGraphUseShape){
252	gDirectory->cd("/");
253	string pathname = Form("TD2ProcessRecoEvent%d", fEv->GetHeader().GetNum());
254	TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
255	path->cd();
256	}
257
258	RecoGlobalData* gdPattReco = fEv->FindGlobalData("PatternRecognition");
259	if (!gdPattReco) {
260	Msg(EsafMsg::Warning) << "Not found Pattern Recognition" << MsgDispatch;
261
262	Int_t totalNumPoints = fEv->GetHeader().GetNumActiveFee();
263	Msg(EsafMsg::Info) <<"TotalNumPoints = "<<totalNumPoints<<MsgDispatch;
264	if ( totalNumPoints <= 0 ) {
265	Msg(EsafMsg::Warning) << "No points in event! Exit from this event." << MsgDispatch;
266	return kFALSE;
267	} else if ( totalNumPoints < fNumPointsMin ) {
268	Msg(EsafMsg::Warning) << "Not enough points in event! Exit from this event." << MsgDispatch;
269	}
270	fPointsId.clear();
271	Bool_t processOnlySignal = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal");
272
273	// fill points with current threshold
274	FillPoints( processOnlySignal );
275
276	// check if there are enough points
277	if (fNumPoints < fNumPointsMin ) {
278	Msg(EsafMsg::Warning) << "Too few pixels selected. Trying to low the threshold." << MsgDispatch;
279	if ( !processOnlySignal ) {
280	for(Int_t i=Abs(fNumHitsMinimum)-1; i>0; i-- ) {
281	fNumHitsMinimum = i * Sign(1, fNumHitsMinimum);
282	FillPoints( processOnlySignal );
283	if ( fNumPoints < fNumPointsMin ) break;
284	}
285	} else {
286	Int_t counter(2);
287	Double_t snratio = fSignalToNoise;
288	while (fNumPoints < fNumPointsMin) {
289	fSignalToNoise = snratio/counter;
290	FillPoints(processOnlySignal);
291	counter++;
292	}
293	}
294	if (fNumPoints < fNumPointsMin ) {
295	Msg(EsafMsg::Warning) << "Too few pixels also with lower thresholds. Event terminated." << MsgDispatch;
296	return kFALSE;
297	}
298	}
299
300	// check if there are too many points
301	if (fNumPoints > fNumPointsMax ) {
302	Msg(EsafMsg::Warning) << "Too many pixels selected. Trying to raise the threshold." << MsgDispatch;
303	Int_t hits = Abs(fNumHitsMinimum) + 1;
304	Int_t counter(2);
305	Double_t snratio = fSignalToNoise;
306	while (fNumPoints > fNumPointsMax) {
307	fNumHitsMinimum = hits * Sign(1, fNumHitsMinimum);
308	fSignalToNoise = snratio * counter;
309	FillPoints( processOnlySignal );
310	hits++;
311	counter++;
312	}
313	if (fNumPoints < fNumPointsMin ) {
314	Msg(EsafMsg::Warning) << "Raising the threshold there are too few pixels. Event terminated." << MsgDispatch;
315	return kFALSE;
316	}
317	}
318
319	if ( processOnlySignal ) {
320	Msg(EsafMsg::Info) << "Num of pixels with S/N ratio >= " << fSignalToNoise << " = " << fNumPoints << MsgDispatch;
321	} else if ( fNumHitsMinimum >= 0 ) {
322	Msg(EsafMsg::Info) << "Num. of pixels with at least " << Abs(fNumHitsMinimum) << " hits = " << fNumPoints << MsgDispatch;
323	} else {
324	Msg(EsafMsg::Info) << "Num. of pixels over " << Abs(fNumHitsMinimum) << " sigma from background mean = " << fNumPoints << MsgDispatch;
325	}
326	} else {
327	if ( (gdPattReco->GetObj("SelectedPixels") == NULL) ) {
328	Msg(EsafMsg::Warning) << "No Pattern Recognition data in event" << MsgDispatch;
329	fNumPoints = 0; fQuality = -1;
330	return kFALSE;
331	}
332	fPointsId = (vector<Int_t>) gdPattReco->GetObj("SelectedPixels");
333	if (fPointsId.size() == 0) {
334	fNumPoints = 0;
335	Msg(EsafMsg::Warning) << "No Pattern Recognition data in event" << MsgDispatch;
336	fQuality = -1;
337	return kFALSE;
338	} else {
339	fNumPoints = fPointsId.size();
340	Msg(EsafMsg::Info) << "Pattern Recognition data found with " << fNumPoints << " points" << MsgDispatch;
341	}
342	}
343
344	if(fNumPoints < fNumPointsMin){
345	fQuality = -1;
346	Msg(EsafMsg::Warning) << "Not enough pixels selected ( minimum is "<< fNumPointsMin <<" ): terminate this event." << MsgDispatch;
347	return kFALSE;
348	} else if(fNumPoints > fNumPointsMax ) {
349	fQuality = -1;
350	Msg(EsafMsg::Warning) << "Too many pixels selected ( maximum is "<< fNumPointsMax <<" ): terminate this event." << MsgDispatch;
351	return kFALSE;
352	} else {
353	fQuality = 0;
354	}
355
356	fNumHits = 0;
357	for(Int_t i=0; i<fNumPoints; i++) {
358	RecoPixelData *pix = fEv->GetRecoPixelData( fPointsId[i] );
359	TVector3 pos(0,0,0);
360	pos.SetMagThetaPhi( 1,pix->GetTheta(), pix->GetPhi()+Pi() );
361	fData.fSpPos.push_back( pos );
362	Int_t pixid = pix->GetPixelId();
363	fData.fPixXYZ.push_back( runpars->PixelCenter(pixid) );
364	fData.fSigmaPhi.push_back( pix->GetPhiSigma() );
365	fData.fSigmaTheta.push_back( pix->GetThetaSigma() );
366	fData.fTime.push_back( 0.01fGtuLength(Double_t)pix->GetGtu() ); // microseconds/100
367	fData.fHits.push_back( pix->GetCounts() );
368	fNumHits += pix->GetCounts();
369	}
370	fData.fNumPoints = fNumPoints;
371	fData.fNumHits = fNumHits;
372	fData.fGtuLength = fGtuLength;
373
374	Msg(EsafMsg::Info) << "Processing " << fNumPoints << " pixel and " << fNumHits << " hits" << MsgDispatch;
375
376
377	fTrueTheta = fEv->GetHeader().GetTrueTheta();
378	fTruePhi = fEv->GetHeader().GetTruePhi();
379	fTrueDir.SetMagThetaPhi( 1., fTrueTheta, fTruePhi);
380	Msg(EsafMsg::Info) << "TRUTH: Theta = " << fTrueThetaRadToDeg() << " deg, Phi = " << fTruePhiRadToDeg() << " deg, Energy = " << header.GetTrueEnergy() << " MeV." << MsgDispatch;
381
382	TVector3 fTrueMaxOldsdr = fEv->GetHeader().GetTrueShowerMaxPos();
383	Double_t xmax = fTrueMaxOldsdr.x()/km;
384	Double_t ymax = fTrueMaxOldsdr.y()/km;
385	Double_t zmax = fHISS-(fTrueMaxOldsdr.z()/km);
386	fTrueMax.SetXYZ(xmax,ymax,zmax);
387	fTrueMax.SetPhi(fTrueMax.Phi()+Pi()); // reconstruiction refsys
388	fTrueNorm = fTrueMax.Cross(fTrueDir);
389	fTrueNorm.SetMag(1.);
390	if ( fTrueNorm.Z() > 0 ) fTrueNorm *= -1;
391	fTrueNorm.SetPhi(fTrueNorm.Phi()+Pi()); // reconstruction refsys
392
393	Bool_t fShapeSelectionDone = kFALSE;
394	if ( fDoHough ) {
395	UseHoughandFindPlane();
396	fShapeSelectionDone = kFALSE;
397	} else if ( !fDoShapeSelection ) {
398	FindPlane();
399	fShapeSelectionDone = kFALSE;
400	} else if ( fDoShapeSelection ) {
401	UseShapeandFindPlane();
402	if( fNumPointsSel < fNumPointsMin ){
403	Msg(EsafMsg::Info) << "Impossible to use selection by shape: only " << fNumPointsSel << " pixels selected!" <<MsgDispatch;
404	fShapeSelectionDone = kFALSE;
405	FindPlane();
406	} else { fShapeSelectionDone = kFALSE; }
407	}
408
409	//Double_t fangleDirNorm = fTrueDir.Angle(fNorm);
410	fDeltaTDP = fTrueNorm.Angle(fNorm);
411
412	// with fNorm these are the fundamental vectors of the TDP
413	fW = (TVector3(0,0,1).Cross(fNorm)).Unit();
414	fU = (fNorm.Cross(fW)).Unit();
415
416	Msg(EsafMsg::Info) << "Found TrackDetectorPlane with error = " << fDeltaTDP*RadToDeg() << " deg" << MsgDispatch;
417
418	/* Msg(EsafMsg::Info) << "True Norm X = " << fTrueNorm.X() << " Y = " << fTrueNorm.Y() << " Z = " << fTrueNorm.Z()
419	<< " Theta = " << fTrueNorm.Theta()RadToDeg() << " Phi = " << fTrueNorm.Phi()RadToDeg() << MsgDispatch;
420	Msg(EsafMsg::Info) << "Reco Norm X = " << fNorm.X() << " Y = " << fNorm.Y() << " Z = " << fNorm.Z()
421	<< " Theta = " << fNorm.Theta()RadToDeg() << " Phi = " << fNorm.Phi()RadToDeg() << MsgDispatch;
422	TVector3 fTrueW = (TVector3(0,0,1).Cross(fTrueNorm)).Unit();
423	Msg(EsafMsg::Info) << "True fW X = " << fTrueW.X() << " Y = " << fTrueW.Y() << " Z = " << fTrueW.Z()
424	<< " Theta = " << fTrueW.Theta()RadToDeg() << " Phi = " << fTrueW.Phi()RadToDeg() << MsgDispatch;
425	Msg(EsafMsg::Info) << "Reco fW X = " << fW.X() << " Y = " << fW.Y() << " Z = " << fW.Z()
426	<< " Theta = " << fW.Theta()RadToDeg() << " Phi = " << fW.Phi()RadToDeg() << MsgDispatch; */
427
428
429	/*if ( !fData.fSpPosSel.size() ) {
430	fQuality = -1;
431	Msg(EsafMsg::Warning) << "No enough selected points. Skip this event. " << MsgDispatch;
432	return kFALSE;
433	}*/
434
435	if( fDoHough \|\| (fUseShapeSelectioninModules && fShapeSelectionDone ) ){
436	fNumPoints = fNumPointsSel;
437	fNumHits = fNumHitsSel;
438	fData.fNumPoints = fNumPoints;
439	fData.fNumHits = fNumHits;
440	fData.fSpPos.clear();
441	fData.fTime.clear();
442	fData.fHits.clear();
443	fData.fSigmaTheta.clear();
444	fData.fSigmaPhi.clear();
445	for(Int_t i=0;i<fNumPointsSel;i++){
446	fData.fSpPos.push_back(fData.fSpPosSel[i]);
447	fData.fTime.push_back(fData.fTimeSel[i]);
448	fData.fHits.push_back(fData.fHitsSel[i]);
449	fData.fSigmaTheta.push_back(fData.fSigmaThetaSel[i]);
450	fData.fSigmaPhi.push_back(fData.fSigmaPhiSel[i]);
451	}
452	}
453
454	switch (fMethodIdentifier) {
455	case kAA1: AA1(); break;
456	case kNE1: NE1(); break;
457	case kAA2: AA2(); break;
458	case kNE2: NE2(); break;
459	case kAE1: AE1(); break;
460	case kALL: All(); break;
461	default: Msg(EsafMsg::Panic) << "Wrong fMethod config value in TrackDirection2Module.cfg" << MsgDispatch; break;
462	}
463
464	fData.Clear();
465
466	//save data
467	MyData()->Add("NormX",fNorm.X());
468	MyData()->Add("NormY",fNorm.Y());
469	MyData()->Add("NormZ",fNorm.Z());
470	MyData()->Add("WAxisX",fW.X());
471	MyData()->Add("WAxisY",fW.Y());
472	MyData()->Add("WAxisZ",fW.Z());
473	MyData()->Add("Theta",fTHETAreco);
474	MyData()->Add("Phi",fPHIreco);
475
476	return kTRUE;
477	}
478
479	//_____________________________________________________________________________
480	void TrackDirection2Module::FillPoints(Bool_t processOnlySignal) {
481	//
482	// Search the points over the threshold in case no pattern recognition
483	// is previously done
484	//
485
486	fPointsId.clear();
487	Int_t totalNumPoints = fEv->GetHeader().GetNumActiveFee();
488	ERunParameters *runpars = fEv->GetHeader().GetRunPars();
489	fNumPoints = 0;
490	for( Int_t i=0; i<totalNumPoints; i++ ) {
491	RecoPixelData pix = (RecoPixelData)fEv->GetRecoPixelData(i);
492	if ( !processOnlySignal ) {
493	Int_t threshold = 0;
494	if ( fNumHitsMinimum >= 0 ) threshold = Abs(fNumHitsMinimum);
495	else {
496	Int_t uid = fEv->GetRecoPixelData(i)->GetPixelId();
497	Double_t rate = runpars->GetNightGlowRateByUId(uid)*fGtuLength;
498	threshold = (Int_t)Ceil(rate + Abs(fNumHitsMinimum)*Sqrt(rate));
499	}
500	if ( pix->GetCounts() >= threshold ) {
501	fPointsId.push_back(i);
502	fNumPoints++;
503	}
504	} else {
505	Int_t signal_counts = pix->GetCountsSignal();
506	Int_t noise_counts = pix->GetCountsNoise();
507	Double_t ratio (-1);
508	if ( noise_counts > 0 ) ratio = 1.e0*signal_counts/Sqrt(noise_counts);
509	else if ( signal_counts > 0 ) ratio = 1.e10;
510	if (ratio >= fSignalToNoise ) {
511	fPointsId.push_back(i);
512	fNumPoints++;
513	}
514	}
515	}
516	}
517
518	//_____________________________________________________________________________
519	Bool_t TrackDirection2Module::PostProcess() {
520	//
521	// Post-processing method
522	//
523
524	if ( fAA1nan ) fQuality = -1;
525	if( fQuality == 0) {
526	fRecoEventsCounter++;
527	vecDeltaTDP.push_back(fDeltaTDP*RadToDeg());
528	if( fMethodIdentifier == kALL ) {
529	vecDeltaEASDirAA1.push_back(fDeltaEASDirAA1*RadToDeg());
530	vecDeltaEASDirAA2.push_back(fDeltaEASDirAA2*RadToDeg());
531	vecDeltaEASDirAE1.push_back(fDeltaEASDirAE1*RadToDeg());
532	vecDeltaEASDirNE1.push_back(fDeltaEASDirNE1*RadToDeg());
533	vecDeltaEASDirNE2.push_back(fDeltaEASDirNE2*RadToDeg());
534	} else {
535	vecDeltaEASDir.push_back(fDeltaEASDir*RadToDeg());
536	}
537	}
538
539	RecoGlobalData* data = fEv->FindCreateGlobalData("TrackDirection");
540	data->CleanMaps();
541	data->SetIssuer(GetName());
542	data->Add("NormX",fNorm.X());
543	data->Add("NormY",fNorm.Y());
544	data->Add("NormZ",fNorm.Z());
545	data->Add("WAxisX",fW.X());
546	data->Add("WAxisY",fW.Y());
547	data->Add("WAxisZ",fW.Z());
548	data->Add("Theta",fTHETAreco);
549	data->Add("Phi",fPHIreco);
550
551	TmpMemoryClean();
552
553	return kTRUE;
554
555	}
556
557	//_____________________________________________________________________________
558	Bool_t TrackDirection2Module::Done() {
559	//
560	// Module done method. Do some statistics about the reconstructed events
561	//
562
563	Msg(EsafMsg::Info) << "Number of reconstructed events = " << fRecoEventsCounter << MsgDispatch;
564	if (fRecoEventsCounter < 2) return kTRUE;
565
566	/*gDirectory->cd("/");
567	string pathname = "TD2Resolution";
568	TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
569	path->cd();
570
571	string aa1 = "AA1";
572	string aa2 = "AA2";
573	string ae1 = "AE1";
574	string ne1 = "NE1";
575	string ne2 = "NE2";
576	string tdp = "TDP";
577
578	DoStat(vecDeltaTDP,tdp);
579	if ( fMethodIdentifier == kALL ){
580	DoStat(vecDeltaEASDirAA1,aa1);
581	DoStat(vecDeltaEASDirAA2,aa2);
582	DoStat(vecDeltaEASDirAE1,ae1);
583	DoStat(vecDeltaEASDirNE1,ne1);
584	DoStat(vecDeltaEASDirNE2,ne2);
585	} else {
586	if ( fMethodIdentifier == kAA1 ) DoStat(vecDeltaEASDir,aa1);
587	if ( fMethodIdentifier == kNE1 ) DoStat(vecDeltaEASDir,ne1);
588	if ( fMethodIdentifier == kAA2 ) DoStat(vecDeltaEASDir,aa2);
589	if ( fMethodIdentifier == kNE2 ) DoStat(vecDeltaEASDir,ne2);
590	if ( fMethodIdentifier == kAE1 ) DoStat(vecDeltaEASDir,ae1);
591	}*/
592
593	vecDeltaEASDir.clear();
594	vecDeltaEASDirAA1.clear();
595	vecDeltaEASDirAA2.clear();
596	vecDeltaEASDirAE1.clear();
597	vecDeltaEASDirNE1.clear();
598	vecDeltaEASDirNE2.clear();
599	vecDeltaTDP.clear();
600	fRecoEventsCounter = 0;
601
602	Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
603	return kTRUE;
604
605	}
606
607	//_____________________________________________________________________________
608	void TrackDirection2Module::DoStat(vector<Double_t> err, string namemethod){
609	//
610	// Statistics of reconstructed events.
611	//
612
613	Int_t best(0);
614	Int_t med(0);
615	Int_t worst(0);
616	Int_t numbins=400;
617	Int_t maxDeltaEAS=20;
618	TH1F *hR=new TH1F("hR","hR",numbins,0,maxDeltaEAS);
619
620	for( Int_t i=0; i<fRecoEventsCounter; i++ ) {
621	hR->Fill(err[i]);
622	if ( err[i] < fStat1 ) best++;
623	if ( err[i] >= fStat1 && err[i] <= fStat2 ) med++;
624	if ( err[i] > fStat2 ) worst++;
625	}
626
627	Bool_t foundresolution = kFALSE;
628	Stat_t under68(0);
629	Float_t frac(0);
630	Double_t resolution(0);
631	for(Int_t i=0;i<numbins;i++){
632	under68 += hR->GetBinContent(i);
633	frac=100*(Float_t)under68/((Float_t)fRecoEventsCounter);
634	if(frac>68){
635	resolution=i*maxDeltaEAS/((Float_t)numbins);
636	Msg(EsafMsg::Info)<<"Eas-Dir reconstruction ("<<namemethod.c_str()<<"):\tResolution(68%)="<<resolution<<" deg"<<MsgDispatch;
637	foundresolution = kTRUE;
638	break;
639	}
640	}
641	if (!foundresolution)
642	Msg(EsafMsg::Info)<<"Eas-Dir reconstruction (" << namemethod.c_str() <<"):\tResolution(68%) not found, too bad!!" << MsgDispatch;
643
644	string titleh = namemethod.c_str();
645	titleh += Form(" Resolution(68%%)=%g deg", resolution);
646	hR->SetTitle(titleh.c_str());
647	hR->GetXaxis()->SetTitle("#Psi (deg)");
648	hR->Write();
649	delete hR;
650
651	Msg(EsafMsg::Info) << "\tEvents with error under " << fStat1<< " deg = " <<best<< " ( " << ((Float_t)best/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
652	Msg(EsafMsg::Info) << "\tEvents with error between "<<fStat1<<" and "<<fStat2<<" deg = " <<med<< " ( " << ((Float_t)med/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
653	Msg(EsafMsg::Info) << "\tEvents with error above "<<fStat2<<" deg = " <<worst<< " ( " << ((Float_t)worst/fRecoEventsCounter*100.) << "% )" << MsgDispatch;
654
655	return;
656	}
657
658	//_____________________________________________________________________________
659	void TrackDirection2Module::UserMemoryClean() {
660	//
661	// User memory clean
662	//
663	}
664
665	//_____________________________________________________________________________
666	void TrackDirection2Module::TmpMemoryClean() {
667	//
668	// Tmp memory clean
669	//
670
671	vector<Int_t> dummy0; dummy0.swap(fPointsId);
672
673	// vector<Double_t> dummy1; dummy1.swap(vecDeltaEASDir);
674	// vector<Double_t> dummy2; dummy2.swap(vecDeltaEASDirAA1);
675	// vector<Double_t> dummy3; dummy3.swap(vecDeltaEASDirAA2);
676	// vector<Double_t> dummy4; dummy4.swap(vecDeltaEASDirNE1);
677	// vector<Double_t> dummy5; dummy5.swap(vecDeltaEASDirNE2);
678	// vector<Double_t> dummy6; dummy6.swap(vecDeltaEASDirAE1);
679	// vector<Double_t> dummy7; dummy7.swap(vecDeltaTDP);
680
681	fData.ClearMemory();
682
683	}
684
685	//_____________________________________________________________________________
686	Bool_t TrackDirection2Module::SaveRootData(RecoRootEvent *fRecoRootEvent) {
687	//
688	// Save data in the reco rootfile
689	//
690
691	fRecoRootEvent->GetRecoTrackDirection2().SetErrorTDP(fDeltaTDP*RadToDeg());
692	fRecoRootEvent->GetRecoTrackDirection2().SetQuality(fQuality); // -1 if reconstruction is not possible in this context, otherwise 0
693
694	if (fMethodIdentifier == kALL) {
695	fRecoRootEvent->GetRecoTrackDirection2().SetAA2(fTHETArecoAA2RadToDeg(), fPHIrecoAA2RadToDeg(), fDeltaEASDirAA2*RadToDeg());
696	fRecoRootEvent->GetRecoTrackDirection2().SetNE1(fTHETArecoNE1RadToDeg(), fPHIrecoNE1RadToDeg(), fDeltaEASDirNE1*RadToDeg());
697	fRecoRootEvent->GetRecoTrackDirection2().SetNE2(fTHETArecoNE2RadToDeg(), fPHIrecoNE2RadToDeg(), fDeltaEASDirNE2*RadToDeg());
698	fRecoRootEvent->GetRecoTrackDirection2().SetAE1(fTHETArecoAE1RadToDeg(), fPHIrecoAE1RadToDeg(), fDeltaEASDirAE1*RadToDeg());
699	}
700	fRecoRootEvent->GetRecoTrackDirection2().SetAA1(fTHETArecoAA1RadToDeg(), fPHIrecoAA1RadToDeg(), fDeltaEASDirAA1*RadToDeg());
701	fRecoRootEvent->GetRecoTrackDirection2().SetTheta(fTHETAreco*RadToDeg());
702	fRecoRootEvent->GetRecoTrackDirection2().SetErrorTheta(fDeltaTheta*RadToDeg());
703	fRecoRootEvent->GetRecoTrackDirection2().SetPhi(fPHIreco*RadToDeg());
704	fRecoRootEvent->GetRecoTrackDirection2().SetErrorPhi(fDeltaPhi*RadToDeg());
705	fRecoRootEvent->GetRecoTrackDirection2().SetErrorDirection(fDeltaEASDir*RadToDeg());
706
707	return kTRUE;
708	}
709
710	//______________________________________________________________________________
711	void TrackDirection2Module::UseHoughandFindPlane(){
712	//
713	// Find track-detector plane (TDP) using Hough Transform
714	//
715
716	Msg(EsafMsg::Info) << "UseHoughandFindPlane()" << MsgDispatch;
717
718	vector<Double_t> x;
719	vector<Double_t> y;
720	vector<Double_t> t;
721	vector<Int_t> c;
722	for( Int_t i=0; i<fNumPoints; i++ ) {
723	TVector3 dummy = fData.fSpPos[i];
724	t.push_back(fData.fTime[i]);
725	x.push_back(dummy.x()/dummy.z());
726	y.push_back(dummy.y()/dummy.z());
727	c.push_back(fData.fHits[i]);
728	}
729
730	// errors for the Hough fit
731	Float_t ex = 0.001;
732	Float_t ey = 0.001;
733	Float_t et = fGTU*0.01;
734
735	// use Hough fit to select track points
736	// x-t Hough fit
737	HoughFit *xtFit = new HoughFit( t, x, c, fOptionSelectionHough,et,ex,fPointsId);
738	Double_t vX = xtFit->GetSlope();
739	Double_t wX = xtFit->GetOffset();
740	vector<Int_t> idSelX = xtFit->GetIdSel();
741	delete xtFit;
742
743	// y-t Hough fit
744	HoughFit *ytFit = new HoughFit( t, y, c, fOptionSelectionHough,et,ey,fPointsId);
745	Double_t vY = ytFit->GetSlope();
746	Double_t wY = ytFit->GetOffset();
747	vector<Int_t> idSelY = ytFit->GetIdSel();
748	delete ytFit;
749
750	// compute the normal vector to the TDP
751	fNorm.SetXYZ(vY,-vX,wYvX-wXvY);
752	fNorm.SetMag(1.);
753	if ( fNorm.Z() > 0 ) fNorm *= -1;
754
755	//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
756	if( fDoGraphUseShape ){
757	Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
758	Double_t xline[fNumPoints],yline[fNumPoints];
759	for( Int_t i=0; i<fNumPoints; i++ ) {
760	TVector3 dummy = fData.fSpPos[i];
761	tprog[i] = fData.fTime[i];
762	xline[i] = vX*fData.fTime[i] + wX;
763	yline[i] = vY*fData.fTime[i]+wY;
764	xprog[i] = dummy.x()/dummy.z();
765	yprog[i] = dummy.y()/dummy.z();
766	}
767	TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
768	gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
769	TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
770	gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
771	TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
772	gypro->SetNameTitle("gypro","yprog vs time");gypro->SetMarkerSize(.5);gypro->SetMarkerStyle(23);
773	TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
774	gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
775	TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
776	gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
777	TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
778	cUSFP->Divide(3,1);
779	cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");
780	cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");
781	cUSFP->cd(3);gxypro->Draw("AP");
782	cUSFP->Write();
783	}
784
785	//cout<<"fPointsId.size()="<<fPointsId.size()<<"\tidSelX.size()="<<idSelX.size()<<"\tidSelY.size()="<<idSelY.size()<<endl;
786	vector<Int_t> idSelectedMerged(idSelY.size()+idSelX.size());
787	merge(idSelX.begin(), idSelX.end(),idSelY.begin(), idSelY.end(),idSelectedMerged.begin());
788	vector<Int_t>::iterator different = unique(idSelectedMerged.begin(),idSelectedMerged.end());
789	idSelectedMerged.erase(different, idSelectedMerged.end());
790
791	fNumPointsSel = idSelectedMerged.size();
792	Msg(EsafMsg::Info) << "Selected " << fNumPointsSel << " points over " << fPointsId.size() << MsgDispatch;
793
794	fNumHitsSel = 0;
795	vector<Int_t> csel;
796	vector<Double_t> xsel;
797	vector<Double_t> ysel;
798	vector<Double_t> errxsel;
799	vector<Double_t> errysel;
800	vector<Double_t> tsel;
801	vector<Double_t> xsel1;
802	vector<Double_t> ysel1;
803	vector<Double_t> tsel1;
804
805	for(Int_t i=0; i<fNumPointsSel; i++) {
806	RecoPixelData *pix = fEv->GetRecoPixelData( idSelectedMerged[i] );
807	TVector3 pos(0,0,0);
808	Double_t phitmp = pix->GetPhi() + Pi();
809	pos.SetMagThetaPhi( 1, pix->GetTheta(), phitmp );
810	xsel.push_back( pos.x()/pos.z() );
811	ysel.push_back( pos.y()/pos.z() );
812	tsel.push_back( 0.01fGtuLengthpix->GetGtu() );
813	csel.push_back( pix->GetCounts() );
814	fNumHitsSel += csel[i];
815	Double_t errphi = pix->GetPhiSigma();
816	Double_t errtheta = pix->GetThetaSigma();
817	fData.fSpPosSel.push_back(pos);
818	fData.fTimeSel.push_back(tsel[i]);
819	fData.fHitsSel.push_back(csel[i]);
820	fData.fSigmaThetaSel.push_back(errtheta);
821	fData.fSigmaPhiSel.push_back(errphi);
822	Double_t dx = DeltaX(pos.Theta(), pos.Phi(), errtheta, errphi);
823	Double_t dy = DeltaY(pos.Theta(), pos.Phi(), errtheta, errphi);
824	errxsel.push_back(dx/csel[i]);
825	errysel.push_back(dy/csel[i]);
826	// add hits for median fit
827	for (Int_t j(0); j<pix->GetCounts(); j++) {
828	xsel1.push_back( pos.x()/pos.z() );
829	ysel1.push_back( pos.y()/pos.z() );
830	tsel1.push_back( 0.01fGtuLengthpix->GetGtu() );
831	}
832	}
833
834	// if there are > 10 points selected recompute the TDP more precisely
835	if( fNumPointsSel > 10 ){
836
837
838	TVector3 normHF, normLF, normMF, normRF;
839	Double_t vX2(0), vY2(0), wX2(0), wY2(0);
840	if ( fDebugInfo \|\| useHFplane ) { // recalculate normal to TDP using Hough fit
841	HoughFit *xhf1 = new HoughFit( tsel, xsel, csel, 0,et,ex);
842	vX2 = xhf1->GetSlope();
843	wX2 = xhf1->GetOffset();
844	delete xhf1;
845	HoughFit *yhf1 = new HoughFit( tsel, ysel, csel, 0,et,ey);
846	vY2 = yhf1->GetSlope();
847	wY2 = yhf1->GetOffset();
848	delete yhf1;
849
850	normHF.SetXYZ(vY2,-vX2,wY2vX2-wX2vY2);
851	normHF.SetMag(1.);
852	if ( normHF.Z() > 0 ) normHF *= -1;
853	Double_t fDeltaTDPhh = fTrueNorm.Angle(normHF);
854	Msg(EsafMsg::Info) << "TDP error (hough selection + hough fit) = " << fDeltaTDPhh*RadToDeg() << MsgDispatch;
855	}
856
857	if ( fDebugInfo \|\| useLFplane ) { // recalculate normal to TDP using least squares fit
858	LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPointsSel, tsel, xsel, errxsel );
859	Double_t vX2ls = xlf1->GetSlope();
860	Double_t wX2ls = xlf1->GetOffset();
861	delete xlf1;
862	LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPointsSel, tsel, ysel, errysel );
863	Double_t vY2ls = ylf1->GetSlope();
864	Double_t wY2ls = ylf1->GetOffset();
865	delete ylf1;
866
867	normLF.SetXYZ(vY2ls,-vX2ls,wY2lsvX2ls-wX2lsvY2ls);
868	normLF.SetMag(1.); if ( normLF.Z() > 0 ) normLF *= -1;
869	Double_t fDeltaTDPhls = fTrueNorm.Angle(normLF);
870	Msg(EsafMsg::Info) << "TDP error (hough selection + linear fit) = " << fDeltaTDPhls*RadToDeg() << MsgDispatch;
871	}
872
873	if ( fDebugInfo \|\| useHFplane ) { // recalculate normal to TDP using median fit
874	MedianFit *xmf1 = new MedianFit( fNumHitsSel, tsel1, xsel1 ); // CHANGED
875	Double_t vX2m = xmf1->GetSlope();
876	Double_t wX2m = xmf1->GetOffset();
877	delete xmf1;
878	MedianFit *ymf1 = new MedianFit( fNumHitsSel, tsel1, ysel1 ); // CHANGED
879	Double_t vY2m = ymf1->GetSlope();
880	Double_t wY2m = ymf1->GetOffset();
881	delete ymf1;
882
883	normMF.SetXYZ(vY2m,-vX2m,wY2mvX2m-wX2mvY2m);
884	normMF.SetMag(1.);
885	if ( normMF.Z() > 0 ) normMF *= -1;
886	Double_t fDeltaTDPhm = fTrueNorm.Angle(normMF);
887	Msg(EsafMsg::Info) << "TDP error (hough selection + median fit) = " << fDeltaTDPhm*RadToDeg() << MsgDispatch;
888	}
889
890	if ( fDebugInfo \|\| useRFplane ) { // recalculate the TDP using the root LTS fit
891	Double_t *t = new Double_t[1];
892	TLinearFitter *xrf1 = new TLinearFitter(1,"1 ++ x");
893	TLinearFitter *yrf1 = new TLinearFitter(1,"1 ++ x");
894	for( Int_t i(0); i<fNumPointsSel; i++ ) {
895	t[0] = tsel[i];
896	xrf1->AddPoint(t,xsel[i],errxsel[i]);
897	yrf1->AddPoint(t,ysel[i],errysel[i]);
898	}
899	xrf1->EvalRobust();
900	Double_t vX2r = xrf1->GetParameter(1);
901	Double_t wX2r = xrf1->GetParameter(0);
902	delete xrf1;
903	yrf1->EvalRobust();
904	Double_t vY2r = yrf1->GetParameter(1);
905	Double_t wY2r = yrf1->GetParameter(0);
906	delete yrf1;
907
908	normRF.SetXYZ(vY2r,-vX2r,wY2rvX2r-wX2rvY2r);
909	normRF.SetMag(1.); if ( normRF.Z() > 0 ) normRF *= -1;
910	Double_t fDeltaTDPhr = fTrueNorm.Angle(normRF);
911	Msg(EsafMsg::Info) << "TDP error (hough selection + root lts fit) = " << fDeltaTDPhr*RadToDeg() << MsgDispatch;
912	}
913
914	if ( useLFplane ) fNorm = normLF;
915	else if ( useMFplane ) fNorm = normMF;
916	else if ( useHFplane ) fNorm = normHF;
917	else if ( useRFplane ) fNorm = normRF;
918	fNorm.SetMag(1.);
919	if ( fNorm.Z() > 0 ) fNorm *= -1;
920
921	//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
922	if( fDoGraphUseShape ){
923	Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel],xline2[fNumPointsSel],yline2[fNumPointsSel];
924	for( Int_t i=0; i<fNumPointsSel; i++ ) {
925	tprog2[i] = tsel[i];
926	xline2[i] = vX2*tsel[i]+wX2;
927	yline2[i] = vY2*tsel[i]+wY2;
928	xprog2[i] = xsel[i];
929	yprog2[i] = ysel[i];
930	}
931	TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
932	gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");gxypro2->SetMarkerSize(.5);gxypro2->SetMarkerStyle(23);
933	TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
934	gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");gxpro2->SetMarkerSize(.5);gxpro2->SetMarkerStyle(23);
935	TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
936	gypro2->SetNameTitle("gypro2","yprog2 vs time2");gypro2->SetMarkerSize(.5);gypro2->SetMarkerStyle(23);
937	TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
938	gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
939	TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
940	gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
941	TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
942	cUSFP2->Divide(3,1);
943	cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");
944	cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");
945	cUSFP2->cd(3);gxypro2->Draw("AP");
946	cUSFP2->Write();
947	}
948	}
949	}
950
951	//______________________________________________________________________________
952	void TrackDirection2Module::FindPlane(){
953	//
954	// Find the track-detector plane (TDP)
955	//
956
957	Msg(EsafMsg::Info) << "FindPlane()" << MsgDispatch;
958
959	vector<Double_t> xpro;
960	vector<Double_t> ypro;
961	vector<Double_t> errxpro;
962	vector<Double_t> errypro;
963	vector<Double_t> tpro;
964	vector<Int_t> cpro;
965	vector<Double_t> xpro1;
966	vector<Double_t> ypro1;
967	vector<Double_t> tpro1;
968
969	for( Int_t i=0; i<fNumPoints; i++ ) {
970	TVector3 dummy = fData.fSpPos[i];
971	Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
972	Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
973	xpro.push_back(dummy.X()/dummy.Z());
974	ypro.push_back(dummy.Y()/dummy.Z());
975	cpro.push_back(fData.fHits[i]);
976	tpro.push_back(fData.fTime[i]);
977	errxpro.push_back(dx/fData.fHits[i]);
978	errypro.push_back(dy/fData.fHits[i]);
979	// add hits for median fit
980	for( Int_t j(0); j<fData.fHits[i]; j++ ) {
981	xpro1.push_back(dummy.X()/dummy.Z());
982	ypro1.push_back(dummy.Y()/dummy.Z());
983	tpro1.push_back(fData.fTime[i]);
984	}
985	}
986
987	TVector3 normHF, normLF, normMF, normRF;
988	if ( fDebugInfo \|\| useHFplane ) { // calculate normal to TDP using Hough fit
989	// errors for the Hough fit
990	Float_t ex = 0.001;
991	Float_t ey = 0.001;
992	Float_t et = fGtuLength*0.01;
993	HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
994	Double_t vX2 = xhf1->GetSlope();
995	Double_t wX2 = xhf1->GetOffset();
996	delete xhf1;
997	HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
998	Double_t vY2 = yhf1->GetSlope();
999	Double_t wY2 = yhf1->GetOffset();
1000	delete yhf1;
1001
1002	normHF.SetXYZ(vY2,-vX2,wY2vX2-wX2vY2);
1003	normHF.SetMag(1.);
1004	if ( normHF.Z() > 0 ) normHF *= -1;
1005	Double_t fDeltaTDPhh = fTrueNorm.Angle(normHF);
1006	Msg(EsafMsg::Info) << "TDP error (hough fit) = " << fDeltaTDPhh*RadToDeg() << MsgDispatch;
1007	}
1008
1009	if ( fDebugInfo \|\| useLFplane ) { // calculate normal to TDP using least squares fit
1010	LeastSquaresFit *xlf1 = new LeastSquaresFit( fNumPoints, tpro, xpro, errxpro );
1011	Double_t vX2ls = xlf1->GetSlope();
1012	Double_t wX2ls = xlf1->GetOffset();
1013	delete xlf1;
1014	LeastSquaresFit *ylf1 = new LeastSquaresFit( fNumPoints, tpro, ypro, errypro );
1015	Double_t vY2ls = ylf1->GetSlope();
1016	Double_t wY2ls = ylf1->GetOffset();
1017	delete ylf1;
1018
1019	normLF.SetXYZ(vY2ls,-vX2ls,wY2lsvX2ls-wX2lsvY2ls);
1020	normLF.SetMag(1.); if ( normLF.Z() > 0 ) normLF *= -1;
1021	Double_t fDeltaTDPhls = fTrueNorm.Angle(normLF);
1022	Msg(EsafMsg::Info) << "TDP error (linear fit) = " << fDeltaTDPhls*RadToDeg() << MsgDispatch;
1023	}
1024
1025	if ( fDebugInfo \|\| useHFplane ) { // recalculate normal to TDP using median fit
1026	MedianFit *xmf1 = new MedianFit( fNumHits, tpro1, xpro1 );
1027	Double_t vX2m = xmf1->GetSlope();
1028	Double_t wX2m = xmf1->GetOffset();
1029	delete xmf1;
1030	MedianFit *ymf1 = new MedianFit( fNumHits, tpro1, ypro1 );
1031	Double_t vY2m = ymf1->GetSlope();
1032	Double_t wY2m = ymf1->GetOffset();
1033	delete ymf1;
1034
1035	normMF.SetXYZ(vY2m,-vX2m,wY2mvX2m-wX2mvY2m);
1036	normMF.SetMag(1.);
1037	if ( normMF.Z() > 0 ) normMF *= -1;
1038	Double_t fDeltaTDPhm = fTrueNorm.Angle(normMF);
1039	Msg(EsafMsg::Info) << "TDP error (median fit) = " << fDeltaTDPhm*RadToDeg() << MsgDispatch;
1040	}
1041
1042	if ( fDebugInfo \|\| useRFplane ) { // recalculate the TDP using the root LTS fit
1043	Double_t *t = new Double_t[1];
1044	TLinearFitter *xrf1 = new TLinearFitter(1,"1 ++ x");
1045	TLinearFitter *yrf1 = new TLinearFitter(1,"1 ++ x");
1046	for( UInt_t i(0); i<tpro.size(); i++ ) {
1047	t[0] = tpro[i];
1048	xrf1->AddPoint(t,xpro[i],errxpro[i]);
1049	yrf1->AddPoint(t,ypro[i],errypro[i]);
1050	}
1051	xrf1->EvalRobust();
1052	Double_t vX2r = xrf1->GetParameter(1);
1053	Double_t wX2r = xrf1->GetParameter(0);
1054	delete xrf1;
1055	yrf1->EvalRobust();
1056	Double_t vY2r = yrf1->GetParameter(1);
1057	Double_t wY2r = yrf1->GetParameter(0);
1058	delete yrf1;
1059
1060	normRF.SetXYZ(vY2r,-vX2r,wY2rvX2r-wX2rvY2r);
1061	normRF.SetMag(1.); if ( normRF.Z() > 0 ) normRF *= -1;
1062	Double_t fDeltaTDPhr = fTrueNorm.Angle(normRF);
1063	Msg(EsafMsg::Info) << "TDP error (root lts fit) = " << fDeltaTDPhr*RadToDeg() << MsgDispatch;
1064	}
1065
1066	if ( useLFplane ) fNorm = normLF;
1067	else if ( useMFplane ) fNorm = normMF;
1068	else if ( useHFplane ) fNorm = normHF;
1069	else if ( useRFplane ) fNorm = normRF;
1070	fNorm.SetMag(1.);
1071	if ( fNorm.Z() > 0 ) fNorm *= -1;
1072
1073	cpro.clear();
1074	xpro.clear();
1075	ypro.clear();
1076	tpro.clear();
1077	errxpro.clear();
1078	errypro.clear();
1079	}
1080
1081	//______________________________________________________________________________
1082	void TrackDirection2Module::UseShapeandFindPlane(){
1083	//
1084	// Find the track-detector plane (TDP) using the shape method
1085	//
1086
1087	vector<Int_t> cpro;
1088	vector<Double_t> xpro;
1089	vector<Double_t> ypro;
1090	vector<Double_t> errxpro;
1091	vector<Double_t> errypro;
1092	vector<Double_t> tpro;
1093	vector<Double_t> xpro1;
1094	vector<Double_t> ypro1;
1095	vector<Double_t> tpro1;
1096
1097	for( Int_t i=0; i<fNumPoints; i++ ) {
1098	TVector3 dummy = fData.fSpPos[i];
1099	Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1100	Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1101	xpro.push_back(dummy.X()/dummy.Z());
1102	ypro.push_back(dummy.Y()/dummy.Z());
1103	tpro.push_back(fData.fTime[i]);
1104	cpro.push_back(fData.fHits[i]);
1105	errxpro.push_back(dx/fData.fHits[i]);
1106	errypro.push_back(dy/fData.fHits[i]);
1107	// add hits for median fit
1108	for(Int_t j(0); j<fData.fHits[i]; i++) {
1109	xpro1.push_back(dummy.X()/dummy.Z());
1110	ypro1.push_back(dummy.Y()/dummy.Z());
1111	tpro1.push_back(fData.fTime[i]);
1112	}
1113	}
1114
1115	// FIRST STEP: median fit of points
1116	MedianFit *xMf1 = new MedianFit(fNumHits, tpro1, xpro1); //CHANGED
1117	Double_t sX = xMf1->GetSlope();
1118	Double_t oX = xMf1->GetOffset();
1119	Double_t aX = xMf1->GetAbsoluteDeviation();
1120	delete xMf1;
1121	MedianFit *yMf1 = new MedianFit(fNumHits, tpro1, ypro1); //CHANGED
1122	Double_t sY = yMf1->GetSlope();
1123	Double_t oY = yMf1->GetOffset();
1124	Double_t aY = yMf1->GetAbsoluteDeviation();
1125	delete yMf1;
1126
1127	//use the following only to debug (to see how the projection on plane Z=1 versus time appear)
1128	if( fDoGraphUseShape ){
1129	Double_t xprog[fNumPoints], yprog[fNumPoints], tprog[fNumPoints];
1130	Double_t xline[fNumPoints], yline[fNumPoints];
1131	Double_t xlineup[fNumPoints], ylineup[fNumPoints];
1132	Double_t xlinedown[fNumPoints], ylinedown[fNumPoints];
1133	for( Int_t i=0; i<fNumPoints; i++ ) {
1134	TVector3 dummy = fData.fSpPos[i];
1135	tprog[i] = fData.fTime[i];
1136	xline[i] = sX * fData.fTime[i] + oX;
1137	xlineup[i] = xline[i] + fMulti1 * aX;
1138	xlinedown[i] = xline[i] - fMulti1 * aX;
1139	yline[i] = sY * fData.fTime[i] + oY;
1140	ylineup[i] = yline[i] + fMulti1 * aY;
1141	ylinedown[i] = yline[i] - fMulti1 * aY;
1142	xprog[i] = dummy.X()/dummy.Z();
1143	yprog[i] = dummy.Y()/dummy.Z();
1144	}
1145	TGraph *gxypro = new TGraph(fNumPoints,yprog,xprog);
1146	gxypro->SetNameTitle("gxypro","xprog vs yprog");gxypro->SetMarkerSize(.5);gxypro->SetMarkerStyle(23);
1147	TGraph *gxpro = new TGraph(fNumPoints,tprog,xprog);
1148	gxpro->SetNameTitle("gxpro","xprog vs time");gxpro->SetMarkerSize(.5);gxpro->SetMarkerStyle(23);
1149	TGraph *gypro = new TGraph(fNumPoints,tprog,yprog);
1150	gypro->SetNameTitle("gypro","yprog vs time");
1151	gypro->SetMarkerSize(.5);
1152	gypro->SetMarkerStyle(23);
1153	TGraph *gxproretta = new TGraph(fNumPoints,tprog,xline);
1154	gxproretta->SetMarkerColor(4);gxproretta->SetLineColor(4);
1155	TGraph *gyproretta = new TGraph(fNumPoints,tprog,yline);
1156	gyproretta->SetMarkerColor(4);gyproretta->SetLineColor(4);
1157	TGraph *gxprorettaup = new TGraph(fNumPoints,tprog,xlineup);
1158	gxprorettaup->SetMarkerColor(2);gxprorettaup->SetLineColor(2);
1159	TGraph *gyprorettaup = new TGraph(fNumPoints,tprog,ylineup);
1160	gyprorettaup->SetMarkerColor(2);gyprorettaup->SetLineColor(2);
1161	TGraph *gxprorettadown = new TGraph(fNumPoints,tprog,xlinedown);
1162	gxprorettadown->SetMarkerColor(2);gxprorettadown->SetLineColor(2);
1163	TGraph *gyprorettadown = new TGraph(fNumPoints,tprog,ylinedown);
1164	gyprorettadown->SetMarkerColor(2);gyprorettadown->SetLineColor(2);
1165	TCanvas *cUSFP=new TCanvas("cUSFP","cUSFP",1);
1166	cUSFP->Divide(3,1);
1167	cUSFP->cd(1);gxpro->Draw("AP");gxproretta->Draw("PLsame");gxprorettaup->Draw("PLsame");gxprorettadown->Draw("PLsame");
1168	cUSFP->cd(2);gypro->Draw("AP");gyproretta->Draw("PLsame");gyprorettaup->Draw("PLsame");gyprorettadown->Draw("PLsame");
1169	cUSFP->cd(3);gxypro->Draw("AP");
1170	cUSFP->Write();
1171	}
1172
1173	// make a shape selection of track points
1174	vector<Double_t> xprosel;
1175	vector<Double_t> yprosel;
1176	vector<Double_t> errxprosel;
1177	vector<Double_t> erryprosel;
1178	vector<Double_t> tprosel;
1179	vector<Double_t> xprosel1;
1180	vector<Double_t> yprosel1;
1181	vector<Double_t> tprosel1;
1182
1183	vector<TVector3> fSpPosSeltmp;
1184	vector<Double_t> fTimeSeltmp;
1185	vector<Int_t> fHitsSeltmp;
1186	vector<Double_t> fSigmaThetaSeltmp;
1187	vector<Double_t> fSigmaPhiSeltmp;
1188
1189	Int_t fNumPointsSeltmp = 0;
1190	Int_t fNumHitsSeltmp = 0;
1191	for( Int_t i=0; i<fNumPoints; i++ ) {
1192	TVector3 dummy = fData.fSpPos[i];
1193	Double_t xprotmp = dummy.X()/dummy.Z();
1194	Double_t yprotmp = dummy.Y()/dummy.Z();
1195	Double_t tprotmp = fData.fTime[i];
1196	// shape selection with multiplier fMulti1
1197	if( Abs(xprotmp-sXtprotmp-oX) < fMulti1aX \|\| Abs(yprotmp-sYtprotmp-oY) < fMulti1aY) {
1198	Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1199	Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1200	fSpPosSeltmp.push_back(dummy);
1201	fTimeSeltmp.push_back(fData.fTime[i]);
1202	fHitsSeltmp.push_back(fData.fHits[i]);
1203	fSigmaThetaSeltmp.push_back(fData.fSigmaTheta[i]);
1204	fSigmaPhiSeltmp.push_back(fData.fSigmaPhi[i]);
1205	xprosel.push_back(xprotmp);
1206	yprosel.push_back(yprotmp);
1207	tprosel.push_back(tprotmp);
1208	errxprosel.push_back(dx/fData.fHits[i]);
1209	erryprosel.push_back(dy/fData.fHits[i]);
1210	fNumPointsSeltmp++;
1211	fNumHitsSeltmp += fData.fHits[i];
1212	//add hits for median fit
1213	for( Int_t j(0); j<fData.fHits[i]; j++) {
1214	xprosel1.push_back(xprotmp);
1215	yprosel1.push_back(yprotmp);
1216	tprosel1.push_back(tprotmp);
1217	}
1218	}
1219	}
1220	Msg(EsafMsg::Info) << "UseShapeandFindPlane(), first step: selected " << fNumPointsSeltmp << " from " << fNumPoints << " pixels" << MsgDispatch;
1221
1222	if ( fNumPointsSeltmp < fNumPointsMin ){ // return if not enough points
1223	fNumPointsSel = fNumPointsSeltmp;
1224	return;
1225	}
1226
1227	// SECOND STEP: new median fit of selected points
1228	MedianFit *xMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, xprosel1); //CHANGED
1229	Double_t sX2 = xMf2->GetSlope();
1230	Double_t oX2 = xMf2->GetOffset();
1231	Double_t aX2 = xMf2->GetAbsoluteDeviation();
1232	delete xMf2;
1233	MedianFit *yMf2 = new MedianFit(fNumHitsSeltmp, tprosel1, yprosel1); //CHANGED
1234	Double_t sY2 = yMf2->GetSlope();
1235	Double_t oY2 = yMf2->GetOffset();
1236	Double_t aY2 = yMf2->GetAbsoluteDeviation();
1237	delete yMf2;
1238
1239	// new shape selection of points from median fit results
1240	vector<Double_t> xprosel2;
1241	vector<Double_t> yprosel2;
1242	vector<Double_t> errxprosel2;
1243	vector<Double_t> erryprosel2;
1244	vector<Double_t> tprosel2;
1245	vector<Int_t> cprosel2;
1246	vector<Double_t> xprosel3;
1247	vector<Double_t> yprosel3;
1248	vector<Double_t> tprosel3;
1249
1250	fNumPointsSel = 0;
1251	fNumHitsSel = 0;
1252
1253	for( Int_t i=0; i<fNumPointsSeltmp; i++ ) {
1254	TVector3 dummy = fSpPosSeltmp[i];
1255	Double_t xprotmp = dummy.X()/dummy.Z();
1256	Double_t yprotmp = dummy.Y()/dummy.Z();
1257	Double_t tprotmp = fTimeSeltmp[i];
1258	// shape selection with multiplier fMulti2
1259	if( Abs(xprotmp-sX2tprotmp-oX2) < fMulti2aX2 \|\| TMath::Abs(yprotmp-sY2tprotmp-oY2) < fMulti2aY2 ){
1260	Double_t dx = DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1261	Double_t dy = DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1262	fData.fSpPosSel.push_back(dummy);
1263	fData.fTimeSel.push_back(fTimeSeltmp[i]);
1264	fData.fHitsSel.push_back(fHitsSeltmp[i]);
1265	fData.fSigmaThetaSel.push_back(fSigmaThetaSeltmp[i]);
1266	fData.fSigmaPhiSel.push_back(fSigmaPhiSeltmp[i]);
1267	cprosel2.push_back(fHitsSeltmp[i]);
1268	xprosel2.push_back(xprotmp);
1269	yprosel2.push_back(yprotmp);
1270	tprosel2.push_back(tprotmp);
1271	errxprosel2.push_back(dx/fHitsSeltmp[i]);
1272	erryprosel2.push_back(dy/fHitsSeltmp[i]);
1273	fNumPointsSel++;
1274	fNumHitsSel += fHitsSeltmp[i];
1275	//add hits for median fit
1276	for (Int_t j(0); j<fHitsSeltmp[i]; j++) {
1277	xprosel3.push_back(xprotmp);
1278	yprosel3.push_back(yprotmp);
1279	tprosel3.push_back(tprotmp);
1280	}
1281	}
1282	}
1283	Msg(EsafMsg::Info) << "UseShapeandFindPlane(), second step: selected " << fNumPointsSel << " from " << fNumPointsSeltmp << " pixels" << MsgDispatch;
1284
1285	if ( fNumPointsSel < fNumPointsMin ) return; // return if not enough points
1286
1287	Double_t vX(0), wX(0), vY(0), wY(0);
1288	// fit of selected points
1289	if ( useLFplane ) {
1290	LeastSquaresFit *xtFit = new LeastSquaresFit(fNumPointsSel, tprosel2, xprosel2, errxprosel2);
1291	vX = xtFit->GetSlope();
1292	wX = xtFit->GetOffset();
1293	delete xtFit;
1294	LeastSquaresFit *ytFit = new LeastSquaresFit(fNumPointsSel, tprosel2, yprosel2, erryprosel2);
1295	vY = ytFit->GetSlope();
1296	wY = ytFit->GetOffset();
1297	delete ytFit;
1298	}
1299
1300	if ( useMFplane ) {
1301	MedianFit *xtFit = new MedianFit(fNumHitsSel, tprosel3, xprosel3);
1302	vX = xtFit->GetSlope();
1303	wX = xtFit->GetOffset();
1304	delete xtFit;
1305	MedianFit *ytFit = new MedianFit(fNumHitsSel, tprosel3, yprosel3);
1306	vY = ytFit->GetSlope();
1307	wY = ytFit->GetOffset();
1308	delete ytFit;
1309	}
1310
1311	if ( useHFplane ) {
1312	Float_t ex = 0.001;
1313	Float_t ey = 0.001;
1314	Float_t et = fGtuLength*0.01;
1315	HoughFit *xhf1 = new HoughFit( tpro, xpro, cpro, 0,et,ex);
1316	vX = xhf1->GetSlope();
1317	wX = xhf1->GetOffset();
1318	delete xhf1;
1319	HoughFit *yhf1 = new HoughFit( tpro, ypro, cpro, 0,et,ey);
1320	vY = yhf1->GetSlope();
1321	wY = yhf1->GetOffset();
1322	delete yhf1;
1323	}
1324
1325	if ( useRFplane ) { // recalculate the TDP using the root LTS fit
1326	Double_t *t = new Double_t[1];
1327	TLinearFitter *xrf1 = new TLinearFitter(1,"1 ++ x");
1328	TLinearFitter *yrf1 = new TLinearFitter(1,"1 ++ x");
1329	for( Int_t i(0); i<fNumPointsSel; i++ ) {
1330	t[0] = tprosel2[i];
1331	xrf1->AddPoint(t,xprosel2[i],errxprosel2[i]);
1332	yrf1->AddPoint(t,yprosel2[i],erryprosel2[i]);
1333	}
1334	xrf1->EvalRobust();
1335	vX = xrf1->GetParameter(1);
1336	wX = xrf1->GetParameter(0);
1337	delete xrf1;
1338	yrf1->EvalRobust();
1339	vY = yrf1->GetParameter(1);
1340	wY = yrf1->GetParameter(0);
1341	delete yrf1;
1342	}
1343
1344	// calculate normal versor to the TDP
1345	fNorm.SetXYZ(vY,-vX,wYvX-wXvY);
1346	fNorm.SetMag(1.);
1347	if ( fNorm.Z() > 0 ) fNorm *= -1;
1348
1349	//use the following only to debug (to see how the projection on plane Z=1 versus time appear after selection)
1350	if( fDoGraphUseShape ) {
1351	Double_t xprog2[fNumPointsSel],yprog2[fNumPointsSel],tprog2[fNumPointsSel];
1352	Double_t xline2[fNumPointsSel],yline2[fNumPointsSel];
1353	//Double_t xline2up[fNumPoints],yline2up[fNumPoints];
1354	//Double_t xline2down[fNumPoints],yline2down[fNumPoints];
1355	for( Int_t i=0; i<fNumPointsSel; i++ ) {
1356	TVector3 dummy = fData.fSpPosSel[i];
1357	//xline2[i] = vX * fData.fTimeSel[i] + wX;
1358	//yline2[i] = vY * fData.fTimeSel[i] + wY;
1359	xline2[i] = vX * tprosel2[i] + wX;
1360	yline2[i] = vY * tprosel2[i] + wY;
1361	//xline2up[i] = xline2[i] + fMulti2 * a2X;
1362	//xline2down[i] = xline2[i] - fMulti2 * a2X;
1363	//yline2up[i] = yline2[i] + fMulti2 * a2Y;
1364	//yline2down[i] = yline2[i] - fMulti2 * a2Y;
1365	xprog2[i] = xprosel2[i]; // xprog2[i] = dummy.X()/dummy.Z();
1366	yprog2[i] = yprosel2[i]; // yprog2[i] = dummy.Y()/dummy.Z();
1367	tprog2[i] = tprosel2[i]; // tprog2[i] = fData.fTimeSel[i];
1368	}
1369	TGraph *gxypro2 = new TGraph(fNumPointsSel,yprog2,xprog2);
1370	gxypro2->SetNameTitle("gxypro2","xprog2 vs yprog2");
1371	gxypro2->SetMarkerSize(.5);
1372	gxypro2->SetMarkerStyle(23);
1373	TGraph *gxpro2 = new TGraph(fNumPointsSel,tprog2,xprog2);
1374	gxpro2->SetNameTitle("gxpro2","xprog2 vs time2");
1375	gxpro2->SetMarkerSize(.5);
1376	gxpro2->SetMarkerStyle(23);
1377	TGraph *gypro2 = new TGraph(fNumPointsSel,tprog2,yprog2);
1378	gypro2->SetNameTitle("gypro2","yprog2 vs time2");
1379	gypro2->SetMarkerSize(.5);
1380	gypro2->SetMarkerStyle(23);
1381	TGraph *gxproretta2 = new TGraph(fNumPointsSel,tprog2,xline2);
1382	gxproretta2->SetMarkerColor(4);gxproretta2->SetLineColor(4);
1383	TGraph *gyproretta2 = new TGraph(fNumPointsSel,tprog2,yline2);
1384	gyproretta2->SetMarkerColor(4);gyproretta2->SetLineColor(4);
1385	/TGraph gxproretta2up = new TGraph(fNumPointsSel,tprog2,xline2up);
1386	gxproretta2up->SetMarkerColor(2);gxproretta2up->SetLineColor(2);
1387	TGraph *gyproretta2up = new TGraph(fNumPointsSel,tprog2,yline2up);
1388	gyproretta2up->SetMarkerColor(2);gyproretta2up->SetLineColor(2);
1389	TGraph *gxproretta2down = new TGraph(fNumPointsSel,tprog2,xline2down);
1390	gxproretta2down->SetMarkerColor(2);gxproretta2down->SetLineColor(2);
1391	TGraph *gyproretta2down = new TGraph(fNumPointsSel,tprog2,yline2down);
1392	gyproretta2down->SetMarkerColor(2);gyproretta2down->SetLineColor(2);*/
1393	TCanvas *cUSFP2=new TCanvas("cUSFP2","cUSFP2",1);
1394	cUSFP2->Divide(3,1);
1395	cUSFP2->cd(1);gxpro2->Draw("AP");gxproretta2->Draw("PLsame");//gxproretta2up->Draw("PLsame");gxproretta2down->Draw("PLsame");
1396	cUSFP2->cd(2);gypro2->Draw("AP");gyproretta2->Draw("PLsame");//gyproretta2up->Draw("PLsame");gyproretta2down->Draw("PLsame");
1397	cUSFP2->cd(3);gxypro2->Draw("AP");
1398	cUSFP2->Write();
1399	}
1400
1401	// clear of containers
1402	cpro.clear(); xpro.clear(); ypro.clear(); tpro.clear();
1403	errxpro.clear(); errypro.clear();
1404	xprosel.clear(); yprosel.clear(); tprosel.clear();
1405	errxprosel.clear(); erryprosel.clear();
1406	cprosel2.clear(); xprosel2.clear(); yprosel2.clear();
1407	errxprosel2.clear(); erryprosel2.clear(); tprosel2.clear();
1408
1409	Msg(EsafMsg::Info) << "UseShapeandFindPlane(): " << fNumPointsSel << " pixel and " << fNumHitsSel << " hits selected " << MsgDispatch;
1410
1411	return;
1412	}
1413
1414	//______________________________________________________________________________
1415	void TrackDirection2Module::AA1() {
1416	//
1417	// ANALITIC APPROXIMATED 1 method.
1418	// Constant angular velocity of the shower approximation.
1419	// This method initialize the parameters for other all fit methods.
1420	//
1421
1422	Msg(EsafMsg::Info) << "Using method AA1()" << MsgDispatch;
1423
1424	if ( fNumPoints < 10 ) { fAA1done = kTRUE; fAA1nan = kTRUE; return; }
1425
1426	fData.fMedDir.SetXYZ(0,0,0);
1427	fData.fMedTime = 0;
1428	fData.fMedAlpha = 0;
1429	vector<Double_t> erralpha;
1430	vector<Int_t> counts;
1431	//vectors of data for median fit
1432	vector<Double_t> alphahits;
1433	vector<Double_t> timehits;
1434	Int_t numhits = 0;
1435
1436	for( Int_t i=0; i<fNumPoints; i++ ) {
1437	TVector3 dummy = fData.fSpPos[i];
1438	Double_t alphatmp = ACos(dummy*fW);
1439	fData.fAlpha.push_back(alphatmp);
1440	//erralpha.push_back(fErrAngle/Sqrt((Double_t)fData.fHits[i]));//FIXME
1441	erralpha.push_back(fGtuLength*0.01);
1442	counts.push_back(fData.fHits[i]);
1443	numhits += fData.fHits[i];
1444	fData.fMedTime += fData.fTime[i]*fData.fHits[i];
1445	fData.fMedDir += fData.fSpPos[i]*fData.fHits[i];
1446	//add hits for median fit
1447	for( Int_t j(0); j<fData.fHits[i]; j++ ) {
1448	alphahits.push_back(alphatmp);
1449	timehits.push_back(fData.fTime[i]);
1450	}
1451	}
1452
1453	fData.fApparentTimeLenght = (Int_t)((fData.fTime[fNumPoints-1]-fData.fTime[0])/fGtuLength);
1454	fData.fMedTime = fData.fMedTime/fNumHits;
1455	fData.fMedDir = (fData.fMedDir*(1/(Double_t)fNumHits)).Unit();
1456	fData.fMedAlpha = ACos(fData.fMedDir*fW);
1457
1458	//fit the angular velocity
1459	// shower angular velocity fit: least squares, median, hough and root lts
1460	Double_t fAngularSpeedLF(0), fQAngularSpeedLF(0), fAngularSpeedtodrawLF(0);
1461	Double_t fAngularSpeedMF(0), fQAngularSpeedMF(0), fAngularSpeedtodrawMF(0);
1462	Double_t fAngularSpeedHF(0), fQAngularSpeedHF(0), fAngularSpeedtodrawHF(0);
1463	Double_t fAngularSpeedRF(0), fQAngularSpeedRF(0), fAngularSpeedtodrawRF(0);
1464
1465	if ( fDebugInfo \|\| fDoGraphUseShape \|\| useLFaa1 ) {
1466	LeastSquaresFit *ASl = new LeastSquaresFit(fNumPoints, fData.fTime, fData.fAlpha, erralpha);
1467	fAngularSpeedLF = ASl->GetSlope()*0.01; // uom fix with C()
1468	fQAngularSpeedLF = ASl->GetOffset();
1469	fAngularSpeedtodrawLF = ASl->GetSlope();
1470	delete ASl;
1471	}
1472
1473	if ( fDebugInfo \|\| fDoGraphUseShape \|\| useMFaa1 ) {
1474	MedianFit *ASm = new MedianFit(numhits, timehits, alphahits); //CHANGED
1475	fAngularSpeedMF = ASm->GetSlope()*0.01; // uom fix with C()
1476	fQAngularSpeedMF = ASm->GetOffset();
1477	fAngularSpeedtodrawMF = ASm->GetSlope();
1478	delete ASm;
1479	}
1480
1481	if ( fDebugInfo \|\| fDoGraphUseShape \|\| useHFaa1 ) {
1482	Float_t et = 0.01; // errors for the hough fit
1483	Float_t ea = fErrAngle;
1484	HoughFit *ASh = new HoughFit( fData.fTime, fData.fAlpha, counts, 0,et,ea);
1485	fAngularSpeedHF = ASh->GetSlope()*0.01; // uom fix with C()
1486	fQAngularSpeedHF = ASh->GetOffset();
1487	fAngularSpeedtodrawHF = ASh->GetSlope();
1488	delete ASh;
1489	}
1490
1491
1492	if ( fDebugInfo \|\| fDoGraphUseShape \|\| useRFaa1 ) { // recalculate the TDP using the root LTS fit
1493	Double_t *t = new Double_t[1];
1494	TLinearFitter *ASr = new TLinearFitter(1,"1 ++ x");
1495	for( Int_t i(0); i<fNumPoints; i++ ) {
1496	t[0] = fData.fTime[i];
1497	ASr->AddPoint(t,fData.fAlpha[i],erralpha[i]);
1498	}
1499	ASr->EvalRobust();
1500	fAngularSpeedRF = ASr->GetParameter(1)*0.01; // uom fix with C()
1501	fQAngularSpeedRF = ASr->GetParameter(0);
1502	fAngularSpeedtodrawRF = ASr->GetParameter(1);
1503
1504	delete ASr;
1505
1506	}
1507
1508
1509	if (fDebugInfo)
1510	Msg(EsafMsg::Info) << "AA1() Angular Speed : LF = " << fAngularSpeedLF << " MF = " << fAngularSpeedMF << " HF = " << fAngularSpeedHF << " RF = " << fAngularSpeedRF << MsgDispatch;
1511
1512	// use the following only to debug (plot angles alphai vs time)
1513	if( fDoGraphUseShape ){
1514	Double_t alphai[fNumPoints],timei[fNumPoints],alphailinelf[fNumPoints],alphailinemf[fNumPoints],
1515	alphailinehf[fNumPoints], alphailinerf[fNumPoints];
1516	for( Int_t i=0; i<fNumPoints; i++ ) {
1517	alphai[i] = fData.fAlpha[i]*RadToDeg();
1518	timei[i] = fData.fTime[i];
1519	alphailinelf[i] = (fData.fTime[i]fAngularSpeedtodrawLF+fQAngularSpeedLF)RadToDeg();
1520	alphailinemf[i] = (fData.fTime[i]fAngularSpeedtodrawMF+fQAngularSpeedMF)RadToDeg();
1521	alphailinehf[i] = (fData.fTime[i]fAngularSpeedtodrawHF+fQAngularSpeedHF)RadToDeg();
1522	alphailinerf[i] = (fData.fTime[i]fAngularSpeedtodrawRF+fQAngularSpeedRF)RadToDeg();
1523	}
1524	TGraph *gang = new TGraph(fNumPoints,timei,alphai);
1525	gang->SetNameTitle("gang","alphai vs time;time;alpha");
1526	gang->SetMarkerSize(.5);gang->SetMarkerStyle(23);
1527	TGraph *ganglinelf = new TGraph(fNumPoints,timei,alphailinelf);ganglinelf->SetLineColor(kRed);
1528	TGraph *ganglinemf = new TGraph(fNumPoints,timei,alphailinemf);ganglinemf->SetLineColor(kBlue);
1529	TGraph *ganglinehf = new TGraph(fNumPoints,timei,alphailinehf);ganglinehf->SetLineColor(kGreen);
1530	TGraph *ganglinerf = new TGraph(fNumPoints,timei,alphailinerf);ganglinerf->SetLineColor(28);
1531	TCanvas *cAA1=new TCanvas("cAA1","cAA1",1);cAA1->Divide(1,1);
1532	cAA1->cd(1);
1533	gang->Draw("AP");
1534	ganglinelf->Draw("PLsame");
1535	ganglinemf->Draw("PLsame");
1536	ganglinehf->Draw("PLsame");
1537	ganglinerf->Draw("PLsame");
1538	TLegend *legend = new TLegend(0.8,0.9,1,0.6);
1539	legend->AddEntry(ganglinelf,"least squares fit","l");
1540	legend->AddEntry(ganglinemf,"median fit","l");
1541	legend->AddEntry(ganglinehf,"hough fit","l");
1542	legend->AddEntry(ganglinerf,"root lts fit","l");
1543	legend->Draw();
1544	cAA1->Write();
1545	delete cAA1;
1546	delete gang;
1547	delete ganglinelf;
1548	delete ganglinemf;
1549	delete ganglinehf;
1550	delete ganglinerf;
1551	delete legend;
1552	}
1553
1554	Double_t angspeedtodraw(0);
1555	// selection of fit method
1556	if ( useLFaa1 ) {
1557	fAngularSpeed = fAngularSpeedLF;
1558	fQAngularSpeed = fQAngularSpeedLF;
1559	angspeedtodraw = fAngularSpeedtodrawLF;
1560	} else if ( useMFaa1 ) {
1561	fAngularSpeed = fAngularSpeedMF;
1562	fQAngularSpeed = fQAngularSpeedMF;
1563	angspeedtodraw = fAngularSpeedtodrawMF;
1564	} else if ( useHFaa1 ) {
1565	fAngularSpeed = fAngularSpeedHF;
1566	fQAngularSpeed = fQAngularSpeedHF;
1567	angspeedtodraw = fAngularSpeedtodrawHF;
1568	} else if ( useRFaa1 ) {
1569	fAngularSpeed = fAngularSpeedRF;
1570	fQAngularSpeed = fQAngularSpeedRF;
1571	angspeedtodraw = fAngularSpeedtodrawRF;
1572	} else {
1573	fAngularSpeed = fAngularSpeedLF;
1574	fQAngularSpeed = fQAngularSpeedLF;
1575	angspeedtodraw = fAngularSpeedtodrawLF;
1576	Msg(EsafMsg::Warning) << "AA1() fit method non specified in config file. Saving least squares fit results." << MsgDispatch;
1577	}
1578
1579	// check if the result of the fit is a 'not a number'
1580	if ( IsNaN(fAngularSpeed) ) {
1581	Msg(EsafMsg::Warning) << "NaN angular speed, skipping this event" << MsgDispatch;
1582	fAA1done = kTRUE;
1583	fAA1nan = kTRUE;
1584	return;
1585	}
1586
1587	// select points and redo the fit
1588	if( !(Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal")) ) {
1589	alphahits.clear(); timehits.clear(); erralpha.clear(); counts.clear();
1590	ContainerData newdata; newdata.Clear();
1591	Double_t *distance = new Double_t[(const Int_t)fNumPoints];
1592	for( Int_t i(0); i<fNumPoints; i++ ) {
1593	distance[i] = Abs(angspeedtodraw*fData.fTime[i] + fQAngularSpeed - fData.fAlpha[i]);
1594	}
1595	Double_t distrms = RMS(fNumPoints, distance);
1596	Int_t numselpoints = 0;
1597	Int_t numselhits = 0;
1598	for( Int_t i(0); i<fNumPoints; i++) {
1599	if( distance[i] < 0.5*distrms ) {
1600	numselpoints++;
1601	erralpha.push_back(fGtuLength*0.01);
1602	counts.push_back(fData.fHits[i]);
1603	newdata.fSpPos.push_back(fData.fSpPos[i]);
1604	//newdata.fPixXYZ.push_back(fData.fPixXYZ[i]);
1605	//newdata.fPos.push_back(fData.fPos[i]);
1606	//newdata.fSpErr.push_back(fData.fSpErr[i]);
1607	newdata.fTime.push_back(fData.fTime[i]);
1608	//newdata.fTimeSel.push_back(fData.fTimeSel[i]);
1609	newdata.fAlpha.push_back(fData.fAlpha[i]);
1610	newdata.fSigmaPhi.push_back(fData.fSigmaPhi[i]);
1611	newdata.fSigmaTheta.push_back(fData.fSigmaTheta[i]);
1612	newdata.fHits.push_back(fData.fHits[i]);
1613	newdata.fMedTime += fData.fTime[i]*fData.fHits[i];
1614	newdata.fMedDir += fData.fSpPos[i]*fData.fHits[i];
1615	numselhits += fData.fHits[i];
1616	for( Int_t j(0); j<fData.fHits[i]; j++ ) {
1617	alphahits.push_back(fData.fAlpha[i]);
1618	timehits.push_back(fData.fTime[i]);
1619	}
1620	}
1621	}
1622
1623	if (numselpoints>10) newdata.fApparentTimeLenght = (Int_t)((newdata.fTime[numselpoints-1]-newdata.fTime[0])/fGtuLength);
1624	newdata.fMedTime = (newdata.fMedTime)/numselhits;
1625	newdata.fMedDir = (newdata.fMedDir*(1/(Double_t)numselhits)).Unit();
1626	newdata.fMedAlpha = ACos(newdata.fMedDir*fW);
1627	newdata.fGtuLength = fData.fGtuLength;
1628	newdata.fNumPoints = numselpoints;
1629	newdata.fNumHits = numselhits;
1630
1631	if (numselpoints>10) {
1632	fData.Clear();
1633	fData = newdata;
1634	}
1635	newdata.ClearMemory();
1636
1637	//redo the fit
1638
1639	if (numselpoints>10 && ( fDebugInfo \|\| fDoGraphUseShape \|\| useLFaa1 ) ) {
1640	LeastSquaresFit *ASl = new LeastSquaresFit(numselpoints, fData.fTime, fData.fAlpha, erralpha);
1641	fAngularSpeedLF = ASl->GetSlope()*0.01; // uom fix with C()
1642	fQAngularSpeedLF = ASl->GetOffset();
1643	fAngularSpeedtodrawLF = ASl->GetSlope();
1644	delete ASl;
1645	}
1646
1647	if (numselpoints>10 && ( fDebugInfo \|\| fDoGraphUseShape \|\| useMFaa1 ) ) {
1648	MedianFit *ASm = new MedianFit(numselhits, timehits, alphahits); //CHANGED
1649	fAngularSpeedMF = ASm->GetSlope()*0.01; // uom fix with C()
1650	fQAngularSpeedMF = ASm->GetOffset();
1651	fAngularSpeedtodrawMF = ASm->GetSlope();
1652	delete ASm;
1653	}
1654
1655	if (numselpoints>10 && ( fDebugInfo \|\| fDoGraphUseShape \|\| useHFaa1 ) ) {
1656	Float_t et = 0.01; // errors for the hough fit
1657	Float_t ea = fErrAngle;
1658	HoughFit *ASh = new HoughFit( fData.fTime, fData.fAlpha, counts, 0,et,ea);
1659	fAngularSpeedHF = ASh->GetSlope()*0.01; // uom fix with C()
1660	fQAngularSpeedHF = ASh->GetOffset();
1661	fAngularSpeedtodrawHF = ASh->GetSlope();
1662	delete ASh;
1663	}
1664
1665
1666	if (numselpoints>10 && ( fDebugInfo \|\| fDoGraphUseShape \|\| useRFaa1 ) ) {
1667	Double_t *t = new Double_t[1];
1668	TLinearFitter *ASr = new TLinearFitter(1,"1 ++ x");
1669	for( Int_t i(0); i<numselpoints; i++ ) {
1670	t[0] = fData.fTime[i];
1671	ASr->AddPoint(t,fData.fAlpha[i],erralpha[i]);
1672	}
1673	ASr->EvalRobust();
1674	fAngularSpeedRF = ASr->GetParameter(1)*0.01; // uom fix with C()
1675	fQAngularSpeedRF = ASr->GetParameter(0);
1676	fAngularSpeedtodrawRF = ASr->GetParameter(1);
1677
1678	delete ASr;
1679
1680	}
1681
1682
1683	if (fDebugInfo) {
1684	Msg(EsafMsg::Info) << "Selected " << numselpoints << " over " << fNumPoints << MsgDispatch;
1685	Msg(EsafMsg::Info) << "AA1() angular speed after selection: LF = " << fAngularSpeedLF << " MF = " << fAngularSpeedMF << " HF = " << fAngularSpeedHF << " RF = " << fAngularSpeedRF << MsgDispatch;
1686	}
1687
1688	// use the following only to debug (plot angles alphai vs time)
1689	if( fDoGraphUseShape ){
1690	Double_t alphai[fData.fNumPoints],timei[fData.fNumPoints],alphailinelf[fData.fNumPoints],alphailinemf[fData.fNumPoints],
1691	alphailinehf[fData.fNumPoints], alphailinerf[fData.fNumPoints];
1692	for( Int_t i=0; i<fData.fNumPoints; i++ ) {
1693	alphai[i] = fData.fAlpha[i]*RadToDeg();
1694	timei[i] = fData.fTime[i];
1695	alphailinelf[i] = (fData.fTime[i]fAngularSpeedtodrawLF+fQAngularSpeedLF)RadToDeg();
1696	alphailinemf[i] = (fData.fTime[i]fAngularSpeedtodrawMF+fQAngularSpeedMF)RadToDeg();
1697	alphailinehf[i] = (fData.fTime[i]fAngularSpeedtodrawHF+fQAngularSpeedHF)RadToDeg();
1698	alphailinerf[i] = (fData.fTime[i]fAngularSpeedtodrawRF+fQAngularSpeedRF)RadToDeg();
1699	}
1700	TGraph *gang = new TGraph(fData.fNumPoints,timei,alphai);
1701	gang->SetNameTitle("gangsel","alphai vs time;time;alpha");
1702	gang->SetMarkerSize(.5);gang->SetMarkerStyle(23);
1703	TGraph *ganglinelf = new TGraph(fData.fNumPoints,timei,alphailinelf);ganglinelf->SetLineColor(kRed);
1704	TGraph *ganglinemf = new TGraph(fData.fNumPoints,timei,alphailinemf);ganglinemf->SetLineColor(kBlue);
1705	TGraph *ganglinehf = new TGraph(fData.fNumPoints,timei,alphailinehf);ganglinehf->SetLineColor(kGreen);
1706	TGraph *ganglinerf = new TGraph(fData.fNumPoints,timei,alphailinerf);ganglinerf->SetLineColor(28);
1707	TCanvas *cAA1=new TCanvas("cAA1sel","cAA1sel",1);cAA1->Divide(1,1);
1708	cAA1->cd(1);
1709	gang->Draw("AP");
1710	ganglinelf->Draw("PLsame");
1711	ganglinemf->Draw("PLsame");
1712	ganglinehf->Draw("PLsame");
1713	ganglinerf->Draw("PLsame");
1714	TLegend *legend = new TLegend(0.8,0.9,1,0.6);
1715	legend->AddEntry(ganglinelf,"least squares fit","l");
1716	legend->AddEntry(ganglinemf,"median fit","l");
1717	legend->AddEntry(ganglinehf,"hough fit","l");
1718	legend->AddEntry(ganglinerf,"root lts fit","l");
1719	legend->Draw();
1720	cAA1->Write();
1721	delete cAA1;
1722	delete gang;
1723	delete ganglinelf;
1724	delete ganglinemf;
1725	delete ganglinehf;
1726	delete ganglinerf;
1727	delete legend;
1728	}
1729
1730	// selection of fit method
1731	if ( useLFaa1 ) {
1732	fAngularSpeed = fAngularSpeedLF;
1733	fQAngularSpeed = fQAngularSpeedLF;
1734	} else if ( useMFaa1 ) {
1735	fAngularSpeed = fAngularSpeedMF;
1736	fQAngularSpeed = fQAngularSpeedMF;
1737	} else if ( useHFaa1 ) {
1738	fAngularSpeed = fAngularSpeedHF;
1739	fQAngularSpeed = fQAngularSpeedHF;
1740	} else if ( useRFaa1 ) {
1741	fAngularSpeed = fAngularSpeedRF;
1742	fQAngularSpeed = fQAngularSpeedRF;
1743	} else {
1744	fAngularSpeed = fAngularSpeedLF;
1745	Msg(EsafMsg::Warning) << "AA1() fit method non specified in config file. Saving least squares fit results." << MsgDispatch;
1746	}
1747	}
1748
1749	// calculate beta and the shower direction
1750	fHmax = 5; // HMax initialized to 5 kilometers
1751	fRmax = CalculateRmax(fHmax);
1752
1753
1754	if ( fDebugInfo ) { // if debug display results for all fith methods
1755	fBetaInit = 2ATan(Clight()/kmmicrosecond/(fAngularSpeedLF*fRmax)) - fData.fMedAlpha;
1756	fBeta = Pi() - fBetaInit;
1757	while(fBeta > 2Pi()) fBeta = fBeta - 2Pi();
1758	CalculatefromBetaEASdir();
1759	Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirLF = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1760
1761	fBetaInit = 2ATan(Clight()/kmmicrosecond/(fAngularSpeedMF*fRmax)) - fData.fMedAlpha;
1762	fBeta = Pi() - fBetaInit;
1763	while(fBeta> 2Pi()) fBeta = fBeta - 2Pi();
1764	CalculatefromBetaEASdir();
1765	Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirMF = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1766
1767	fBetaInit = 2ATan(Clight()/kmmicrosecond/(fAngularSpeedHF*fRmax)) - fData.fMedAlpha;
1768	fBeta = Pi() - fBetaInit;
1769	while(fBeta > 2Pi()) fBeta = fBeta - 2Pi();
1770	CalculatefromBetaEASdir();
1771	Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirHF = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1772
1773	fBetaInit = 2ATan(Clight()/kmmicrosecond/(fAngularSpeedRF*fRmax)) - fData.fMedAlpha;
1774	fBeta = Pi() - fBetaInit;
1775	while(fBeta > 2Pi()) fBeta = fBeta - 2Pi();
1776	CalculatefromBetaEASdir();
1777	Msg(EsafMsg::Info) << "AA1(): fDeltaEASDirRF = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1778	}
1779	// first value of beta
1780	fBetaInit = 2ATan(Clight()/kmmicrosecond/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
1781	fBeta = Pi() - fBetaInit;
1782	while(fBeta > 2Pi()) fBeta = fBeta - 2Pi();
1783	CalculatefromBetaEASdir();
1784	if ( IsNaN(fTHETAreco) \|\| IsNaN(fPHIreco) ) { fAA1done = kTRUE; fAA1nan = kTRUE; return; }
1785
1786	// recalculate Hmax
1787	fHmax = FindHmax( fTHETAreco );
1788	fRmax = CalculateRmax( fHmax );
1789	TVector3 nmax = fData.fMedDir;
1790	fData.fVectorRmax.SetMagThetaPhi(fRmax,nmax.Theta(),nmax.Phi());
1791
1792	// recalculate beta
1793	fBetaInit = 2atan(Clight()/kmmicrosecond/(fAngularSpeed*fRmax)) - fData.fMedAlpha;
1794	fBeta = Pi() - fBetaInit;
1795	while(fBeta > 2Pi()) fBeta = fBeta - 2Pi();
1796	CalculatefromBetaEASdir();
1797	if ( IsNaN(fTHETAreco) \|\| IsNaN(fPHIreco) ) { fAA1done = kTRUE; fAA1nan = kTRUE; return; }
1798
1799	fHmax = FindHmax(fTHETAreco);
1800	Msg(EsafMsg::Info) << "method AA1(): fDeltaEASDir = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1801
1802	// always save results of this method in rootfile
1803	fTHETArecoAA1 = fTHETAreco;
1804	fPHIrecoAA1 = fPHIreco;
1805	fDeltaEASDirAA1 = fDeltaEASDir;
1806
1807	fAA1done = kTRUE;
1808	}
1809
1810	//______________________________________________________________________________
1811	void TrackDirection2Module::AA2() {
1812	//
1813	// ANALYTIC APPROXIMATED 2 method.
1814	// Approximation: Shower velocity on a plane perpendicular to the detector
1815	// axis is constant
1816	//
1817
1818	Msg(EsafMsg::Info) << "Using method AA2()" << MsgDispatch;
1819	if ( !fAA1done ) AA1();
1820
1821	if ( fAA1nan ) return;
1822
1823	vector<Double_t> xprov,yprov,errxyprov;
1824	vector<Double_t> errxprov,erryprov;
1825	for(Int_t i=0;i<fData.fNumPoints;i++){
1826	TVector3 dummy = fData.fSpPos[i];
1827	xprov.push_back(dummy.x()*(fHISS-fHmax)/dummy.z());
1828	yprov.push_back(dummy.y()*(fHISS-fHmax)/dummy.z());
1829	Double_t dx = (fHISS-fHmax) * DeltaX(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1830	Double_t dy = (fHISS-fHmax) * DeltaY(dummy.Theta(), dummy.Phi(), fData.fSigmaTheta[i], fData.fSigmaPhi[i]);
1831	errxyprov.push_back(1./Sqrt((Double_t)fData.fHits[i]));
1832	errxprov.push_back( dx / fData.fHits[i]);
1833	erryprov.push_back( dy / fData.fHits[i]);
1834	}
1835
1836	Double_t speedx,speedy;
1837	LeastSquaresFit *fitvx = new LeastSquaresFit(fData.fNumPoints, fData.fTime, xprov, errxprov);
1838	speedx = fitvx->GetSlope()*0.01; // why?
1839	delete fitvx;
1840	LeastSquaresFit *fitvy = new LeastSquaresFit(fData.fNumPoints, fData.fTime, yprov, erryprov);
1841	speedy = fitvy->GetSlope()*0.01; // why?
1842	delete fitvy;
1843
1844	Double_t speed = Sqrt(speedxspeedx + speedyspeedy);
1845	Double_t gthetamax = fData.fMedAlpha;
1846	Double_t betacontrol = fBeta;
1847	Double_t gammaV = 2ATan((speed/Clight()km/microsecond) * Sin(gthetamax));
1848	Double_t betaV(0);
1849
1850	if(betacontrolRadToDeg()>(90-gthetamaxRadToDeg()) && betacontrol*RadToDeg()<90){
1851	betaV = -gammaV + gthetamax;
1852	} else if (betacontrolRadToDeg() > (90-gthetamaxRadToDeg()) && betacontrol*RadToDeg() > 90 ) {
1853	betaV = gthetamax + gammaV;
1854	}
1855
1856	fBeta = betaV;
1857	while(fBeta > 2Pi()) fBeta=fBeta-2Pi();
1858	CalculatefromBetaEASdir();
1859	fHmax = FindHmax(fTHETAreco);
1860	Msg(EsafMsg::Info) << "method AA2(): fDeltaEASDir = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1861	}
1862
1863	//______________________________________________________________________________
1864	void TrackDirection2Module::NE1() {
1865	//
1866	// NUMERICAL EXACT 1 method
1867	// Chi-square minimization of the difference between arrival times of photons
1868	// measured and teoretically computed
1869	//
1870
1871	Msg(EsafMsg::Info) << "using method NE1()" << MsgDispatch;
1872	if( !fAA1done ) AA1();
1873	if ( fAA1nan ) return;
1874
1875	TMinuit *minuitNE1 = new TMinuit(3);
1876	minuitNE1->SetObjectFit(&fData);
1877	minuitNE1->SetPrintLevel(fMinuitOutputLevel);
1878
1879	TString namebeta="beta", nametmax="tmax", namehmax="hmax";
1880	Int_t ierflg(0);
1881	Double_t fcnout,edm,errdef;
1882	Int_t nvpar,nparx,icstat;
1883	Double_t value, errv, vlow,vup;
1884	Int_t gmaxcall = 1000;
1885	Double_t deltabeta,deltatmax;
1886	Double_t betaStart = Pi() - fBeta;
1887	Double_t hmaxStart = fHmax;
1888	Double_t tmaxStart = fData.fMedTime*100;
1889	Double_t step[3] = {0.2,fGTU,0.3}; //about 10 deg (beta), 2500ns (tmax), 0.3 km (hmax)
1890	Double_t arglist[10];
1891
1892	minuitNE1->SetFCN(ChiSquareNE1);
1893	minuitNE1->mnparm(0, "beta",betaStart , step[0], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER beta.");
1894	minuitNE1->mnparm(1, "tmax",tmaxStart , step[1], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
1895	minuitNE1->mnparm(2, "hmax",hmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER hmax.");
1896
1897	if( fFixTmaxNumeric ) {
1898	arglist[0] = 2;
1899	minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
1900	}
1901	arglist[0] = 3;
1902	minuitNE1->mnexcm("FIX",arglist,1,ierflg); if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
1903	arglist[0] = gmaxcall;
1904	minuitNE1->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) Printf(" ......UNABLE to minimize with MIGRAD");
1905
1906	minuitNE1->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
1907	minuitNE1->mnpout(0,namebeta, value,errv,vlow,vup,ierflg);
1908	fBeta = Pi() - value;
1909	deltabeta = errv;
1910
1911	minuitNE1->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
1912	fTmaxFit = value;
1913	deltatmax = errv;
1914	CalculatefromBetaEASdir();
1915	fHmax=FindHmax(fTHETAreco);
1916	Msg(EsafMsg::Info) << "method NE1(): fDeltaEASDir = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
1917	}
1918
1919	//______________________________________________________________________________
1920	void ChiSquareNE1(Int_t &npar, Double_t gin, Double_t &chisqnorm, Double_t par, Int_t iflag) {
1921	//
1922	// Function FCN called by minuit for the algorithm NE1
1923	//
1924
1925	chisqnorm = 0;
1926	ContainerData fDataMin = (ContainerData) gMinuit->GetObjectFit();
1927
1928	Int_t numhits(0);
1929	Double_t chisq(0);
1930	Double_t gerrort = fDataMin->fGtuLength;
1931	Double_t Rm = (fDataMin->fVectorRmax).Mag();
1932	Double_t Am = fDataMin->fMedAlpha;
1933	for(Int_t i=0;i<fDataMin->fNumPoints;i++){
1934	Double_t TimeExpected = ( par[1] - Rm * ( ( Sin(Am-fDataMin->fAlpha[i]) + Sin(fDataMin->fAlpha[i]+par[0]) - Sin(Am+par[0]) )/ Sin(fDataMin->fAlpha[i]+par[0]) )/Clight()*km/microsecond);
1935	chisq += Power((fDataMin->fTime[i]100-TimeExpected)/gerrort,2)fDataMin->fHits[i];
1936	numhits += fDataMin->fHits[i];
1937	}
1938
1939	chisqnorm = (Double_t)chisq/numhits;
1940	}
1941
1942
1943	//______________________________________________________________________________
1944	void TrackDirection2Module::NE2() {
1945	//
1946	// NUMERICAL EXACT 2 method
1947	// Chi-square minimization of angle between the versors of pixels in FOV
1948	// and the corresponding vectors from points of the track and the detector
1949	//
1950
1951	Msg(EsafMsg::Info) << "using method NE2()" << MsgDispatch;
1952	if( !fAA1done ) AA1();
1953	if ( fAA1nan ) return;
1954
1955	TMinuit *minuitNE2 = new TMinuit(3);
1956	minuitNE2->SetObjectFit(&fData);
1957	minuitNE2->SetPrintLevel(fMinuitOutputLevel);
1958
1959	TString nametheta="theta",namephi="phi",nametmax="tmax";
1960	Int_t ierflg(0);
1961	Double_t fcnout,edm,errdef;
1962	Int_t nvpar,nparx,icstat;
1963	Double_t value, errv, vlow,vup;
1964	Int_t gmaxcall = 1000;
1965	Double_t deltatheta,deltaphi,deltatmax;
1966
1967	Double_t thetastart = fTHETAreco;
1968	Double_t phistart = fPHIreco;
1969	Double_t tmaxStart = fData.fMedTime*100; //microseconds
1970	Double_t step[3] = {0.2,0.2,fGTU};
1971	Double_t arglist[10];
1972
1973	minuitNE2->SetFCN(ChiSquareNE2);
1974	minuitNE2->mnparm(0, "theta",thetastart , step[0], 0,0,ierflg); if (ierflg) Printf("......UNABLE TO DEFINE PARAMETER theta.");
1975	minuitNE2->mnparm(1, "phi",phistart , step[1], 0,0,ierflg); if (ierflg) Printf(".....UNABLE TO DEFINE PARAMETER phi.");
1976	minuitNE2->mnparm(2, "tmax",tmaxStart , step[2], 0,0,ierflg); if (ierflg) Printf(" ........UNABLE TO DEFINE PARAMETER tmax.");
1977
1978	if( fFixTmaxNumeric ) {
1979	arglist[0] = 3;
1980	minuitNE2->mnexcm("FIX",arglist,1,ierflg);if (ierflg) Printf(" .....UNABLE to fix parameter tmax");
1981	}
1982	arglist[0] = gmaxcall;
1983	minuitNE2->mnexcm("MIGRAD",arglist ,1 , ierflg);if (ierflg) Printf(" .......UNABLE to minimize with migrad");
1984
1985	minuitNE2->mnstat(fcnout,edm,errdef,nvpar,nparx,icstat);
1986	minuitNE2->mnpout(0,nametheta, value,errv,vlow,vup,ierflg);
1987	fTHETAreco = value;
1988	deltatheta = errv;
1989
1990	minuitNE2->mnpout(1,namephi, value,errv,vlow,vup,ierflg);
1991	fPHIreco = value;
1992	deltaphi = errv;
1993
1994	minuitNE2->mnpout(2,nametmax, value,errv,vlow,vup,ierflg);
1995	fTmaxFit = value;
1996	deltatmax = errv;
1997
1998	fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco+Pi()); // main euso system
1999	fTHETAreco = fEASDir.Theta();
2000	fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*Pi()) );
2001	fDeltaTheta = fTHETAreco-fTrueTheta;
2002	fDeltaPhi = fPHIreco-fTruePhi;
2003	fDeltaEASDir = fEASDir.Angle(fTrueDir);
2004
2005	fHmax = FindHmax(fTHETAreco);
2006	Msg(EsafMsg::Info) << "method NE2(): fDeltaEASDir = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
2007	}
2008
2009	//______________________________________________________________________
2010	void ChiSquareNE2(Int_t &npar, Double_t gin, Double_t &chisqnorm, Double_t par, Int_t iflag) {
2011	//
2012	// Function FCN called by minuit for the algorithm NE2
2013	//
2014
2015	ContainerData fDataMin = (ContainerData) gMinuit->GetObjectFit();
2016
2017	Int_t numhits(0);
2018	Double_t PidotRi(0), modRi(0), Li(0), betaprimo(0);
2019	TVector3 dummyRmax = fDataMin->fVectorRmax;
2020	Double_t Xm = dummyRmax.x();
2021	Double_t Ym = dummyRmax.y();
2022	Double_t Zm = dummyRmax.z();
2023	Double_t Rm = dummyRmax.Mag();
2024	Double_t gerrorpix = 0.0018; //rad, about 0.1 deg
2025	Double_t sctmp = Sin(par[0])*Cos(par[1]-Pi());
2026	Double_t sstmp = Sin(par[0])*Sin(par[1]-Pi());
2027	Double_t ctmp = -Cos(par[0]);
2028	Double_t thetaprimotmp = ACos((-Xmsctmp-Ymsstmp-Zm*ctmp)/Rm);
2029
2030	chisqnorm=0;
2031	for(Int_t i=0; i<fDataMin->fNumPoints; i++){
2032	betaprimo = 2ATan(1/((1/Tan(thetaprimotmp/2))+Clight()/kmmicrosecond(fDataMin->fTime[i]100-par[2])/(Rm*Sin(thetaprimotmp))));
2033	Li = Rm*(Cos(thetaprimotmp) - Sin(thetaprimotmp) / Tan(betaprimo));
2034	PidotRi = fDataMin->fSpPos[i].x()(Xm+Lisctmp) + fDataMin->fSpPos[i].y()(Ym+Lisstmp) + fDataMin->fSpPos[i].z()(Zm+Lictmp);
2035	modRi = Power(Xm+Lisctmp,2) + Power(Ym+Lisstmp,2) + Power(Zm+Li*ctmp,2);
2036	chisqnorm += Power((ACos(PidotRi/Sqrt(modRi)))/gerrorpix,2)*fDataMin->fHits[i];
2037	numhits += fDataMin->fHits[i];
2038	}
2039	chisqnorm = (Double_t)chisqnorm/numhits;
2040	}
2041
2042	//______________________________________________________________________________
2043	void TrackDirection2Module::AE1() {
2044	//
2045	// ANALYTICAL EXACT 1 method
2046	// Fit using exact relations between pixel directions in FOV and photons
2047	// arrival times. This method doesn't require the knowledge of the TDP.
2048	//
2049
2050	Msg(EsafMsg::Info) << "Using method AE1()" << MsgDispatch;
2051	if( !fAA1done ) AA1();
2052	if ( fAA1nan ) return;
2053
2054	vector<Double_t> alpha;
2055	vector<Double_t> Ri;
2056	vector<Double_t> dLi;
2057	vector<Double_t> timestart;
2058	vector<Double_t> nh;
2059	vector<Double_t> x;
2060	vector<Double_t> y;
2061	vector<Double_t> z;
2062	vector<Double_t> ex;
2063	vector<Double_t> ey;
2064	vector<Double_t> ez;
2065	Double_t gnthetamax = fData.fMedDir.Theta();
2066	Double_t gnphimax = fData.fMedDir.Phi();
2067	Double_t timemax = fData.fMedTime*100; // microseconds
2068
2069	Int_t pointsAE1(0),hitsAE1(0);
2070	Double_t Ritmp(0), dLitmp(0);
2071	for ( Int_t i=0 ; i<fData.fNumPoints; i++ ){
2072	TVector3 dummy = fData.fSpPos[i];
2073	Double_t dummyTime = fData.fTime[i]*100.;//mus
2074	Double_t cost = Cos(dummy.Theta())Cos(gnthetamax) + Sin(dummy.Theta())Sin(gnthetamax)*Cos(dummy.Phi()-gnphimax);
2075	if( Abs(dummyTime-timemax)/fGtuLength>1 && Abs(dummyTime-timemax)/fGtuLength>fData.fApparentTimeLenght/10 && fData.fHits[i]>2 ){
2076	if( dummyTime > timemax){
2077	Ritmp = (2Clight()/kmmicrosecond(dummyTime-timemax)fRmax + Power(Clight()/kmmicrosecond(dummyTime-timemax),2))/(2(fRmax(1-cost)
2078	+ Clight()/kmmicrosecond(dummyTime-timemax)));
2079	dLitmp = -Sqrt( RitmpRitmp + fRmaxfRmax - 2RitmpfRmax*cost );
2080	}
2081	if(dummyTime < timemax){
2082	Ritmp = (2Clight()/kmmicrosecond(-dummyTime+timemax)fRmax - Power(Clight()/kmmicrosecond(dummyTime-timemax),2))/(2(fRmax(-1+cost)
2083	+ Clight()/kmmicrosecond(-dummyTime+timemax)));
2084	dLitmp = Sqrt( RitmpRitmp + fRmaxfRmax - 2RitmpfRmax*cost );
2085	}
2086	alpha.push_back(fData.fAlpha[i]);
2087	Ri.push_back(Ritmp);
2088	dLi.push_back(dLitmp);
2089	timestart.push_back(-dLitmp/Clight()/km*microsecond);
2090	nh.push_back(fData.fHits[i]);
2091	x.push_back(RitmpSin(dummy.Theta())Cos(dummy.Phi()));
2092	ex.push_back(0.1/fData.fHits[i]);
2093	y.push_back(RitmpSin(dummy.Theta())Sin(dummy.Phi()));
2094	ey.push_back(0.1/fData.fHits[i]);
2095	z.push_back(Ritmp*Cos(dummy.Theta()));
2096	ez.push_back(0.1/fData.fHits[i]);
2097	hitsAE1 += fData.fHits[i];
2098	pointsAE1++;
2099	}
2100	}
2101	/*alpha.push_back(fData.fMedAlpha);
2102	dLi.push_back(0);
2103	timestart.push_back(0);
2104	Ri.push_back(fRmax);
2105	x.push_back(fData.fVectorRmax.x());
2106	y.push_back(fData.fVectorRmax.y());
2107	z.push_back(fData.fVectorRmax.z());*/
2108	Double_t vx,vy,vz;
2109	Double_t qx,qy,qz;
2110
2111	Double_t *t = new Double_t[1];
2112	TLinearFitter *ASx = new TLinearFitter(1,"1 ++ x");
2113	TLinearFitter *ASy = new TLinearFitter(1,"1 ++ x");
2114	TLinearFitter *ASz = new TLinearFitter(1,"1 ++ x");
2115	for( Int_t i(0); i<pointsAE1; i++ ) {
2116	t[0] = timestart[i];
2117	ASx->AddPoint(t,x[i],ex[i]);
2118	ASy->AddPoint(t,y[i],ey[i]);
2119	ASz->AddPoint(t,z[i],ez[i]);
2120	}
2121	(pointsAE1>10 ? ASx->EvalRobust() : ASx->Eval()); vx = ASx->GetParameter(1); qx = ASx->GetParameter(0);
2122	(pointsAE1>10 ? ASy->EvalRobust() : ASy->Eval()); vy = ASy->GetParameter(1); qy = ASy->GetParameter(0);
2123	(pointsAE1>10 ? ASz->EvalRobust() : ASz->Eval()); vz = ASz->GetParameter(1); qz = ASz->GetParameter(0);
2124	delete ASx; delete ASy; delete ASz;
2125
2126	/LeastSquaresFit fitvx = new LeastSquaresFit(pointsAE1, timestart, x, ex);
2127	vx = fitvx->GetSlope();
2128	delete fitvx;
2129	LeastSquaresFit *fitvz = new LeastSquaresFit(pointsAE1, timestart, z, ez);
2130	vz = fitvz->GetSlope();
2131	delete fitvz;
2132	LeastSquaresFit *fitvy = new LeastSquaresFit(pointsAE1, timestart, y, ey);
2133	vy = fitvy->GetSlope();
2134	delete fitvy;*/
2135
2136	if ( IsNaN(vx) \|\|IsNaN(vy) \|\| IsNaN(vz) ) {
2137	Msg(EsafMsg::Warning) << "AE1() fit returns a NaN, skipping" << MsgDispatch;
2138	return;
2139	}
2140
2141	// use the following only to debug (plot angles alphai vs time)
2142	if( fDoGraphUseShape ){
2143	Double_t tg[pointsAE1],xg[pointsAE1],yg[pointsAE1],zg[pointsAE1];
2144	Double_t xr[pointsAE1],yr[pointsAE1],zr[pointsAE1];
2145	for( Int_t i=0; i<pointsAE1; i++ ) {
2146	tg[i] = timestart[i];
2147	xg[i] = x[i]; yg[i]=y[i]; zg[i]=z[i];
2148	xr[i] = timestart[i]*vx + qx;
2149	yr[i] = timestart[i]*vy + qy;
2150	zr[i] = timestart[i]*vz + qz;
2151	}
2152	TGraph *grx = new TGraph(pointsAE1,tg,xg);
2153	grx->SetNameTitle("grx","x vs time");
2154	grx->SetMarkerSize(.5); grx->SetMarkerStyle(23);
2155	TGraph *gry = new TGraph(pointsAE1,tg,yg);
2156	gry->SetNameTitle("gry","y vs time");
2157	gry->SetMarkerSize(.5); gry->SetMarkerStyle(23);
2158	TGraph *grz = new TGraph(pointsAE1,tg,zg);
2159	grz->SetNameTitle("grz","z vs time");
2160	grz->SetMarkerSize(.5); grz->SetMarkerStyle(23);
2161
2162	TGraph *linex = new TGraph(pointsAE1,tg,xr); linex->SetLineColor(kRed);
2163	TGraph *liney = new TGraph(pointsAE1,tg,yr); liney->SetLineColor(kRed);
2164	TGraph *linez = new TGraph(pointsAE1,tg,zr); linez->SetLineColor(kRed);
2165	TCanvas *cAE1=new TCanvas("cAE1","cAE1",1); cAE1->Divide(2,2);
2166	cAE1->cd(1);
2167	grx->Draw("AP");
2168	linex->Draw("PLsame");
2169	cAE1->cd(2);
2170	gry->Draw("AP");
2171	liney->Draw("PLsame");
2172	cAE1->cd(3);
2173	grz->Draw("AP");
2174	linez->Draw("PLsame");
2175	cAE1->Write();
2176	delete cAE1;
2177	delete grx; delete gry; delete grz;
2178	delete linex; delete liney; delete linez;
2179	}
2180
2181	TVector3 rumenta(vx,vy,vz); rumenta.SetMag(1.);
2182	fTHETAreco = rumenta.Theta(); fPHIreco = rumenta.Phi();
2183
2184	fEASDir.SetMagThetaPhi(1,fTHETAreco,fPHIreco+Pi()); // main euso system
2185	fTHETAreco = fEASDir.Theta();
2186	fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*Pi()) );
2187	fDeltaTheta = fTHETAreco - fTrueTheta;
2188	fDeltaPhi = fPHIreco - fTruePhi;
2189	fDeltaEASDir = fEASDir.Angle(fTrueDir);
2190
2191	fHmax = FindHmax(fTHETAreco);
2192	Msg(EsafMsg::Info) << "method AE1(): fDeltaEASDir = " <<fDeltaEASDirRadToDeg() << " deg ( dT = " << fDeltaThetaRadToDeg() << " , dP = " << fDeltaPhi*RadToDeg() << " )" << MsgDispatch;
2193
2194	}
2195	//______________________________________________________________________________
2196	void TrackDirection2Module::All() {
2197	//
2198	// Execute all methods for reconstructing the shower direction.
2199	//
2200
2201	Msg(EsafMsg::Info) << "Using all methods .. " << MsgDispatch;
2202
2203	AA1(); if (fAA1nan) return;
2204	fTHETArecoAA1 = fTHETAreco;
2205	fPHIrecoAA1 = fPHIreco;
2206	fDeltaEASDirAA1 = fDeltaEASDir;
2207	AA2();
2208	fTHETArecoAA2 = fTHETAreco;
2209	fPHIrecoAA2 = fPHIreco;
2210	fDeltaEASDirAA2 = fDeltaEASDir;
2211	NE1();
2212	fTHETArecoNE1 = fTHETAreco;
2213	fPHIrecoNE1 = fPHIreco;
2214	fDeltaEASDirNE1 = fDeltaEASDir;
2215	NE2();
2216	fTHETArecoNE2 = fTHETAreco;
2217	fPHIrecoNE2 = fPHIreco;
2218	fDeltaEASDirNE2 = fDeltaEASDir;
2219	AE1();
2220	fTHETArecoAE1 = fTHETAreco;
2221	fPHIrecoAE1 = fPHIreco;
2222	fDeltaEASDirAE1 = fDeltaEASDir;
2223
2224	//if all is chosen are saved only AA1 results
2225	fTHETAreco = fTHETArecoAA1;
2226	fPHIreco = fPHIrecoAA1;
2227	fDeltaTheta = fTHETAreco - fTrueTheta;
2228	fDeltaPhi = fPHIreco - fTruePhi;
2229	fEASDir.SetMagThetaPhi(1., fTHETAreco, fPHIreco);
2230	fDeltaEASDir = fEASDir.Angle(fTrueDir);
2231
2232	}
2233
2234	//______________________________________________________________________________
2235	void TrackDirection2Module::CalculatefromBetaEASdir() {
2236	//
2237	// Calculate vector of EAS direction using angle fBeta
2238	// and the equation of TrackDirectionPlane
2239	//
2240
2241	fEASDir = Cos(fBeta) * fW + Sin(fBeta) * fU;
2242	fEASDir.SetPhi(fEASDir.Phi()+Pi()); // main euso system
2243	fTHETAreco = fEASDir.Theta();
2244	fPHIreco = ( fEASDir.Phi()>0 ? fEASDir.Phi() : fEASDir.Phi()+(2*Pi()) );
2245	fDeltaTheta = fTHETAreco - fTrueTheta;
2246	fDeltaPhi = fPHIreco - fTruePhi;
2247	fDeltaEASDir = fEASDir.Angle(fTrueDir);
2248	}
2249
2250	//______________________________________________________________________________
2251	Double_t TrackDirection2Module::FindHmax( Double_t theta ) {
2252	//
2253	// Find Hmax with the Linsley parametrization of the atmosphere
2254	// depending on the zenith angle of the shower and the value of Xmax
2255	// (fixed value is only a first approximation).
2256	//
2257
2258	Double_t Xmax = 831; //g/cm^2
2259	Double_t thetaloc = theta; //it should be zenith angle in local reference frame, to improve.
2260	Double_t Xv = Xmax*Cos(thetaloc);
2261
2262	Double_t X0 = 1036.1; // g/cm^2
2263	Double_t X4 = 631.1; // g/cm^2
2264	Double_t X10 = 271.1; // g/cm^2
2265	Double_t al(0), bl(0), cl(0);
2266	Double_t val(0);
2267	if( Xv<X0 && Xv>X4 ){
2268	cl = 9.9418638; // km
2269	bl = 1222.6562; // g/cm2
2270	al = -186.5562; // g/cm2
2271	} else if( Xv<X4 && Xv>X10 ) {
2272	cl = 8.7815355; // km
2273	bl = 1144.9069; // g/cm2
2274	al = -94.9199; // g/cm2
2275	} else if( Xv<X10 ) {
2276	cl = 6.3614304; // km
2277	bl = 1305.5948; // g/cm2
2278	al = 0.61289; // g/cm2
2279	} else {
2280	Msg(EsafMsg::Warning) <<"...something was wrong in TrackDirection2Module::FindHmax"<< MsgDispatch;
2281	return 0;
2282	}
2283	val = -cl * Log((Xmax * Cos(thetaloc)-al)/bl);
2284	return val;
2285	}
2286
2287	//______________________________________________________________________________
2288	Double_t TrackDirection2Module::CalculateRmax( Double_t hmax ) {
2289	//
2290	// Method to geometrically find Rmax using Hmax and the versor pointing
2291	// to the maximum of the shower
2292	//
2293
2294	Double_t thmax = fData.fMedDir.Theta();
2295	Double_t rmax = (EarthRadius()/km + fHISS)*Cos(thmax) - Sqrt(Power(EarthRadius()/km+hmax,2)
2296	- Power((EarthRadius()/km+fHISS)*Sin(thmax),2));
2297	return rmax;
2298	}
2299
2300	//______________________________________________________________________________
2301	Double_t TrackDirection2Module::DeltaX( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
2302	//
2303	// Calculate error on X projection on the focal surface of a given point
2304	// on the unitary sphere
2305	//
2306
2307	return Sqrt( Power((Cos(phi)Cos(theta)+Sin(theta))errtheta/(Cos(theta)*Cos(theta)),2.) +
2308	Power(Sin(theta)Sin(phi)errphi/Cos(theta),2.));
2309	}
2310
2311	//______________________________________________________________________________
2312	Double_t TrackDirection2Module::DeltaY( Double_t theta, Double_t phi, Double_t errtheta, Double_t errphi ) {
2313	//
2314	// Calculate error on Y projection on the focal surface of a given point
2315	// on the unitary sphere
2316	//
2317
2318	return Sqrt( Power((Sin(phi)Cos(theta)+Sin(theta))errtheta/(Cos(theta)*Cos(theta)),2.) +
2319	Power(Sin(theta)Cos(phi)errphi/Cos(theta),2.));
2320	}
2321
2322	//___________________________________________________
2323	void ContainerData::Clear(Option_t*) {
2324	//
2325	// Clear of the container data class
2326	//
2327
2328	fPixXYZ.clear();
2329	fPos.clear();
2330	fErr.clear();
2331	fSpPos.clear();
2332	fSpErr.clear();
2333	fSpPosSel.clear();
2334	fTime.clear();
2335	fTimeSel.clear();
2336	fAlpha.clear();
2337	fSigmaPhi.clear();
2338	fSigmaTheta.clear();
2339	fSigmaPhiSel.clear();
2340	fSigmaThetaSel.clear();
2341	fHits.clear();
2342	fHitsSel.clear();
2343	fNumPoints = 0;
2344	fNumHits = 0;
2345	fCentroid.SetXYZ(0,0,0);
2346	fMedDir.SetXYZ(0,0,0);
2347	fVectorRmax.SetXYZ(0,0,0);
2348	fMedTime = 0;
2349	fMedAlpha = 0;
2350	fApparentTimeLenght = 0;
2351	fGtuLength = 0;
2352	}
2353	//___________________________________________________
2354	void ContainerData::ClearMemory() {
2355	//
2356	// Clear of the container data clas and free the memory held in the buffers
2357	//
2358
2359	vector<TVector3> dummyV1; dummyV1.swap(fPixXYZ);
2360	vector<TVector3> dummyV2; dummyV2.swap(fPos);
2361	vector<TVector3> dummyV3; dummyV3.swap(fErr);
2362	vector<TVector3> dummyV4; dummyV4.swap(fSpPos);
2363	vector<TVector3> dummyV5; dummyV5.swap(fSpErr);
2364	vector<TVector3> dummyV6; dummyV6.swap(fSpPosSel);
2365
2366	vector<Double_t> dummyD1; dummyD1.swap(fTime);
2367	vector<Double_t> dummyD2; dummyD2.swap(fTimeSel);
2368	vector<Double_t> dummyD3; dummyD3.swap(fAlpha);
2369	vector<Double_t> dummyD4; dummyD4.swap(fSigmaPhi);
2370	vector<Double_t> dummyD5; dummyD5.swap(fSigmaTheta);
2371	vector<Double_t> dummyD6; dummyD6.swap(fSigmaPhiSel);
2372	vector<Double_t> dummyD7; dummyD7.swap(fSigmaThetaSel);
2373
2374	vector<Int_t> dummyI1; dummyI1.swap(fHits);
2375	vector<Int_t> dummyI2; dummyI2.swap(fHitsSel);
2376	}

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

Context Navigation

source: JEM-EUSO/esaf_lal/branches/camille/packages/reconstruction/modules/shower/fitting/src/TrackDirection2Module.cc @ 181

Download in other formats: