Changeset 17914 for trunk/Mars/fact
- Timestamp:
- 07/15/14 15:44:29 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/fact/analysis/mc/callisto.C
r17870 r17914 32 32 #include "MHSectorVsTime.h" 33 33 #include "MHCamEvent.h" 34 #include "MExtract TimeAndChargeSpline.h"34 #include "MExtractFACT.h" 35 35 #include "MFillH.h" 36 36 #include "MDrsCalibApply.h" … … 75 75 */ 76 76 77 int callisto(const TString drsfile="test300samples .drs.fits",77 int callisto(const TString drsfile="test300samples2.drs.fits.gz", 78 78 const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits", 79 79 const TString datfile="00000003.387_D_MonteCarlo010_Events.fits", 80 TString out file= "",80 TString outpath = "", 81 81 TString displayfile = "", TString displaytitle = "") 82 82 { … … 88 88 89 89 FileStat_t fstat; 90 int rc = gSystem->GetPathInfo(out file, fstat);90 int rc = gSystem->GetPathInfo(outpath, fstat); 91 91 bool isdir = !rc || R_ISDIR(fstat.fMode); 92 92 93 const char *buf = gSystem->ConcatFileName(outfile, "callisto.root"); 94 outfile = buf; 93 TString filename = datfile + "_callisto.root"; 94 filename.Replace(0, filename.Last('/')+1, ""); 95 const char *buf = gSystem->ConcatFileName(outpath, filename); 96 TString outfile = buf; 95 97 delete [] buf; 96 98 … … 119 121 120 122 // map file to use (get that from La Palma!) 121 const char *pmap = usemap ? "TestForThomas/FACT/FACTmap111030.txt" : NULL; 122 123 Bool_t maximum = kTRUE; 124 125 //const char *lp_template = maximum ? 126 // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" : 127 // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root"; 128 129 const char *pulse_template = "TestForThomas/FACT/template-pulse.root"; 123 const char *pmap = usemap ? "/home/isdc/toscanos/FACT/Mars_svn/resources/FACTmap111030.txt" : NULL; 130 124 131 125 // ------------------------------------------------------ 132 133 // Calib: 51 / 90 / 197 (20% TH)134 // Data: 52 / 64 / 104 (20% TH)135 136 // Extraction range in slices. It will always(!) contain the full range137 // of integration138 const int first_slice = 20; // 10ns139 const int last_slice = 250; // 125ns140 141 // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:142 //143 // IntegralFixed:144 // * fRiseTime: Number of slices left from arrival time145 // * fFallTime: Number of slices right from arrival time146 // IntegralRel:147 // * fRiseTime: Number of slices left from maximum time148 // * fFallTime: Number of slices right from maximum time149 //150 const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns151 const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns152 153 const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns154 const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns155 156 // Extraction type: Extract integral and half leading edge157 158 const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);159 //const int type = MExtralgoSpline::kIntegralFixed;160 161 162 const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }163 164 Long_t max = 0; // All165 Long_t max3 = max; // Pedestal Rndm166 Long_t max4 = max; // Pedestal Ext167 168 // ======================================================169 170 if (pmap && gSystem->AccessPathName(pmap, kFileExists))171 {172 gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;173 return 1;174 }175 176 gLog.Separator("Callisto");177 gLog << all;178 gLog << "Data File: " << datfile << '\n';179 gLog << "DRS calib 300: " << drsfile << endl;;180 181 MDrsCalibration drscalib300;182 if (!drscalib300.ReadFits(drsfile.Data())) {183 gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;184 return 5;185 }186 gLog << all;187 gLog << "Pedestal file: " << pedfile << '\n';188 gLog << "Output file: " << outfile << '\n';189 gLog << "Display file: " << displayfile << '\n';190 gLog << "Display title: " << displaytitle << endl;191 192 // ------------------------------------------------------193 MStatusArray arrt, arrp;194 195 // TFile ft(lp_template);196 // if (arrt.Read()<=0)197 // {198 // gLog << err << "ERROR - Reading LP template from " << lp_template << endl;199 // return 100;200 // }201 202 // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");203 // if (!lpref)204 // {205 // gLog << err << "ERROR - LP Template not found in " << lp_template << endl;206 // return 101;207 // }208 // lpref->SetDirectory(0);209 210 // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");211 // if (!gain)212 // {213 // gLog << err << "ERROR - Gain not found in " << lp_template << endl;214 // return 101;215 // }216 // gain->SetDirectory(0);217 218 TFile fp(pulse_template);219 if (arrp.Read()<=0)220 {221 gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;222 return 102;223 }224 225 TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");226 if (!hpulse)227 {228 gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;229 return 103;230 }231 hpulse->SetDirectory(0);232 // ======================================================233 126 234 127 MStatusDisplay *d = new MStatusDisplay; … … 236 129 MBadPixelsCam badpixels; 237 130 badpixels.InitSize(1440); 238 /*239 131 badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); 240 132 badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); … … 243 135 badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); 244 136 badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); 245 */246 137 // Twin pixel 247 138 // 113 … … 251 142 // 1195 252 143 // 1393 144 145 // ------------------------------------------------------ 146 147 // Calib: 51 / 90 / 197 (20% TH) 148 // Data: 52 / 64 / 104 (20% TH) 149 150 // Extraction range in slices. It will always(!) contain the full range 151 // of integration 152 const int first_slice = 25; // 10ns 153 const int last_slice = 225; // 125ns 154 155 Long_t max = 0; // All 156 Long_t max3 = max; // Pedestal Rndm 157 Long_t max4 = max; // Pedestal Ext 158 159 // ====================================================== 160 161 //double scale = 0.1; 162 double scale = 0.1024; 163 164 // ====================================================== 165 166 if (pmap && gSystem->AccessPathName(pmap, kFileExists)) 167 { 168 gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl; 169 return 1; 170 } 171 172 gLog.Separator("Callisto"); 173 gLog << all; 174 gLog << "Data File: " << datfile << '\n'; 175 gLog << "DRS calib 300: " << drsfile << endl;; 176 177 MDrsCalibration drscalib300; 178 if (!drscalib300.ReadFits(drsfile.Data())) { 179 gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; 180 return 5; 181 } 182 gLog << all; 183 gLog << "Pedestal file: " << pedfile << '\n'; 184 gLog << "Output file: " << outfile << '\n'; 185 gLog << "Display file: " << displayfile << '\n'; 186 gLog << "Display title: " << displaytitle << endl; 187 188 // ------------------------------------------------------ 189 /* 190 MStatusArray arrt, arrp; 191 192 TFile ft(lp_template); 193 if (arrt.Read()<=0) 194 { 195 gLog << err << "ERROR - Reading LP template from " << lp_template << endl; 196 return 100; 197 } 198 199 MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); 200 if (!lpref) 201 { 202 gLog << err << "ERROR - LP Template not found in " << lp_template << endl; 203 return 101; 204 } 205 lpref->SetDirectory(0); 206 207 MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); 208 if (!gain) 209 { 210 gLog << err << "ERROR - Gain not found in " << lp_template << endl; 211 return 101; 212 } 213 gain->SetDirectory(0); 214 215 TFile fp(pulse_template); 216 if (arrp.Read()<=0) 217 { 218 gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; 219 return 102; 220 } 221 222 TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); 223 if (!hpulse) 224 { 225 gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; 226 return 103; 227 } 228 hpulse->SetDirectory(0); 229 */ 230 // ====================================================== 253 231 254 232 MDrsCalibrationTime timecam; … … 268 246 269 247 // ========================= Result ================================== 270 248 /* 271 249 //~ Double_t avgS = evt1f.GetHist()->GetMean(); 272 250 //~ Double_t medS = evt1f.GetHist()->GetMedian(); … … 340 318 //hcalco.Scale(scale); 341 319 hcalco.DrawCopy(); 342 320 */ 343 321 // ====================================================== 344 322 … … 372 350 MDrsCalibApply drsapply3; 373 351 352 MFilterData filterdata5; 353 374 354 //--- 375 355 376 MExtract TimeAndChargeSplineextractor3;356 MExtractFACT extractor3; 377 357 extractor3.SetRange(first_slice, last_slice); 378 extractor3.SetRiseTimeHiGain(rise_time_dat);379 extractor3.SetFallTimeHiGain(fall_time_dat);380 extractor3.SetHeightTm(heighttm);381 extractor3.SetChargeType(type);382 extractor3.SetSaturationLimit(600000);383 358 extractor3.SetNoiseCalculation(kTRUE); 384 385 // MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S");386 387 // MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm");388 389 // Use this for data, but not for calibration events390 // evt2a.SetErrorSpread(kFALSE);391 392 /*393 MCalibrateData conv3;394 conv3.SetCalibrationMode(MCalibrateData::kNone);395 conv3.SetPedestalFlag(MCalibrateData::kNo);396 conv3.SetCalibConvMinLimit(0);397 conv3.SetCalibConvMaxLimit(10000);398 conv3.SetScaleFactor(scale);399 */400 359 401 360 MCalibrateFact conv3; 402 361 conv3.SetScale(scale); 403 conv3.SetCalibConst(calib);362 //conv3.SetCalibConst(calib); 404 363 405 364 MBadPixelsTreat treat3; … … 420 379 tlist3.AddToList(&drsapply3); 421 380 tlist3.AddToList(&cont3); 381 tlist3.AddToList(&filterdata3); 422 382 tlist3.AddToList(&extractor3); 423 383 // tlist3.AddToList(&fill3a); … … 461 421 MDrsCalibApply drsapply4; 462 422 463 MExtractTimeAndChargeSpline extractor4; 423 MFilterData filterdata4; 424 425 MExtractFACT extractor4; 464 426 extractor4.SetRange(first_slice, last_slice); 465 extractor4.SetRiseTimeHiGain(rise_time_dat);466 extractor4.SetFallTimeHiGain(fall_time_dat);467 extractor4.SetHeightTm(heighttm);468 extractor4.SetChargeType(type);469 extractor4.SetSaturationLimit(600000);470 427 extractor4.SetNoiseCalculation(kFALSE); 471 428 472 // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S");473 474 // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt");475 476 // Use this for data, but not for calibration events477 // evt3a.SetErrorSpread(kFALSE);478 /*479 MCalibrateData conv4;480 conv4.SetCalibrationMode(MCalibrateData::kNone);481 conv4.SetPedestalFlag(MCalibrateData::kNo);482 conv4.SetCalibConvMinLimit(0);483 conv4.SetCalibConvMaxLimit(10000);484 conv4.SetScaleFactor(scale);485 */486 429 MCalibrateFact conv4; 487 430 conv4.SetScale(scale); 488 conv4.SetCalibConst(calib);431 //conv4.SetCalibConst(calib); 489 432 490 433 MBadPixelsTreat treat4; … … 504 447 tlist4.AddToList(&drsapply4); 505 448 tlist4.AddToList(&cont4); 449 tlist4.AddToList(&filterdata4); 506 450 tlist4.AddToList(&extractor4); 507 // tlist4.AddToList(&fill4a);508 451 tlist4.AddToList(&conv4); 509 452 tlist4.AddToList(&treat4); … … 550 493 MDrsCalibApply drsapply5; 551 494 495 MTreatSaturation treatsat5; 496 497 MFilterData filterdata5; 498 552 499 MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat"); 553 500 MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal"); … … 559 506 // --- 560 507 561 MExtract TimeAndChargeSplineextractor5dat;508 MExtractFACT extractor5dat; 562 509 extractor5dat.SetRange(first_slice, last_slice); 563 extractor5dat.SetRiseTimeHiGain(rise_time_dat);564 extractor5dat.SetFallTimeHiGain(fall_time_dat);565 extractor5dat.SetHeightTm(heighttm);566 extractor5dat.SetChargeType(type);567 extractor5dat.SetSaturationLimit(600000);568 510 extractor5dat.SetNoiseCalculation(kFALSE); 569 511 570 MExtract TimeAndChargeSplineextractor5cal;512 MExtractFACT extractor5cal; 571 513 extractor5cal.SetRange(first_slice, last_slice); 572 extractor5cal.SetRiseTimeHiGain(rise_time_cal);573 extractor5cal.SetFallTimeHiGain(fall_time_cal);574 extractor5cal.SetHeightTm(heighttm);575 extractor5cal.SetChargeType(type);576 extractor5cal.SetSaturationLimit(600000);577 514 extractor5cal.SetNoiseCalculation(kFALSE); 578 515 579 MExtract TimeAndChargeSplineextractor5tm("ExtractTM");516 MExtractFACT extractor5tm("ExtractTM"); 580 517 extractor5tm.SetRange(last_slice, 294); 581 extractor5tm.SetRiseTimeHiGain(1);582 extractor5tm.SetFallTimeHiGain(1);583 extractor5tm.SetHeightTm(0.5);584 extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);585 extractor5tm.SetSaturationLimit(600000);586 518 extractor5tm.SetNoiseCalculation(kFALSE); 587 519 extractor5tm.SetNameSignalCam("TimeMarkerAmplitude"); … … 593 525 594 526 // --- 595 /*596 MCalibrateData conv5;597 conv5.SetCalibrationMode(MCalibrateData::kNone);598 conv5.SetPedestalFlag(MCalibrateData::kNo);599 conv5.SetCalibConvMinLimit(0);600 conv5.SetCalibConvMaxLimit(10000);601 conv5.SetScaleFactor(scale);602 */603 527 MCalibrateFact conv5; 604 528 conv5.SetScale(scale); 605 conv5.SetCalibConst(calib);529 //conv5.SetCalibConst(calib); 606 530 607 531 MCalibrateDrsTimes calctm5; … … 689 613 tlist5.AddToList(&filterped); 690 614 tlist5.AddToList(&fill5a); 615 tlist5.AddToList(&treatsat5); 616 tlist5.AddToList(&filterdata5); 691 617 tlist5.AddToList(&extractor5dat); 692 618 tlist5.AddToList(&extractor5cal);
Note:
See TracChangeset
for help on using the changeset viewer.