Changeset 19951 for trunk


Ignore:
Timestamp:
04/29/20 21:48:05 (5 years ago)
Author:
giangdo
Message:
script to synchronise the data from HAWC's Eye and HAWC
Location:
trunk/Mars/hawc
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/Mars/hawc/callisto.C

    r19825 r19951  
    2727
    2828 To run the macro from the command line (assuming you are in a directory
    29  Mars/build where you have built your Mars environment) ou can do
     29 Mars/build where you have built your Mars environment) you can do
    3030
    3131    root ../hawc/callisto.C\(\"20191002_018.fits.fz\",\"20191002_013.drs.fits\"\)
    … …  
    5454
    5555    // mapping file (found in Mars/hawc)
    56     const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
     56    const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
    5757
    5858    //const char *pulse_template = "template-pulse.root";
    … …  
    7676    // Calibration constant (for the moment a single constant to
    7777    // convert the extracted charge to photo-electrons)
    78     double scale = 1./22.553;//0.2;
     78    //double scale = 1./22.553;//0.2;
     79    double scale = 1./17.557;
    7980
    8081    // ------------------------------------------------------
    … …  
    317318    write5.AddContainer("MSimSourcePos",       "Events", kFALSE);
    318319    write5.AddContainer("MTime",               "Events", kFALSE);
     320    write5.AddContainer("MRawBoardsFACT",      "Events", kFALSE);
    319321    write5.AddContainer("MSrcPosCam",          "Events", kFALSE);
    320322

  • trunk/Mars/hawc/extract_pulse.C

    r19876 r19951  
    6060};
    6161
    62 //Storage class to kKeep a list of the extracted single
     62//Storage class to keep a list of the extracted single
    6363class MSingles : public MParContainer, public MCamEvent
    6464{
    … …  
    6868
    6969public:
    70     MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
     70    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(40)
    7171    {
    7272        fName = name ? name : "MSingles";
    73         fName = title ? title : "Storeage for vector of singles";
     73        fName = title ? title : "Storage for vector of singles";
    7474    }
    7575
    … …  
    109109
    110110public:
    111     MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
     111    MHBaseline() : fPedestal(0), fSkipStart(5), fSkipEnd(10)
    112112    {
    113113        fName = "MHBaseline";
    … …  
    185185        for (int x=0; x<64; x++)
    186186        {
     187            //exclude the pixels without cones
     188            if (x == 61)
     189               continue;
     190            if (x == 62)
     191               continue;
     192            if (x == 63)
     193               continue;
     194
    187195            // Get the corresponding slice from the histogram
    188196            TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1);
    … …  
    286294class MHSingles : public MH
    287295{
    288     TH2F fSignalAll;
     296    TH2F fSignals; // all signals, without normalization
     297    TH2F fOsicllation;
     298
     299    TH2F fSignalAll; // Histogram of all signals
    289300    TH2F fSignalOut;
    290     TH2F fSignal1;
    291     TH2F fSignal2;
     301    TH2F fSignal1; // Histogram of signals with 1pe
     302    TH2F fSignal2; // Histogram of signals with 2pe
    292303    TH2F fSignal3;
    293304    TH2F fSignal4;
    … …  
    320331
    321332        // Setup histograms
     333
     334        fSignals.SetName("Signals");
     335        fSignals.SetTitle("");
     336        fSignals.SetXTitle("Time [ns]");
     337        fSignals.SetYTitle("Amplitude [mV]");
     338        fSignals.GetXaxis()->CenterTitle();
     339        fSignals.GetXaxis()->SetTitleSize(0.045);
     340        fSignals.GetYaxis()->SetTitleSize(0.045);
     341        fSignals.SetDirectory(NULL);
     342        fSignals.SetStats(kFALSE);
     343
     344        fOsicllation.SetName("Oscillations");
     345        fOsicllation.SetTitle("");
     346        fOsicllation.SetXTitle("Oscillation's Amplitude [mV]");
     347        fOsicllation.SetYTitle("Pulse's Amplitude [mV]");
     348        fOsicllation.GetXaxis()->CenterTitle();
     349        fOsicllation.GetXaxis()->SetTitleSize(0.045);
     350        fOsicllation.GetYaxis()->SetTitleSize(0.045);
     351        fOsicllation.SetDirectory(NULL);
     352        fOsicllation.SetStats(kFALSE);
     353
    322354        fSignalAll.SetName("SignalAll");
    323355        fSignalAll.SetTitle("");
    … …  
    498530        //const Int_t w = fSingles->GetIntegrationWindow();
    499531
    500         MBinning binsx, binsy;
     532        MBinning binsx, binsy, bin_time, bin_pulse, binosc, binpulse;
    501533        binsx.SetEdges(1024,  -512, 512);
    502534        binsy.SetEdges( 800,  -20,   20);
     535
     536        bin_time.SetEdges(1024, -512, 512);
     537        bin_pulse.SetEdges(240, -60, 60); //bin width is about 1 dac count = 0.5mV
     538
     539        binosc.SetEdges(140, -65, 5);
     540        binpulse.SetEdges(140, -5, 65);
     541
     542        MH::SetBinning(fSignals, bin_time, bin_pulse);
     543        MH::SetBinning(fOsicllation, binosc, binpulse);
    503544
    504545        MH::SetBinning(fSignalAll, binsx, binsy);
    … …  
    579620                continue;
    580621
    581             if (idx!=0)
     622            if (idx!=43)
    582623                continue;
    583624
    584             const double gain  = hgain->GetBinContent(idx+1)/30;
     625            const double gain  = hgain->GetBinContent(idx+1)/40;
    585626            const double base  = hbase->GetBinContent(idx+1);
    586627
    … …  
    594635            // Fill pulse shape
    595636            const double tm  = result[0].fTime;
    596             const double sig = result[0].fSignal/30 - base;
     637            const double sig = result[0].fSignal/40 - base;
    597638
    598639            bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5;
    599640
    600             double mean = 0;
    601 
    602             for (int s=20; s<roi-20; s++)
     641            double corr = -0.25;
     642
     643            // Make histograms to check correlation between signals and oscillations
     644            for(int s=0; s<roi-10; s++)
     645            {   double time = (s-tm)/2;
     646                double pulse = ptr[s]-ped-base; // corrected pulse without normalization
     647
     648                // Fill all extracted signals into one histogram
     649                fSignals.Fill(time, pulse);
     650            }
     651
     652            int t_arrival = (int)tm;
     653            double max_pulse = ptr[t_arrival]-ped-base;
     654            int max_pulse_pos = 0; // position of the maximal pulse in sample
     655            // Find pulse's amplitude
     656            for(int s = tm; s < tm + 15; s++)
    603657            {
    604                 double x = (s - tm)/2; // [ns]
     658              if (ptr[s]-ped-base > max_pulse)
     659                 max_pulse = ptr[s]-ped-base;
     660                 max_pulse_pos = s;
     661            }
     662            //cout << "________Check if the filled histograms are correct________" << endl;
     663            //cout << "Arrival time: " << tm << ",Amplitude: " << ptr[t_arrival]-ped-base << endl;
     664            //cout << "Pulse's position: "<< max_pulse_pos << ",Max. Pulse's Amplitude: " << ptr[max_pulse_pos]-ped-base << endl;
     665
     666            // Find maximal oscillation's amplitude
     667            double max_osc = 0;
     668            for (int osc = max_pulse_pos; osc < max_pulse_pos+40; osc++)
     669            {
     670                if (ptr[osc]-ped-base < max_osc)
     671                {
     672                  max_osc = ptr[osc]-ped-base;
     673                }
     674                else
     675                  continue;
     676            }
     677            fOsicllation.Fill(max_osc, max_pulse);
     678
     679            // Set conditions for 1pe signals
     680            for (int s=5; s<roi-10; s++)
     681            {
     682                double x = (s - tm)/2; // return the time in [ns]
    605683                double y = (ptr[s]-ped-base)/gain;
    606684
    607                 if (x>-100 && x<-10)
    608                     mean += y/90;
    609 
    610                 if (y<-3)
     685                //if (x>-100 && x<-10)
     686                //    mean += y/90;
     687
     688                if (x>5 && y>2.5+x*-0.01)
     689                    ok = false;
     690                if (y<-1.5)
    611691                    ok = false;
    612692                if (x>5 && x<40 && y<=x/40*-0.5)
    … …  
    614694                if (x>-150 && x<-5 && y>1.5)
    615695                    ok = false;
    616                 if (y>2.5)
     696                if (x<0 && y > 1.2)
    617697                    ok = false;
     698
    618699            }
    619700
    620             // Nothing at or bleow zero from 5 to 40machen
    621 
     701            // Fill nPE signals into histograms
     702
     703            // Set conditions for signals with n PE
     704            // Nothing at or below zero from 5 to 40machen
    622705            bool ok1 = sig>gain*0.5 && sig<gain*1.5;
    623706            bool ok2 = sig>gain*1.5 && sig<gain*2.5;
    … …  
    630713
    631714
    632             for (int s=20; s<roi-20; s++)
     715            for (int s=5; s<roi-10; s++)
    633716            {
    634717                double x = (s - tm)/2; // [ns]
    635                 double y = (ptr[s]-ped-base)/gain -  mean;
     718                double y = (ptr[s]-ped-base)/gain - corr;
    636719
    637720                int ix = TMath::Nint((x-x0)/wx*nx);
    … …  
    642725
    643726                int bin = (iy+1)*(nx+2)+(ix+1);
     727
    644728
    645729                if (ok)
    … …  
    724808
    725809    // Getter for histograms
     810    TH2 *GetSignals() { return &fSignals; }
     811    TH2 *GetOscillation() { return &fOsicllation; }
    726812    TH2 *GetSignalAll() { return &fSignalAll; }
    727813    TH2 *GetSignalOut() { return &fSignalOut; }
    … …  
    774860        gPad->SetFrameBorderMode(0);
    775861        fSignal8.Draw("colz");
    776         //fProf8.Draw("same");
     862        fProf8.Draw("same");
    777863    }
    778864
    … …  
    836922    MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
    837923        fExtractionRange(64),
    838         fNumAverage(10), fStartSkip(20), fEndSkip(10),
    839         fIntegrationSize(3*10), fMaxSearchWindow(30)
     924        fNumAverage(10), fStartSkip(10), fEndSkip(10),
     925        fIntegrationSize(4*10), fMaxSearchWindow(20)
    840926    {
    841927    }
    … …  
    877963
    878964        const size_t navg              = fNumAverage;
    879         const float  threshold         = fThreshold; 
     965        const float  threshold         = fThreshold;
    880966        const UInt_t start_skip        = fStartSkip;
    881967        const UInt_t end_skip          = fEndSkip;
    … …  
    9341020                if (k_max == i+5 || k_max == i + max_search_window-1)
    9351021                    continue;
     1022
     1023                //Make sure the pulse is isolated (area after the maximum is empty)
     1024                UInt_t k_falling = k_max+200;
     1025                for (UInt_t k_after = k_falling; k_after < k_max +250; k_after++)
     1026                {
     1027                    if (val[k_after] > val[k_falling] &&
     1028                        val[k_after + fNumAverage/2] > val[k_falling])
     1029                    continue;
     1030                }
    9361031
    9371032                //search for half maximum before maximum
    … …  
    9461041                    }
    9471042                }
     1043
    9481044                // Check if the finding makes sense
    9491045                if (k_half_max+integration_size > roi-navg-end_skip)
    … …  
    9541050                    continue;
    9551051
    956                 // FIXME: Evaluate arrival time more precisely!!!
    957                 // FIXME: Get a better integration window
     1052                //cout << "________Check if the finding makes sense________" << endl;
     1053                //cout << "Pulse's position: " << k_max << ", Max. Pulse's Amplitude: " << val[k_max] << endl;
     1054                //cout << "Half maximum: " << k_half_max << ", Pulse at half maximum: " << val[k_half_max] << endl;
     1055
    9581056
    9591057                // Compile "single"
    … …  
    9701068                // Add single to list
    9711069                result.push_back(single);
     1070                break; //accept only one pulse per event
    9721071
    9731072                // skip falling edge
    … …  
    9851084                    Int_t        max_dist   = 14,
    9861085                    Double_t     threshold  =  5,
    987                     const char  *runfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_",
     1086                    const char  *runfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_",
    9881087                    int          firstRunID =  6,
    9891088                    int          lastRunID  =  6,
    990                     const char  *drsfile    = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_003.drs.fits",
     1089                    const char  *drsfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_004.drs.fits",
    9911090                    const char  *outpath    = "."
    9921091                   )
    9931092{
    9941093    // ======================================================
    995     if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
    996     {
    997         gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl;
     1094    if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt"))
     1095    {
     1096        gLog << err << "HAWCsEyemap181214.txt not found." << endl;
    9981097        return 1;
    9991098    }
    … …  
    10051104
    10061105    // built output filename and path
    1007     TString outfile = "20191027_006_006-pulse.root";
    1008     TString infile  = "20190227_006_006-fit.root";
     1106    TString outfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_pulse.root";
     1107    TString infile  = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root";
    10091108
    10101109    //if (!gSystem->AccessPathName(outfile))
    … …  
    10481147
    10491148    // map file to use (get that from La Palma!)
    1050     const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
     1149    const char *map = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
    10511150
    10521151    // ------------------------------------------------------
    10531152
    1054     Long_t max0 = 0;
    1055     Long_t max1 = 0;
     1153    Long_t max0 = 10000;
     1154    Long_t max1 = 10000;
    10561155
    10571156    // ======================================================
    … …  
    10971196    MBadPixelsCam badpixels;
    10981197    badpixels.InitSize(64);
    1099     badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    1100     badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
     1198    badpixels[43].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    11011199
    11021200    MH::SetPalette();
    … …  
    12191317    gStyle->SetOptStat(10);
    12201318
     1319    TH2 *hsignals = h->GetSignals();
     1320    TH2 *hoscillation = h->GetOscillation();
     1321
     1322    TH2 *hall = h->GetSignalAll();
     1323
    12211324    TH1 *h1 = h->GetSignal1();
    12221325    TH1 *p1 = h->GetProf1();
    12231326
     1327    cout << "Entries:\t" << (p1->GetEntries()) << endl;
     1328
    12241329    TH1 *h2 = h->GetSignal2();
    12251330    TH1 *p2 = h->GetProf2();
    … …  
    12431348    TH1 *p8 = h->GetProf8();
    12441349
    1245     TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
    1246     pulse.SetLineColor(kBlack);
    1247     pulse.SetParameters(0, 1./0.094, 1./0.05);
     1350    // This is the old fit
     1351    //TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300);
     1352    //pulse.SetLineColor(kBlack);
     1353    //pulse.SetParameters(0, 1./0.094, 1./0.05);
     1354
     1355    TF1 pulse("pulse", "[0]*(1-1/(1+exp((x-[1])/[2])))*exp(-(x-[1])/[3])", 0, 250);
     1356    pulse.SetLineColor(kRed);
     1357    pulse.SetParLimits(1, 1.0, 3);
     1358    pulse.SetParLimits(2, 0.5, 2.0);
     1359    pulse.SetParLimits(3, 50, 110);
    12481360
    12491361    // =====
    … …  
    12591371
    12601372    // =====
    1261 
     1373    d->AddTab("Extracted Signals");
     1374    gPad->SetGrid();
     1375    gPad->SetLeftMargin(0.1);
     1376    gPad->SetTopMargin(0.1);
     1377    gPad->SetRightMargin(0.15);
     1378    hsignals->DrawCopy("colz");
     1379
     1380    d->AddTab("OscillationAmpl");
     1381    gPad->SetGrid();
     1382    gPad->SetLeftMargin(0.1);
     1383    gPad->SetTopMargin(0.1);
     1384    gPad->SetRightMargin(0.15);
     1385    hoscillation->DrawCopy("colz");
     1386
     1387
     1388    d->AddTab("All Signals");
     1389    gPad->SetGrid();
     1390    gPad->SetLeftMargin(0.1);
     1391    gPad->SetTopMargin(0.1);
     1392    gPad->SetRightMargin(0.15);
     1393    hall->DrawCopy("colz");
     1394
     1395
     1396    // =====
    12621397    d->AddTab("1");
    12631398    gPad->SetGrid();
    1264     gPad->SetLeftMargin(0.07);
    1265     gPad->SetTopMargin(0.01);
     1399    gPad->SetLeftMargin(0.1);
     1400    gPad->SetTopMargin(0.1);
    12661401    gPad->SetRightMargin(0.15);
    12671402    h1->DrawCopy("colz");
    12681403    p1->DrawCopy("same");
    12691404
    1270     pulse.SetParameter(0, 1*1.3);
    1271     p1->Fit(&pulse, "N0", "", 1.5, 25);
    1272     TF1 *f1 = pulse.DrawCopy("same");
     1405    //pulse.FixParameter(0, 1*1.5);
     1406    //p1->Fit(&pulse, "N0", "", 1, 70);
     1407    //TF1 *f1 = pulse.DrawCopy("same");
    12731408
    12741409    // =====
    … …  
    12761411    d->AddTab("2");
    12771412    gPad->SetGrid();
    1278     gPad->SetLeftMargin(0.07);
    1279     gPad->SetTopMargin(0.01);
     1413    gPad->SetLeftMargin(0.1);
     1414    gPad->SetTopMargin(0.1);
    12801415    gPad->SetRightMargin(0.15);
    12811416    h2->DrawCopy("colz");
    … …  
    12901425    d->AddTab("3");
    12911426    gPad->SetGrid();
    1292     gPad->SetLeftMargin(0.07);
    1293     gPad->SetTopMargin(0.01);
     1427    gPad->SetLeftMargin(0.1);
     1428    gPad->SetTopMargin(0.1);
    12941429    gPad->SetRightMargin(0.15);
    12951430    h3->DrawCopy("colz");
    … …  
    13041439    d->AddTab("4");
    13051440    gPad->SetGrid();
    1306     gPad->SetLeftMargin(0.07);
    1307     gPad->SetTopMargin(0.01);
     1441    gPad->SetLeftMargin(0.1);
     1442    gPad->SetTopMargin(0.1);
    13081443    gPad->SetRightMargin(0.15);
    13091444    h4->DrawCopy("colz");
    … …  
    13181453    d->AddTab("5");
    13191454    gPad->SetGrid();
    1320     gPad->SetLeftMargin(0.07);
    1321     gPad->SetTopMargin(0.01);
     1455    gPad->SetLeftMargin(0.1);
     1456    gPad->SetTopMargin(0.1);
    13221457    gPad->SetRightMargin(0.15);
    13231458    h5->DrawCopy("colz");
    … …  
    13321467    d->AddTab("6");
    13331468    gPad->SetGrid();
    1334     gPad->SetLeftMargin(0.07);
    1335     gPad->SetTopMargin(0.01);
     1469    gPad->SetLeftMargin(0.1);
     1470    gPad->SetTopMargin(0.1);
    13361471    gPad->SetRightMargin(0.15);
    13371472    h6->DrawCopy("colz");
    … …  
    13461481    d->AddTab("7");
    13471482    gPad->SetGrid();
    1348     gPad->SetLeftMargin(0.07);
    1349     gPad->SetTopMargin(0.01);
     1483    gPad->SetLeftMargin(0.1);
     1484    gPad->SetTopMargin(0.1);
    13501485    gPad->SetRightMargin(0.15);
    13511486    h7->DrawCopy("colz");
    … …  
    13601495    d->AddTab("8");
    13611496    gPad->SetGrid();
    1362     gPad->SetLeftMargin(0.07);
    1363     gPad->SetTopMargin(0.01);
     1497    gPad->SetLeftMargin(0.1);
     1498    gPad->SetTopMargin(0.1);
    13641499    gPad->SetRightMargin(0.15);
    13651500    h8->DrawCopy("colz");
    … …  
    13741509    d->AddTab("Pulse");
    13751510    gPad->SetGrid();
    1376     gPad->SetLeftMargin(0.07);
    1377     gPad->SetTopMargin(0.01);
     1511    gPad->SetLeftMargin(0.1);
     1512    gPad->SetTopMargin(0.1);
    13781513    gPad->SetRightMargin(0.15);
    13791514    p8->SetStats(kFALSE);
    … …  
    13941529    p1->DrawCopy("same");
    13951530
    1396     f1->Draw("same");
     1531    //f1->Draw("same");
    13971532    f2->Draw("same");
    13981533    f3->Draw("same");
    … …  
    14061541
    14071542
     1543
     1544    d->SaveAs(outfile);
     1545
     1546    TFile f(outfile, "UPDATE");
    14081547    /*
    1409     d->SaveAs(outfile);
    1410 
    1411     TFile f(outfile, "UPDATE");
    1412 
    14131548    MParameterI par("NumEvents");
    14141549    par.SetVal(fill.GetNumExecutions());

  • trunk/Mars/hawc/extract_singles.C

    r19857 r19951  
    5757{
    5858    Int_t fIntegrationWindow;
    59 
    6059    vector<vector<Single>> fData;
    6160
    6261public:
    63     MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30)
     62    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(40)
    6463    {
    6564        fName = name ? name : "MSingles";
    … …  
    8887    ClassDef(MSingles, 1)
    8988};
    90 
    9189// Histogram class to extract the baseline
    9290class MHBaseline : public MH
    … …  
    102100
    103101public:
    104     MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10)
     102    MHBaseline() : fPedestal(0), fSkipStart(10), fSkipEnd(10)
    105103    {
    106104        fName = "MHBaseline";
    … …  
    108106        // Setup the histogram
    109107        fBaseline.SetName("Baseline");
    110         fBaseline.SetTitle("Median spectrum");
     108        fBaseline.SetTitle("Baseline Spectrum");
    111109        fBaseline.SetXTitle("Pixel [idx]");
    112         fBaseline.SetYTitle("Median baseline [mV]");
     110        fBaseline.SetYTitle("Amplitude [mV]");
    113111        fBaseline.SetDirectory(NULL);
    114112    }
    … …  
    282280    TH2F       fTime;
    283281    TProfile2D fPulse;
     282    TH2F       fPulse2D; //plot the extracted singles in a 2D plot vs samples(time)
    284283
    285284    UInt_t fNumEntries;
    … …  
    296295        // Setup histograms
    297296        fSignal.SetName("Signal");
    298         fSignal.SetTitle("pulse integral spectrum");
     297        fSignal.SetTitle("Pulse Integral Spectrum");
    299298        fSignal.SetXTitle("pixel [SoftID]");
    300299        fSignal.SetYTitle("time [a.u.]");
    … …  
    302301
    303302        fTime.SetName("Time");
    304         fTime.SetTitle("arival time spectrum");
     303        fTime.SetTitle("Arrival Time Spectrum");
    305304        fTime.SetXTitle("pixel [SoftID]");
    306305        fTime.SetYTitle("time [a.u.]");
    … …  
    308307
    309308        fPulse.SetName("Pulse");
    310         fPulse.SetTitle("average pulse shape spectrum");
     309        fPulse.SetTitle("Average Pulse Shape Spectrum");
    311310        fPulse.SetXTitle("pixel [SoftID]");
    312311        fPulse.SetYTitle("time [a.u.]");
    313312        fPulse.SetDirectory(NULL);
     313
     314        fPulse2D.SetName("Pulse2D");
     315        fPulse2D.SetTitle("Extracted Pulses");
     316        fPulse2D.SetXTitle("time [a.u.]");
     317        fPulse2D.SetYTitle("amplitude [a.u.]");
     318        fPulse2D.SetDirectory(NULL);
    314319    }
    315320
    … …  
    363368
    364369        // Correct for DRS timing!!
     370
     371        //Setup binning for the average time spectrum
    365372        MBinning binst(roi, -0.5, roi-0.5);
    366373        MH::SetBinning(fTime, binsx, binst);
    367374
     375        //Setup binning for the average pulse spectrum
    368376        MBinning binspy(2*roi, -roi-0.5, roi-0.5);
     377        //MBinning binspy(1024, -522-0.5, 522-0.5);
    369378        MH::SetBinning(fPulse, binsx, binspy);
     379
     380        //Setup binning for the 2D plot amplitude vs time
     381        MBinning bins_time(2*roi, -roi-0.5, roi-0.5);
     382        //MBinning bins_time(1024, -522-0.5, 522-0.5);
     383        MBinning bins_ampl(160, -30, 50); //N_bin is chosen such that 1 DAC = 0.5mV
     384        MH::SetBinning(fPulse2D, bins_time, bins_ampl);
    370385
    371386        return kTRUE;
    … …  
    402417                for (int s=0; s<roi; s++)
    403418                    fPulse.Fill(i, s-tm, ptr[s]);
     419
     420                for (int s=0; s<roi; s++)
     421                    fPulse2D.Fill(s-tm, ptr[s]);
    404422            }
    405423
    … …  
    416434    const TH2 *GetTime() const   { return &fTime; }
    417435    const TH2 *GetPulse() const  { return &fPulse; }
     436    const TH2F *GetPulse2D() const {return &fPulse2D;}
    418437
    419438    UInt_t GetNumEntries() const { return fNumEntries; }
    … …  
    441460        gPad->SetFrameBorderMode(0);
    442461        fPulse.Draw("colz");
     462
     463        pad->cd(4);
     464        gPad->SetBorderMode(0);
     465        gPad->SetFrameBorderMode(0);
     466        fPulse2D.Draw("colz");
    443467    }
    444468
    … …  
    465489        gPad->SetFrameBorderMode(0);
    466490        fPulse.DrawCopy("colz");
     491
     492        pad->cd(4);
     493        gPad->SetBorderMode(0);
     494        gPad->SetFrameBorderMode(0);
     495        fPulse2D.Draw("colz");
    467496    }
    468497
    … …  
    497526    MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0),
    498527        fExtractionRange(64),
    499         fNumAverage(10), fStartSkip(20), fEndSkip(10),
    500         fIntegrationSize(3*10), fMaxSearchWindow(30)
     528        fNumAverage(10), fStartSkip(10), fEndSkip(10),
     529        fIntegrationSize(4*10), fMaxSearchWindow(20)
    501530    {
    502531    }
    … …  
    538567
    539568        const size_t navg              = fNumAverage;
    540         const float  threshold         = fThreshold; 
     569        const float  threshold         = fThreshold;
    541570        const UInt_t start_skip        = fStartSkip;
    542571        const UInt_t end_skip          = fEndSkip;
    … …  
    595624                if (k_max == i+5 || k_max == i + max_search_window-1)
    596625                    continue;
     626
     627                //Make sure the pulse is isolated (area after the maximum is empty)
     628                UInt_t k_falling = k_max+200;
     629                for (UInt_t k_after = k_falling; k_after < k_max +250; k_after++)
     630                {
     631                    if (val[k_after] > val[k_falling] &&
     632                        val[k_after + fNumAverage/2] > val[k_falling])
     633                    continue;
     634                }
    597635
    598636                //search for half maximum before maximum
    … …  
    623661                single.fTime   = k_half_max;
    624662
    625                 // Crossing of the threshold is at 5
     663                // Crossing of the threshold is at 7
    626664                for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
    627665                    single.fSignal += ptr[j];
    … …  
    631669                // Add single to list
    632670                result.push_back(single);
     671                break; //accept only one pulse per event
    633672
    634673                // skip falling edge
    … …  
    644683
    645684int extract_singles(
    646                     Int_t        max_dist   = 14,
    647                     Double_t     threshold  = 5,
    648                     const char  *runfile    = "raw/2019/10/27/20191027_",
    649                     int          firstRunID =  6,
    650                     int          lastRunID  =  6,
    651                     const char  *drsfile    = "raw/2019/10/27/20191027_003.drs.fits",
    652                     const char  *outpath    = "."
     685                      Int_t        max_dist   = 14,
     686                      Double_t     threshold  = 5,
     687                      const char  *runfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_",
     688                      int          firstRunID =  6,
     689                      int          lastRunID  =  6,
     690                      const char  *drsfile    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_004.drs.fits",
     691                      const char  *outpath    = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted"
    653692                   )
    654693{
    … …  
    688727
    689728    // map file to use (get that from La Palma!)
    690     const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
     729    const char *map = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
    691730
    692731    // ------------------------------------------------------
    … …  
    856895    const TH2 *htime   = h->GetTime();
    857896    const TH2 *hpulse  = h->GetPulse();
     897    const TH2F *ppulse2d = h->GetPulse2D();
    858898
    859899    d->AddTab("Time");
    … …  
    869909    hp->Draw();
    870910
     911    d->AddTab("2DPulse");
     912    ppulse2d->DrawCopy("colz");
     913
    871914    d->SaveAs(outfile);
    872915

  • trunk/Mars/hawc/fit_singles.C

    r19875 r19951  
    9393
    9494// Print function parameters
    95 void PrintFunc(TF1 &f, float integration_window=30)
     95void PrintFunc(TF1 &f, float integration_window=40)
    9696{
    9797    //cout << "--------------------" << endl;
    … …  
    113113// The parameters for the function are the filenam from the gain extraction
    114114// and the output filename
    115 int fit_singles(const char *filename = "20191027_006_006.root",
    116                 const char *outfile  = "20190227_006_006-fit.root",
     115int fit_singles(const char *filename = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006.root",
     116                const char *outfile  = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root",
    117117                bool fixednoise=false)
    118118{
    119119    // Read the mapping between bias channels and hardware channels
    120120    PixelMap pmap;
    121     if (!pmap.Read("../hawc/FAMOUSmap171218.txt"))
    122     {
    123         cout << "FAMOUSmap171218.txt not found." << endl;
     121    if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt"))
     122    {
     123        cout << "HAWCsEyemap181214.txt not found." << endl;
    124124        return 1;
    125125    }
    126126
    127127    // It is more correct to fit the integral, but this is super
    128     // slow, especially for 1440 pixel. To get that, one would have
     128    // slow. To get that, one would have
    129129    // to analytically integrate and fit this function.
    130     // (Currently the integral is switched off for the 1440 individual
    131     // spectra in all cases)
    132130    bool fast = false; // Switch off using integral
    133131
    134132    // Values which should be read from the file but not available atm
    135     Int_t integration_window = 30;
     133    Int_t integration_window = 40;
    136134
    137135    // Map for which pixel shall be plotted and which not
    … …  
    266264    TH1F hCoeffR2    ("CoeffR2",    "Coefficient R",          90, -1,     2);
    267265
    268     // Histigram for the sum of all spectrums
    269     TH1F hSum("Sum1", "Signal spectrum of all pixels",
     266    // Histogram for the sum of all spectrums
     267    TH1F hSum("Sum1", "Signal spectrum of all pixels (excl. blind pixels)",
    270268              hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
    271269    hSum.SetXTitle("pulse integral [mV sample]");
    … …  
    275273
    276274    // Histogram for the sum of all pixels excluding the ones with faulty fits
    277     TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)",
     275    TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (excl. faulty and blind pixels)",
    278276                    hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
    279277    hSumClear1.SetXTitle("pulse integral [mV sample]");
    … …  
    303301    hPulse.SetStats(false);
    304302
    305     // Spektrum for the normalized individual spectra
    306     TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)",  120, 0.05, 12.05);
     303    // Spectrum for the normalized individual spectra
     304    TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (excl. blind pixels)",  120, 0.05, 12.05);
    307305    hSumScale1.SetXTitle("pulse integral [pe]");
    308306    hSumScale1.SetYTitle("Rate");
    … …  
    332330        if (usePixel[pixel]==0)
    333331            continue;
     332
     333        //exclude the pixels without cones
     334        if (pixel == 61)
     335           continue;
     336        if (pixel == 62)
     337           continue;
     338        if (pixel == 63)
     339           continue;
     340
     341
     342
    334343
    335344        TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
    … …  
    434443        if (usePixel[pixel]==0)
    435444            continue;
     445
     446        // Exclude pixels without cones
     447        if (pixel == 61)
     448           continue;
     449        if (pixel == 62)
     450           continue;
     451        if (pixel == 63)
     452           continue;
    436453
    437454        //Projectipon of a certain Pixel out of Ampl.Spectrum

  • trunk/Mars/hawc/plot_trace.C

    r19846 r19951  
     1#include <TStyle.h>
     2#include <TCanvas.h>
     3#include <TSystem.h>
     4#include <TF1.h>
     5#include <TProfile.h>
     6#include <TProfile2D.h>
     7#include <TMath.h>
     8#include <TGraph.h>
     9#include <TLine.h>
     10#include <TFitResultPtr.h>
     11#include <TFitResult.h>
     12#include <TFile.h>
     13#include <TLine.h>
     14
     15#include <cstdio>
     16#include <stdio.h>
     17#include <stdint.h>
     18
     19#include "Getline.h"
     20#include "MH.h"
     21#include "MArrayI.h"
     22#include "MDirIter.h"
     23#include "MFillH.h"
     24#include "MEvtLoop.h"
     25#include "MCamEvent.h"
     26#include "MGeomApply.h"
     27#include "MTaskList.h"
     28#include "MParList.h"
     29#include "MContinue.h"
     30#include "MBinning.h"
     31#include "MDrsCalibApply.h"
     32#include "MDrsCalibration.h"
     33#include "MRawFitsRead.h"
     34#include "MStatusDisplay.h"
     35#include "MTaskInteractive.h"
     36#include "MPedestalSubtractedEvt.h"
     37#include "MGeomCamFAMOUS.h"
     38#include "MRawRunHeader.h"
     39#include "MPedestalCam.h"
     40#include "MPedestalPix.h"
     41#include "MParameters.h"
     42
    143// ==========================================================================
    244// ============ see plot_callisto function at the end of the file ===========
    … …  
    2163
    2264// Create a histogram for the trace
    23 TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75);
     65//TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75); //range = [-0.25,511.75]ns
     66TH1D hist("Trace", "Waveform (Pixel 5, Event 1)", 1024, -0.5, 1024-0.5);
    2467
    2568// All 'data members' that are required globally
    … …  
    3477{
    3578    fEvt = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
    36     if (!fEvt)
    37     {
    38         gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
    39         return kALSE;
    40     }
     79
     80    //if (!fEvt)
     81    //{
     82    //    gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
     83    //    return kFALSE;
     84    //}
    4185
    4286    // Create a canvas and plot the histogram into it
    43     c = new TCanvas;
     87    c = new TCanvas();
    4488    hist.SetStats(kFALSE);
    4589    hist.Draw();
    46 
    47     TF1 constant("zero", "[0]", -0.25, 511.75);
     90    hist.SetXTitle("Samples");
     91    hist.SetYTitle("Amplitude [a.u.]");
     92
     93    //TF1 constant("zero", "[0]", -0.25, 511.75);
     94    TF1 constant("zero", "[0]", -0.5, 1024);
    4895    constant.SetParameter(0, 0);
    4996    constant.SetLineColor(kBlack);
    … …  
    59106    constant.DrawCopy("same");
    60107
     108
    61109    return kTRUE;
    62110}
    … …  
    64112Int_t Process()
    65113{
    66     int pixel_index = 0;
     114    int pixel_index = 6; //plot trace only of this pixel
     115    int skipstart = 5;
     116    int skipend = 10;
    67117
    68118    // Number of samples per pixel (usually 1024)
    … …  
    73123    hist.Reset();
    74124    for (int i=0; i<numsamples; i++)
    75         hist.Fill(i*0.5, cal[i]); // Convert sample to nano-second
     125        hist.Fill(i, cal[i]); // 0.5ns = 1 sample
    76126
    77127    // Reset min/max range
    78     hist.SetMinimum();
    79     hist.SetMaximum();
    80 
    81     // Set a reasonable range
    82     if (hist.GetMinimum()>-50)
    83         hist.SetMinimum(-50);
    84     if (hist.GetMaximum()<250)
    85         hist.SetMaximum(250);
     128    hist.SetMinimum(-50);
     129    hist.SetMaximum(50);
     130
     131    //// Set a reasonable range
     132    //if (hist.GetMinimum()>-30)
     133    //    hist.SetMinimum(-20);
     134    //if (hist.GetMaximum()<60)
     135    //    hist.SetMaximum(50);
     136
     137
     138
     139    TLine *line = new TLine(skipstart, 50, skipstart, -50);
     140    line->SetLineStyle(kSolid);
     141    line->SetLineColor(kRed);
     142    line->Draw();
     143
     144    TLine *line2 = new TLine(1024-skipend, 50, 1024-skipend, -50);
     145    line2->SetLineStyle(kSolid);
     146    line2->SetLineColor(kRed);
     147    line2->Draw();
    86148
    87149    // Signal root to update the canvas/pad
    88150    c->Modified();
    89151    c->Update();
     152    //c->SaveAs("/home/giangdo/Documents/Master/MA/pulse/data_single_pe/traces/pixel3/event_0.png");
     153
    90154
    91155    bool rc = HandleInput();
    … …  
    131195
    132196    // mapping file (found in Mars/hawc)
    133     const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL;
     197    const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;
    134198
    135199    // ======================================================
    … …  
    193257    MEvtLoop loop(gSystem->BaseName(datafile));
    194258    loop.SetParList(&plist);
    195  
    196     if (!loop.Eventloop())
     259
     260    if (!loop.Eventloop(100))
    197261        return 4;
    198262

  • trunk/Mars/hawc/star.C

    r19821 r19951  
    1818
    1919 To run the macro from the command line (assuming you are in a directory
    20  Mars/build where you have built your Mars environment) ou can do
     20 Mars/build where you have built your Mars environment) you can do
    2121
    2222    root ../hawc/star.C\(\"00000015.003_Y_MonteCarlo003_Events.root\"\)
    … …  
    206206    write.AddContainer("MPointingPos",        "Events", kFALSE);
    207207    write.AddContainer("MTime",               "Events", kFALSE);
     208    write.AddContainer("MRawBoardsFACT",      "Events", kFALSE);
    208209    write.AddContainer("MSrcPosCam",          "Events", kFALSE);
    209210
Note: See TracChangeset for help on using the changeset viewer.