Changeset 17999


Ignore:
Timestamp:
10/28/14 14:35:42 (10 years ago)
Author:
toscano
Message:
Callisto macro for MC with the saturation treatement and new calibration
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/MarsISDCBranchBasedOn17887/fact/analysis/mc/callisto.C

    r17870 r17999  
    3232#include "MHSectorVsTime.h"
    3333#include "MHCamEvent.h"
    34 #include "MExtractTimeAndChargeSpline.h"
     34#include "MExtractFACT.h"
    3535#include "MFillH.h"
    3636#include "MDrsCalibApply.h"
    … …  
    7575 */
    7676
    77 int callisto(const TString drsfile="test300samples.drs.fits",
     77int callisto(const TString drsfile="test300samples2.drs.fits.gz",
    7878             const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits",
    7979             const TString datfile="00000003.387_D_MonteCarlo010_Events.fits",
    80              TString outfile = "",
     80             TString outpath = "",
    8181             TString displayfile = "", TString displaytitle = "")
    8282{
    … …  
    8888
    8989    FileStat_t fstat;
    90     int rc = gSystem->GetPathInfo(outfile, fstat);
     90    int rc = gSystem->GetPathInfo(outpath, fstat);
    9191    bool isdir = !rc || R_ISDIR(fstat.fMode);
    9292
    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;
    9597    delete [] buf;
    9698
    … …  
    119121
    120122    // 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;
    130124
    131125    // ------------------------------------------------------
    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 range
    137     // of integration
    138     const int first_slice =  20; //  10ns
    139     const int last_slice  = 250; // 125ns
    140 
    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 time
    145     //    * fFallTime: Number of slices right from arrival time
    146     //  IntegralRel:
    147     //    * fRiseTime: Number of slices left  from maximum time
    148     //    * fFallTime: Number of slices right from maximum time
    149     //
    150     const int rise_time_cal = maximum ?  40 :  10; // was 13;   5ns
    151     const int fall_time_cal = maximum ? 120 : 160; // was 23;  80ns
    152 
    153     const int rise_time_dat = maximum ?  10 :   2; // was 13; was 10;   1ns
    154     const int fall_time_dat = maximum ?  40 :  48; // was 23; was 40;  24ns
    155 
    156     // Extraction type: Extract integral and half leading edge
    157 
    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;  // All
    165     Long_t max3 =  max;  // Pedestal Rndm
    166     Long_t max4 =  max;  // Pedestal Ext
    167 
    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     // ======================================================
    233126
    234127    MStatusDisplay *d = new MStatusDisplay;
    … …  
    236129    MBadPixelsCam badpixels;
    237130    badpixels.InitSize(1440);
    238     /*
    239131    badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    240132    badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    … …  
    243135    badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    244136    badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    245     */
    246137    //  Twin pixel
    247138    //     113
    … …  
    251142    //    1195
    252143    //    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    // ======================================================
    253231
    254232    MDrsCalibrationTime timecam;
    … …  
    268246
    269247    // ========================= Result ==================================
    270 
     248/*
    271249    //~ Double_t avgS = evt1f.GetHist()->GetMean();
    272250    //~ Double_t medS = evt1f.GetHist()->GetMedian();
    … …  
    340318    //hcalco.Scale(scale);
    341319    hcalco.DrawCopy();
    342 
     320*/
    343321    // ======================================================
    344322
    … …  
    372350    MDrsCalibApply drsapply3;
    373351
     352    MFilterData filterdata3;
     353
    374354    //---
    375355
    376     MExtractTimeAndChargeSpline extractor3;
     356    MExtractFACT extractor3;
    377357    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);
    383358    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 events
    390 //    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      */
    400359
    401360    MCalibrateFact conv3;
    402361    conv3.SetScale(scale);
    403     conv3.SetCalibConst(calib);
     362    //conv3.SetCalibConst(calib);
    404363
    405364    MBadPixelsTreat treat3;
    … …  
    420379    tlist3.AddToList(&drsapply3);
    421380    tlist3.AddToList(&cont3);
     381    tlist3.AddToList(&filterdata3);
    422382    tlist3.AddToList(&extractor3);
    423383//    tlist3.AddToList(&fill3a);
    … …  
    461421    MDrsCalibApply drsapply4;
    462422
    463     MExtractTimeAndChargeSpline extractor4;
     423    MFilterData filterdata4;
     424
     425    MExtractFACT extractor4;
    464426    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);
    470427    extractor4.SetNoiseCalculation(kFALSE);
    471428
    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 events
    477 //    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 */
    486429    MCalibrateFact conv4;
    487430    conv4.SetScale(scale);
    488     conv4.SetCalibConst(calib);
     431    //conv4.SetCalibConst(calib);
    489432
    490433    MBadPixelsTreat treat4;
    … …  
    504447    tlist4.AddToList(&drsapply4);
    505448    tlist4.AddToList(&cont4);
     449    tlist4.AddToList(&filterdata4);
    506450    tlist4.AddToList(&extractor4);
    507 //    tlist4.AddToList(&fill4a);
    508451    tlist4.AddToList(&conv4);
    509452    tlist4.AddToList(&treat4);
    … …  
    550493    MDrsCalibApply drsapply5;
    551494
     495    MTreatSaturation treatsat5;
     496
     497    MFilterData filterdata5;
     498
    552499    MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0",     "SelectDat");
    553500    MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
    … …  
    559506    // ---
    560507
    561     MExtractTimeAndChargeSpline extractor5dat;
     508    MExtractFACT extractor5dat;
    562509    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);
    568510    extractor5dat.SetNoiseCalculation(kFALSE);
    569511
    570     MExtractTimeAndChargeSpline extractor5cal;
     512    MExtractFACT extractor5cal;
    571513    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);
    577514    extractor5cal.SetNoiseCalculation(kFALSE);
    578515
    579     MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
     516    MExtractFACT extractor5tm("ExtractTM");
    580517    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);
    586518    extractor5tm.SetNoiseCalculation(kFALSE);
    587519    extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
    … …  
    593525
    594526    // ---
    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 */
    603527    MCalibrateFact conv5;
    604528    conv5.SetScale(scale);
    605     conv5.SetCalibConst(calib);
     529    //conv5.SetCalibConst(calib);
    606530
    607531    MCalibrateDrsTimes calctm5;
    … …  
    689613    tlist5.AddToList(&filterped);
    690614    tlist5.AddToList(&fill5a);
     615    tlist5.AddToList(&treatsat5);
     616    tlist5.AddToList(&filterdata5);
    691617    tlist5.AddToList(&extractor5dat);
    692618    tlist5.AddToList(&extractor5cal);
Note: See TracChangeset for help on using the changeset viewer.