source: branches/Mars_use_drstimefiles/mtemp/mifae/macros/OptimizeDisp.C@ 19727

Last change on this file since 19727 was 5920, checked in by domingo, 20 years ago
*** empty log message ***
File size: 18.9 KB
Line 
1//************************************************************************
2//
3// Authors : Eva Domingo, 12/2004 <mailto:domingo@ifae.es>
4// Wolfgang Wittek, 12/2004 <mailto:wittek@mpppmu.mpg.de>
5//
6//
7// Macro for estimating the DISP parameter
8// ---------------------------------------
9//
10// DISP is the distance between the center of gravity of the shower image
11// and the estimated source position, assumed to lie on the major
12// axis of the shower
13//
14// In order to get an estimate of DISP
15// - one assumes a parametrization DISP = f(p[0],p[1],p[2],... ; s1,s2,s3,...)
16// where s1,s2,s3,... are measurable quantities like
17// image parameters, (theta, phi), etc.
18// and p[0],p[1],p[2], ... are free parameters
19//
20// - and determines the free parameters p[0],p[1],p[2]... from MC gamma data
21// by minimizing a parameter which measures the deviation of the estimated
22// source position (using DISP) and the true source position
23//
24// The following classes are used :
25//
26// MDispParameters container holding the parameters p[0],p[1],p[2],...
27//
28// MDispCalc task calculating DISP with the parameters stored in
29// MDispParameters
30//
31// ::Calc member function where the DISP parameterization is
32// defined. It calculates DISP for a given event;
33// it is called by MDispCalc::Process()
34//
35// MFindDisp class with several functions :
36//
37// ::DefineTrainMatrix \ member functions which generate the training
38// ::DefineTestMatrix | and test samples (in the form of matrices)
39// ::DefineTrainTestMatrix / from the MC gamma data
40//
41// ::FindParams member function steering the minimization
42// (definition of the fcn function for Minuit, setting up the
43// event loop to be executed in each minimization step,
44// call to Minuit);
45// for the minimization the training matrices are used
46//
47// ::TestParams member function testing the quality of the DISP estimate;
48// for the test the test matrices are used
49//
50// MHDisp container for histograms which are useful for judging
51// the quality of the DISP estimate
52// Also task calculating the Minimization parameter
53// contribution of an event (in Fill())
54// and computing the final Minimization parameter of the
55// whole event sample (in Finalize())
56//
57// MFDisp filter to select sample for the optimization
58//
59//
60// The classes are stored in the CVS directory Mars/mtemp/mifae/library
61// and Mars/mtemp/mifae/macros
62//
63//
64//************************************************************************
65
66void OptimizeDisp()
67{
68 //************************************************************************
69
70 Bool_t CMatrix = kFALSE; // Create training and test matrices
71 Bool_t WOptimize = kFALSE; // Do optimization using the training sample
72 // and write Disp parameter values
73 // onto the file parDispfile
74 // Optimize Disp with :
75 //TString typeOpt = "Data";
76 TString typeOpt = "MC";
77 Bool_t RTest = kFALSE; // Test the quality of the Disp estimate
78 // using the test matrix
79 Bool_t WDisp = kTRUE; // Make Disp plots for the data of type :
80 //TString typeInput = "ON";
81 //TString typeInput = "OFF";
82 TString typeInput = "MC";
83
84 //************************************************************************
85
86 gLog.SetNoColors();
87
88 if (gRandom)
89 delete gRandom;
90 gRandom = new TRandom3(0);
91
92 //-----------------------------------------------
93 //names of files to be read for generating the training and test matrices
94 const char *filetrain = "...";
95 const char *filetest = "...";
96
97 //-----------------------------------------------
98 // path for input for Mars
99 TString inPathON =
100 "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_16/";
101 TString inPathOFF =
102 "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_17/";
103 TString inPathMC =
104 //"/discpepe/root_0.73mirror/wuerzburg/";
105 "~domingo/MAGIC/Disp05/";
106
107 // path for output from Mars
108 TString outPath = "~domingo/MAGIC/Disp05/DispOptimization/";
109
110 //-----------------------------------------------
111 // names of files for which Disp plots should be made
112 const char *offfile = "*OffCrabLow42*";
113 const char *onfile = "*CrabLowEn42*";
114 const char *mcfile = "All_gammas0.73wuerzburg_zbin0to12_cleaned_3020_PeakSearch6_RunHeaders";
115 //-----------------------------------------------
116
117
118 //---------------------------------------------------------------------
119 gLog << "=====================================================" << endl;
120 gLog << "Macro OptimizeDisp : Start " << endl;
121
122 gLog << "" << endl;
123 gLog << "Macro OptimizeDisp : CMatrix, WOptimize, RTest, WDisp = "
124 << (CMatrix ? "kTRUE" : "kFALSE") << ", "
125 << (WOptimize ? "kTRUE" : "kFALSE") << ", "
126 << (RTest ? "kTRUE" : "kFALSE") << ", "
127 << (WDisp ? "kTRUE" : "kFALSE") << endl;
128
129
130 //--------------------------------------------
131 // files to contain the matrices (generated from filenameTrain and
132 // filenameTest)
133 //
134 // for the training
135 TString fileMatrixTrain = outPath;
136 fileMatrixTrain += "MatrixTrainDisp";
137 fileMatrixTrain += ".root";
138 gLog << "" << endl;
139 gLog << "Files containing the Training and Test matrices :" << endl;
140 gLog << " fileMatrixTrain = " << fileMatrixTrain << endl;
141
142 // for testing
143 TString fileMatrixTest = outPath;
144 fileMatrixTest += "MatrixTestDisp";
145 fileMatrixTest += ".root";
146 gLog << " fileMatrixTest = " << fileMatrixTest << endl;
147 gLog << "" << endl;
148
149
150 //---------------
151 // file to contain the optimum Disp parameter values
152 TString parDispFile = outPath;
153 // parDispFile += "parDispMC.root";
154 parDispFile += "parDispMC_withcuts.root";
155
156 gLog << "" << endl;
157 gLog << "File containing the optimum Disp parameter values : parDispFile = "
158 << parDispFile << endl;
159
160
161 // Set the filter cuts to select a sample of events for the Disp optimization
162 // (islandsmin,islandsmax,usedpixelsmin,usedpixelsmax,corepixelsmin,
163 // corepixelsmax,sizemin,sizemax,leakage1min,leakage1max,leakage2min,
164 // leakage2max,lengthmin,widthmin);
165 MFDisp *fdisp = NULL;
166 fdisp = new MFDisp;
167 fdisp->SetCuts(0,2,7,600,0,600,0.,3000.,0.,0.,0.,0.,0.,0.);
168 fdisp->SetName("FilterSelector2");
169
170
171 // Create the MFindDisp object and set the file to store optimium Disp parameters
172 MFindDisp finddisp(fdisp);
173 finddisp.SetFilenameParam(parDispFile);
174
175
176
177 //======================================================================
178 // Create matrices and write them onto files
179 //======================================================================
180
181 if (CMatrix)
182 {
183 gLog << "-----------------------------------------------------" << endl;
184 gLog << "Generate the training and test samples" << endl;
185
186 //--------------------------------------------
187 // files to be read for optimizing the Disp parameters
188 //
189 // for the training
190 TString filenameTrain = inPathMC;
191 filenameTrain += mcfile;
192 filenameTrain += ".root";
193 Int_t howManyTrain = 30000;
194 gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = "
195 << howManyTrain << endl;
196
197 // for testing
198 TString filenameTest = inPathMC;
199 filenameTest += mcfile;
200 filenameTest += ".root";
201 Int_t howManyTest = 30000;
202 gLog << "filenameTest = " << filenameTest << ", howManyTest = "
203 << howManyTest << endl;
204
205
206 //--------------------------------------------
207 // For the events in the matrices define requested distribution
208 // in some quantities; this may be a 1-dim, 2-dim or 3-dim distribution
209
210 /*
211 // Define the cos(theta) distribution
212 TString mname("costheta");
213 MBinning bincost("Binning"+mname);
214 bincost.SetEdges(10, 0., 1.0);
215
216 //MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)");
217 MH3 mh3("cos(MMcEvt.fTelescopeTheta)");
218 mh3.SetName(mname);
219 MH::SetBinning(&mh3.GetHist(), &bincost);
220
221 for (Int_t i=1; i<=mh3.GetNbins(); i++)
222 mh3.GetHist().SetBinContent(i, 1.0);
223 */
224
225 // Define the log10(Etrue) distribution
226 TString mh3name("log10Etrue");
227 MBinning binslogE("Binning"+mh3name);
228 binslogE.SetEdges(50, 1., 3.);
229
230 MH3 mh3("log10(MMcEvt.fEnergy)");
231 mh3.SetName(mh3name);
232 MH::SetBinning(&mh3.GetHist(), &binslogE);
233
234 // Flat energy distribution
235 // for (Int_t i=1; i<=mh3.GetNbins(); i++)
236 // mh3.GetHist().SetBinContent(i, 1.0);
237
238 // Power law energy distribution
239 // for (Int_t i=1; i<=mh3.GetNbins(); i++)
240 // {
241 // Double_t weight = pow((Double_t)i, -1.7);
242 // mh3.GetHist().SetBinContent(i, weight);
243 // }
244
245 /*
246 // define the 'cos(Theta) vs. log10(Etrue)' distribution
247 TString mh3name("cosThetaVslog10Etrue");
248 MBinning binslogE("Binning"+mh3name+"X");
249 binslogE.SetEdges(18, 1.5, 2.2);
250 MBinning binscost("Binning"+mh3name+"Y");
251 binscost.SetEdges(10, 0., 1.0);
252
253 // MH3 mh3("log10(MMcEvt.fEnergy)", "cos((MPointingPos.fZd)/kRad2Deg)");
254 MH3 mh3("log10(MMcEvt.fEnergy)", "cos(MMcEvt.fTelescopeTheta)");
255 mh3.SetName(mh3name);
256 MH::SetBinning((TH2*)&mh3.GetHist(), &binslogE, &binscost);
257
258 for (Int_t i=1; i<=((TH2*)mh3.GetHist()).GetNbinsX(); i++)
259 {
260 // Double_t weight = pow((Double_t)i, -1.7);
261 for (Int_t j=1; j<=((TH2*)mh3.GetHist()).GetNbinsY(); j++)
262 {
263 // mh3.GetHist().SetBinContent(i, j, weight);
264 mh3.GetHist().SetBinContent(i, j, 1.);
265 }
266 }
267 */
268
269 //--------------------------
270 // Training and test samples are generated from the same input files
271 if (filenameTrain == filenameTest)
272 {
273 if ( !finddisp.DefineTrainTestMatrix(
274 filenameTrain, mh3,
275 howManyTrain, howManyTest,
276 fileMatrixTrain, fileMatrixTest, 0) )
277 {
278 *fLog << "OptimizeDisp.C : DefineTrainTestMatrix failed" << endl;
279 return;
280 }
281 }
282
283 //--------------------------
284 // Training and test samples are generated from different input files
285 else
286 {
287 if ( !finddisp.DefineTrainMatrix(filenameTrain, mh3,
288 howManyTrain, fileMatrixTrain, 0) )
289 {
290 *fLog << "OptimizeDisp.C : DefineTrainMatrix failed" << endl;
291 return;
292 }
293
294 if ( !finddisp.DefineTestMatrix( filenameTest, mh3,
295 howManyTest, fileMatrixTest, 0) )
296 {
297 *fLog << "OptimizeDisp.C : DefineTestMatrix failed" << endl;
298 return;
299 }
300 }
301
302 gLog << "Generation of training and test samples finished" << endl;
303 gLog << "-----------------------------------------------------" << endl;
304 }
305
306
307
308 //======================================================================
309 // Optimize Disp parameters using the training sample
310 //
311 // the initial values of the parameters are taken
312 // - from the file parDispInit (if != "")
313 // - or from the arrays params and steps (if their sizes are != 0)
314 // - or from the MDispParameters constructor
315 //======================================================================
316
317 if (WOptimize)
318 {
319 gLog << "-----------------------------------------------------" << endl;
320 gLog << "Optimize the Disp parameters over a "
321 << typeOpt << " sample, using the training matrices" << endl;
322
323 // Read training matrices from file
324 finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest);
325
326 //--------------------------------------------
327 // file which contains the initial values for the Disp parameters;
328 // if parDispInit =""
329 // - the initial values are taken from the arrays params and steps
330 // if their sizes are different from 0
331 // - otherwise they are taken from the MDispParameters constructor
332
333 TString parDispInit = outPath;
334 //parDispInit += "parDispInit.root";
335 parDispInit = "";
336
337 //--------------------------------------------
338
339 TArrayD params(0);
340 TArrayD steps(0);
341
342 if (parDispInit == "")
343 {
344 Double_t vparams[5] = {
345 // p[0], p[1], p[2], p[3], p[4], p[5]
346 1.0, 0.6, -0.8, -0.8, -1.2/*, 0.5*/};
347 // 0.5, 0., 0.03, 0., 0./*, 0.5*/};
348 // 0.8, 0., 0., 0., 0./*, 0.5*/};
349
350 Double_t vsteps[5] = {
351 // dp[0], dp[1], dp[2], dp[3], dp[4], dp[5]
352 0.01, 0.005, 0.005, 0.005, 0.01/*, 0.001*/};
353 // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/};
354 // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/};
355
356 params.Set(5, vparams);
357 steps.Set (5, vsteps );
358 }
359
360
361 Bool_t rf;
362 rf = finddisp.FindParams(parDispInit, params, steps);
363
364 if (!rf)
365 {
366 gLog << "OptimizeDisp.C : optimization of the Disp parameters failed" << endl;
367 return;
368 }
369
370 gLog << "Optimization of Disp parameters finished" << endl;
371 gLog << "-----------------------------------------------------" << endl;
372 }
373
374
375
376 //======================================================================
377 // Test the Disp parameters on the test sample
378 //======================================================================
379
380 if (RTest)
381 {
382 // Read test matrices from file
383 finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest);
384
385 gLog << "-----------------------------------------------------" << endl;
386 gLog << "Test the Disp parameters using the test matrices" << endl;
387 Bool_t rt = finddisp.TestParams();
388 if (!rt)
389 {
390 gLog << "Test of the Disp parameters on the test matrices failed"
391 << endl;
392 }
393 gLog << "Test of the Disp parameters finished" << endl;
394 gLog << "-----------------------------------------------------" << endl;
395 }
396
397
398
399
400 //======================================================================
401 // Make Disp plots
402 //======================================================================
403
404 if (WDisp)
405 {
406 gLog << "" << endl;
407 gLog << "-----------------------------------------------------" << endl;
408 gLog << "Make plots for Disp for data of the type " << typeInput << endl;
409
410 //--------------------------------------------
411 // type of data to be read (ON, OFF or MC)
412 if (typeInput == "ON")
413 TString file(onfile);
414 else if (typeInput == "OFF")
415 TString file(offfile);
416 else if (typeInput == "MC")
417 TString file(mcfile);
418
419 // name of input root file
420 TString filenameData = inPathMC;
421 filenameData += mcfile;
422 filenameData += ".root";
423 gLog << "Input file '" << filenameData << endl;
424
425 //----------------------------------------------------
426 // read in optimum Disp parameter values
427
428 TFile inparam(parDispFile);
429 MDispParameters dispin;
430 dispin.Read("MDispParameters");
431 inparam.Close();
432
433 gLog << "Disp parameter values were read in from file '"
434 << parDispFile << "'" << endl;
435
436 TArrayD dispPar;
437 dispPar = dispin.GetParameters();
438
439 TArrayD dispStep;
440 dispStep = dispin.GetStepsizes();
441
442 gLog << "optimum parameter values for calculating Disp : " << endl;
443 for (Int_t i=0; i<dispPar.GetSize(); i++)
444 {
445 gLog << dispPar[i] << ", ";
446 }
447 gLog << endl;
448
449
450 //----------------------------------------------------
451 MTaskList tliston;
452 MParList pliston;
453
454 MReadMarsFile read("Events", filenameData);
455 read.DisableAutoScheme();
456
457 // set cuts to select an event sample to apply Disp
458 MContinue contdisp(fdisp);
459
460 // calculate Disp
461 MDispCalc dispcalc;
462
463 // make Disp plots
464 // SelectedPos = 1 means choose the right source position
465 // 2 wrong
466 // 3 the position according to M3Long
467 // 4 the position according to Asym
468 MHDisp hdisp1;
469 hdisp1.SetName("MHDispCorr");
470 hdisp1.SetSelectedPos(1);
471 MFillH filldisp1("MHDispCorr[MHDisp]", "");
472
473 MHDisp hdisp2;
474 hdisp2.SetName("MHDispWrong");
475 hdisp2.SetSelectedPos(2);
476 MFillH filldisp2("MHDispWrong[MHDisp]", "");
477
478 MHDisp hdisp3;
479 hdisp3.SetName("MHDispM3Long");
480 hdisp3.SetSelectedPos(3);
481 MFillH filldisp3("MHDispM3Long[MHDisp]", "");
482
483 MHDisp hdisp4;
484 hdisp4.SetName("MHDispAsym");
485 hdisp4.SetSelectedPos(4);
486 MFillH filldisp4("MHDispAsym[MHDisp]", "");
487
488
489 //*****************************
490 // entries in MParList
491
492 pliston.AddToList(&tliston);
493 pliston.AddToList(&dispin);
494 pliston.AddToList(&hdisp1);
495 pliston.AddToList(&hdisp2);
496 pliston.AddToList(&hdisp3);
497 pliston.AddToList(&hdisp4);
498
499 //*****************************
500 // entries in MTaskList
501
502 tliston.AddToList(&read);
503 if (fdisp != NULL)
504 tliston.AddToList(&contdisp);
505 tliston.AddToList(&dispcalc);
506 tliston.AddToList(&filldisp1);
507 tliston.AddToList(&filldisp2);
508 tliston.AddToList(&filldisp3);
509 tliston.AddToList(&filldisp4);
510
511 //*****************************
512
513 //-------------------------------------------
514 // Execute event loop
515 //
516 MProgressBar bar;
517 MEvtLoop evtloop;
518 evtloop.SetParList(&pliston);
519 evtloop.SetProgressBar(&bar);
520
521 Int_t maxevents = -1;
522 if ( !evtloop.Eventloop(maxevents) )
523 return;
524
525 tliston.PrintStatistics(0, kTRUE);
526
527 //-------------------------------------------
528 // Display the histograms
529 //
530 hdisp1.Draw();
531 hdisp2.Draw();
532 hdisp3.Draw();
533 hdisp4.Draw();
534
535 //-------------------------------------------
536
537 gLog << "Disp plots were made for file '" << filenameData << endl;
538 gLog << "-----------------------------------------------------" << endl;
539 }
540
541 delete fdisp;
542
543 gLog << "Macro OptimizeDisp.C : End " << endl;
544 gLog << "=======================================================" << endl;
545
546}
547
548
549
550
551
552
553
554
555
Note: See TracBrowser for help on using the repository browser.