Index: /trunk/Mars/fact/analysis/gain/extract_singles.C =================================================================== --- /trunk/Mars/fact/analysis/gain/extract_singles.C (revision 17033) +++ /trunk/Mars/fact/analysis/gain/extract_singles.C (revision 17033) @@ -0,0 +1,892 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "MH.h" +#include "MArrayI.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MDirIter.h" +#include "MFillH.h" +#include "MEvtLoop.h" +#include "MCamEvent.h" +#include "MHCamEvent.h" +#include "MGeomApply.h" +#include "MTaskList.h" +#include "MParList.h" +#include "MContinue.h" +#include "MBinning.h" +#include "MDrsCalibApply.h" +#include "MDrsCalibration.h" +#include "MRawFitsRead.h" +#include "MBadPixelsCam.h" +#include "MStatusDisplay.h" +#include "MTaskInteractive.h" +#include "MPedestalSubtractedEvt.h" +#include "MHCamera.h" +#include "MGeomCamFACT.h" +#include "MRawRunHeader.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MParameters.h" + +using namespace std; +using namespace TMath; + +// Structure to store the extracted value and the extracted time +struct Single +{ + float fSignal; + float fTime; +}; + +//Storage class to kKeep a list of the extracted single +class MSingles : public MParContainer, public MCamEvent +{ + Int_t fIntegrationWindow; + + vector> fData; + +public: + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30) + { + fName = name ? name : "MSingles"; + fName = title ? title : "Storeage for vector of singles"; + } + + void InitSize(const UInt_t i) + { + fData.resize(i); + } + + vector &operator[](UInt_t i) { return fData[i]; } + vector &GetVector(UInt_t i) { return fData[i]; } + + UInt_t GetNumPixels() const { return fData.size(); } + + void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; } + Int_t GetIntegrationWindow() const { return fIntegrationWindow; } + + Bool_t GetPixelContent(Double_t &, Int_t , const MGeomCam &, Int_t) const + { + return kTRUE; + } + void DrawPixelContent(Int_t) const { } + + ClassDef(MSingles, 1) +}; + +// Histogram class to extract the baseline +class MHBaseline : public MH +{ + TH2F fBaseline; + + MPedestalCam *fPedestal; + + // The baseline is only extracted where also the signal is extracted + // FIXME: Make sure this is consistent with MExtractSingles + UShort_t fSkipStart; + UShort_t fSkipEnd; + +public: + MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10) + { + fName = "MHBaseline"; + + // Setup the histogram + fBaseline.SetName("Baseline"); + fBaseline.SetTitle("Median spectrum"); + fBaseline.SetXTitle("Pixel [idx]"); + fBaseline.SetYTitle("Median baseline [mV]"); + fBaseline.SetDirectory(NULL); + } + + Bool_t ReInit(MParList *plist) + { + fPedestal = (MPedestalCam*)plist->FindCreateObj("MPedestalCam"); + if (!fPedestal) + return kFALSE; + + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << err << "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + // Bin width should be around 1 dac count which is about 0.5mV + MBinning binsx, binsy; + binsx.SetEdges(header->GetNumNormalPixels(), -0.5, header->GetNumNormalPixels()-0.5); + binsy.SetEdges(100, -20.5, 29.5); + + // Setup binnning + MH::SetBinning(fBaseline, binsx, binsy); + + return kTRUE; + } + + // Fill the samples into the histogram + Int_t Fill(const MParContainer *par, const Stat_t) + { + const MPedestalSubtractedEvt *evt = dynamic_cast(par); + + const Int_t n = evt->GetNumSamples()-fSkipStart-fSkipEnd; + + // Loop over all pixels + for (int pix=0; pixGetNumPixels(); pix++) + { + // Get samples for each pixel + const Float_t *ptr = evt->GetSamples(pix); + + // Average two consecutive samples + for (int i=0; iInitSize(fBaseline.GetNbinsX()); + fPedestal->SetNumEvents(GetNumExecutions()); + + Int_t cnt = 0; + Double_t avg = 0; + Double_t rms = 0; + + // Loop over all 'pixels' + for (int x=0; xGetMaximumBin(); + + // Calculate a parabola through this and the surrounding points + // on logarithmic values (that's a gaussian) + + // + // Quadratic interpolation + // + // calculate the parameters of a parabula such that + // y(i) = a + b*x(i) + c*x(i)^2 + // y'(i) = b + 2*c*x(i) + // + // + + // ------------------------------------------------------------------------- + // a = y2; + // b = (y3-y1)/2; + // c = (y3+y1)/2 - y2; + + const Double_t v1 = hist->GetBinContent(bin-1); + const Double_t v2 = hist->GetBinContent(bin); + const Double_t v3 = hist->GetBinContent(bin+1); + if (v1<=0 || v2<=0 || v3<=0) + continue; + + const Double_t y1 = TMath::Log(v1); + const Double_t y2 = TMath::Log(v2); + const Double_t y3 = TMath::Log(v3); + + //Double_t a = y2; + const Double_t b = (y3-y1)/2; + const Double_t c = (y1+y3)/2 - y2; + if (c>=0) + continue; + + const Double_t w = -1./(2*c); + const Double_t dx = b*w; + + if (dx<-1 || dx>1) + continue; + + // y = exp( - (x-k)^2 / s^2 / 2 ) + // + // -2*s^2 * log(y) = x^2 - 2*k*x + k^2 + // + // a = (k/s0)^2/2 + // b = k/s^2 + // c = -1/(2s^2) --> s = sqrt(-1/2c) + + const Double_t binx = hist->GetBinCenter(bin); + const Double_t binw = hist->GetBinWidth(bin); + + //const Double_t p = hist->GetBinCenter(hist->GetMaximumBin()); + const Double_t p = binx + dx*binw; + + avg += p; + rms += p*p; + + cnt++; + + // Store baseline and sigma + MPedestalPix &pix = (*fPedestal)[x]; + + pix.SetPedestal(p); + pix.SetPedestalRms(TMath::Sqrt(w)*binw); + + delete hist; + } + + avg /= cnt; + rms /= cnt; + + cout << "Baseline(" << cnt << "): " << avg << " +- " << sqrt(rms-avg*avg) << endl; + + return kTRUE; + } + + // Draw histogram + void Draw(Option_t *) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->SetBorderMode(0); + pad->SetFrameBorderMode(0); + fBaseline.Draw("colz"); + } + + + ClassDef(MHBaseline, 1); +}; + +// Histogram class for the signal and time distribution as +// well as the pulse shape +class MHSingles : public MH +{ + TH2F fSignal; + TH2F fTime; + TProfile2D fPulse; + + UInt_t fNumEntries; + + MSingles *fSingles; //! + MPedestalSubtractedEvt *fEvent; //! + MBadPixelsCam *fBadPix; //! + +public: + MHSingles() : fNumEntries(0), fSingles(0), fEvent(0) + { + fName = "MHSingles"; + + // Setup histograms + fSignal.SetName("Signal"); + fSignal.SetTitle("pulse integral spectrum"); + fSignal.SetXTitle("pixel [SoftID]"); + fSignal.SetYTitle("time [a.u.]"); + fSignal.SetDirectory(NULL); + + fTime.SetName("Time"); + fTime.SetTitle("arival time spectrum"); + fTime.SetXTitle("pixel [SoftID]"); + fTime.SetYTitle("time [a.u.]"); + fTime.SetDirectory(NULL); + + fPulse.SetName("Pulse"); + fPulse.SetTitle("average pulse shape spectrum"); + fPulse.SetXTitle("pixel [SoftID]"); + fPulse.SetYTitle("time [a.u.]"); + fPulse.SetDirectory(NULL); + } + + Bool_t SetupFill(const MParList *plist) + { + fSingles = (MSingles*)plist->FindObject("MSingles"); + if (!fSingles) + { + *fLog << err << "MSingles not found... abort." << endl; + return kFALSE; + } + + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fBadPix = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam"); + if (!fBadPix) + { + *fLog << err << "MBadPixelsCam not found... abort." << endl; + return kFALSE; + } + + fNumEntries = 0; + + return kTRUE; + } + + Bool_t ReInit(MParList *plist) + { + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << err << "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + // Setup binning + const Int_t w = fSingles->GetIntegrationWindow(); + + MBinning binsx, binsy; + binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5); + binsy.SetEdges(22*w, -10*w, 100*w); + + MH::SetBinning(fSignal, binsx, binsy); + + const UShort_t roi = header->GetNumSamples(); + + // Correct for DRS timing!! + MBinning binst(roi, -0.5, roi-0.5); + MH::SetBinning(fTime, binsx, binst); + + MBinning binspy(2*roi, -roi-0.5, roi-0.5); + MH::SetBinning(fPulse, binsx, binspy); + + return kTRUE; + } + + // Fill singles into histograms + Int_t Fill(const MParContainer *, const Stat_t) + { + // Get pointer to samples to fill samples + const Float_t *ptr = fEvent->GetSamples(0); + const Short_t roi = fEvent->GetNumSamples(); + + // Loop over all pixels + for (unsigned int i=0; iGetNumPixels(); i++) + { + if ((*fBadPix)[i].IsUnsuitable()) + continue; + + // loop over all singles + const vector &result = fSingles->GetVector(i); + + for (unsigned int j=0; jDivide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal.Draw("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fTime.Draw("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse.Draw("colz"); + } + + void DrawCopy(Option_t *) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->Divide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal.DrawCopy("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fTime.DrawCopy("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse.DrawCopy("colz"); + } + + ClassDef(MHSingles, 1) +}; + +// Task to extract the singles +class MExtractSingles : public MTask +{ + MSingles *fSingles; + MPedestalCam *fPedestal; + MPedestalSubtractedEvt *fEvent; + + // On time for each pixel in samples + MArrayI fExtractionRange; + + // Number of samples for sliding average + size_t fNumAverage; + + // The range in which the singles have to fit in + // FIXME: Make sure this is consistent with MHBaseline + UInt_t fStartSkip; + UInt_t fEndSkip; + + UInt_t fIntegrationSize; + UInt_t fMaxSearchWindow; + + Int_t fMaxDist; + Double_t fThreshold; + +public: + MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0), + fExtractionRange(1440), + fNumAverage(10), fStartSkip(20), fEndSkip(10), + fIntegrationSize(3*10), fMaxSearchWindow(30) + { + } + + void SetMaxDist(Int_t m) { fMaxDist = m; } + void SetThreshold(Int_t th) { fThreshold = th; } + + UInt_t GetIntegrationSize() const { return fIntegrationSize; } + + const MArrayI &GetExtractionRange() const { return fExtractionRange; } + + + Int_t PreProcess(MParList *plist) + { + fSingles = (MSingles*)plist->FindCreateObj("MSingles"); + if (!fSingles) + return kFALSE; + + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << err << "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fPedestal = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPedestal) + { + *fLog << err << "MPedestalCam not found... abort." << endl; + return kFALSE; + } + + return kTRUE; + } + + Int_t Process() + { + const UInt_t roi = fEvent->GetNumSamples(); + + const size_t navg = fNumAverage; + const float threshold = fThreshold; + const UInt_t start_skip = fStartSkip; + const UInt_t end_skip = fEndSkip; + const UInt_t integration_size = fIntegrationSize; + const UInt_t max_search_window = fMaxSearchWindow; + const UInt_t max_dist = fMaxDist; + + vector val(roi-navg);//result of first sliding average + for (int pix=0; pixGetNumPixels(); pix++) + { + // Total number of samples as 'on time' + fExtractionRange[pix] += roi-start_skip-navg-end_skip-integration_size; + + // Clear singles for this pixel + vector &result = fSingles->GetVector(pix); + result.clear(); + + // Get pointer to samples + const Float_t *ptr = fEvent->GetSamples(pix); + + // Get Baseline + const Double_t ped = (*fPedestal)[pix].GetPedestal(); + + // Subtract baseline and do a sliding average over + // the samples to remove noise (mainly the 200MHz noise) + for (size_t i=0; ithreshold || + val[i+4]>threshold || + + val[i+5] val[k_max]) + k_max = k; + } + + // Check if the findings make sense + if (k_max == i+5 || k_max == i + max_search_window-1) + continue; + + //search for half maximum before maximum + UInt_t k_half_max = 0; + for (UInt_t k=k_max; k>k_max-25; k--) + { + if (val[k-1] < val[k_max]/2 && + val[k] >= val[k_max]/2 ) + { + k_half_max = k; + break; + } + } + // Check if the finding makes sense + if (k_half_max+integration_size > roi-navg-end_skip) + continue; + if (k_half_max == 0) + continue; + if (k_max - k_half_max > max_dist) + continue; + + // FIXME: Evaluate arrival time more precisely!!! + // FIXME: Get a better integration window + + // Compile "single" + Single single; + single.fSignal = 0; + single.fTime = k_half_max; + + // Crossing of the threshold is at 5 + for (UInt_t j=k_half_max; jimax ? 5+integration_size-(i-imax) : 5+integration_size; + } + } + return kTRUE; + } +}; + +int extract_singles( + Int_t max_dist = 14, + Double_t threshold = 5, + const char *runfile = "/fact/raw/2012/07/23/20120723_", + int firstRunID = 6, + int lastRunID = 6, + const char *drsfile = "/fact/raw/2012/07/23/20120723_004.drs.fits.gz", + const char *outpath = "." + ) +{ + // ====================================================== + + MDirIter iter; + + // built outpu filename and path + TString outfile = Form("%s/", outpath); + outfile += gSystem->BaseName(runfile); + outfile += Form("%03d_%03d.root", firstRunID, lastRunID); + + + if (!gSystem->AccessPathName(outfile)) + return 0; + + // add input files to directory iterator + for (Int_t fileNr = firstRunID; fileNr <= lastRunID; fileNr++) + { + TString filename = runfile; + filename += Form("%03d.fits", fileNr); + iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".fz")); + iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(filename+".gz")); + } + + // ====================================================== + + // true: Display correctly mapped pixels in the camera displays + // but the value-vs-index plot is in software/spiral indices + // false: Display pixels in hardware/linear indices, + // but the order is the camera display is distorted + bool usemap = true; + + // map file to use (get that from La Palma!) + const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL; + + // ------------------------------------------------------ + + Long_t max0 = 0; + Long_t max1 = 0; + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + if (gSystem->AccessPathName(drsfile, kFileExists)) + { + gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl; + return 2; + } + + // -------------------------------------------------------------------------------- + + gLog.Separator("Singles"); + gLog << "Extract singles with '" << drsfile << "'" << endl; + gLog << endl; + + // ------------------------------------------------------ + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile)) + return 3; + + // ------------------------------------------------------ + + iter.Sort(); + iter.Print(); + + TString title; + title = iter.Next(); + title += "; "; + title += max1; + iter.Reset(); + + // ====================================================== + + MStatusDisplay *d = new MStatusDisplay; + + MBadPixelsCam badpixels; + badpixels.InitSize(1440); + badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); + + // Twin pixel + // 113 + // 115 + // 354 + // 423 + // 1195 + // 1393 + + MH::SetPalette(); + + MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); + + MGeomApply apply; + + MDrsCalibApply drsapply; + + MRawFitsRead read; + read.LoadMap(map); + read.AddFiles(iter); + + // ====================================================== + + gLog << endl; + gLog.Separator("Calculating baseline"); + + MTaskList tlist0; + + MRawRunHeader header; + MPedestalCam pedcam; + + MParList plist; + plist.AddToList(&tlist0); + plist.AddToList(&drscalib300); + plist.AddToList(&header); + plist.AddToList(&pedcam); + + // ------------------ Setup the tasks --------------- + + MFillH fill0("MHBaseline", "MPedestalSubtractedEvt"); + + drsapply.SetRemoveSpikes(3); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist0.AddToList(&read); + tlist0.AddToList(&apply); + tlist0.AddToList(&drsapply); + tlist0.AddToList(&fill0); + + // ------------------ Setup and run eventloop --------------- + + MEvtLoop loop0(title); + loop0.SetDisplay(d); + loop0.SetParList(&plist); + + if (!loop0.Eventloop(max0)) + return 4; + + if (!loop0.GetDisplay()) + return 5; + + // ---------------------------------------------------------------- + + MGeomCamFACT fact; + MHCamera hped(fact); + hped.SetCamContent(pedcam); + hped.SetCamError(pedcam, 4); + hped.SetAllUsed(); + MHCamEvent display; + display.SetHist(hped); + + d->AddTab("Baseline"); + display.Clone()->Draw(); + + // ====================================================== + + gLog << endl; + gLog.Separator("Extracting singles"); + + MTaskList tlist1; + + MSingles singles; + + plist.Replace(&tlist1); + plist.AddToList(&badpixels); + plist.AddToList(&singles); + + // ------------------ Setup the tasks --------------- + + MExtractSingles extract; + extract.SetMaxDist(max_dist); + extract.SetThreshold(threshold); + + MFillH fill("MHSingles"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&read); + tlist1.AddToList(&apply); + tlist1.AddToList(&drsapply); + tlist1.AddToList(&extract); + tlist1.AddToList(&fill); + + // ------------------ Setup and run eventloop --------------- + + MEvtLoop loop1(title); + loop1.SetDisplay(d); + loop1.SetParList(&plist); + + if (!loop1.Eventloop(max1)) + return 6; + + if (!loop1.GetDisplay()) + return 7; + + // ============================================================= + + gStyle->SetOptFit(kTRUE); + + MHSingles* h = (MHSingles*) plist.FindObject("MHSingles"); + if (!h) + return 8; + + const TH2 *htime = h->GetTime(); + const TH2 *hpulse = h->GetPulse(); + + d->AddTab("Time"); + TH1 *ht = htime->ProjectionY(); + ht->SetYTitle("Counts [a.u.]"); + ht->Scale(1./1440); + ht->Draw(); + + d->AddTab("Pulse"); + TH1 *hp = hpulse->ProjectionY(); + hp->SetYTitle("Counts [a.u.]"); + hp->Scale(1./1440); + hp->Draw(); + + d->SaveAs(outfile); + + TFile f(outfile, "UPDATE"); + + MParameterI par("NumEvents"); + par.SetVal(fill.GetNumExecutions()); + par.Write(); + + MParameterI win("IntegrationWindow"); + win.SetVal(extract.GetIntegrationSize()); + win.Write(); + + extract.GetExtractionRange().Write("ExtractionRange"); + + return 0; +} Index: /trunk/Mars/fact/analysis/gain/fit.sh =================================================================== --- /trunk/Mars/fact/analysis/gain/fit.sh (revision 17033) +++ /trunk/Mars/fact/analysis/gain/fit.sh (revision 17033) @@ -0,0 +1,28 @@ +inpath=$PWD/$1 +outpath=${inpath}/fit +macro=${PWD}/fit_spectra.C + +mkdir -p $outpath + +files_out=`ls -1 ${inpath} | grep \.root` +files_fit=`ls -1 ${outpath} | grep \.root` + +for file in $files_out +do + infile=$inpath"/"${file} + outfile=${outpath}"/"${file} + + if [ -e $outfile ] + then + continue + fi + + cmd="${ROOTSYS}/root -b -q -l ${macro}+\(\\\"${infile}\\\",\\\"${outfile}\\\"\)" + + cd ~/Mars + + echo ${cmd} | qsub -b yes -q test -cwd -e ${outfile}".err" -o ${outfile}".log" + + cd - +done + Index: /trunk/Mars/fact/analysis/gain/fit_spectra.C =================================================================== --- /trunk/Mars/fact/analysis/gain/fit_spectra.C (revision 17033) +++ /trunk/Mars/fact/analysis/gain/fit_spectra.C (revision 17033) @@ -0,0 +1,751 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "MStatusArray.h" +#include "MStatusDisplay.h" +#include "MHCamera.h" +#include "MGeomCamFACT.h" +#include "MParameters.h" +#include "MArrayI.h" + +using namespace std; + +// -------------------------------------------------------------------------- + +// Fit function for a single pe spectrum +Double_t fcn_g(Double_t *xx, Double_t *par) +{ + const Double_t ampl = par[0]; + const Double_t gain = par[1]; + const Double_t sigma = par[2]*gain; + const Double_t cross = par[3]; + const Double_t shift = par[4]; + const Double_t noise = par[5]; + const Double_t expo = par[6]; + + Double_t y = 0; + for (int N=1; N<14; N++) + { + const Double_t muN = N*gain + shift; + const Double_t sigN = TMath::Sqrt(N*sigma*sigma + noise*noise); + + const Double_t p = TMath::Power(cross, N-1) * TMath::Power(N, -expo); + + y += TMath::Gaus(xx[0], muN, sigN) * p / sigN; + } + + const Double_t sig1 = TMath::Sqrt(sigma*sigma + noise*noise); + return ampl*sig1*y; +} + +// Calculate the crosstalk from the function parameters +Double_t xtalk(TF1 &f) +{ + Double_t cross = f.GetParameter(3); + Double_t expo = f.GetParameter(6); + + Double_t y = 0; + for (int N=2; N<14; N++) + y += TMath::Power(cross, N-1) * TMath::Power(N, -expo); + + return y / (y + 1); +} + +// calculate the integral in units per millisecond +double integral(TF1 &func, TH1 &hist) +{ + const Double_t sigma = func.GetParameter(2)*func.GetParameter(1); + const Double_t cross = func.GetParameter(3); + const Double_t noise = func.GetParameter(5); + const Double_t expo = func.GetParameter(6); + + Double_t sum = 0; + for (int N=1; N<14; N++) + sum += TMath::Power(cross, N-1) * TMath::Power(N, -expo); + + const Double_t scale = hist.GetBinWidth(1); + const Double_t sig1 = TMath::Sqrt(sigma*sigma + noise*noise); + const Double_t integ = func.GetParameter(0)*sum*sig1*sqrt(TMath::TwoPi())/scale; + + return integ/1e-9/1e6; +} + + +// -------------------------------------------------------------------------- + +// Print function parameters +void PrintFunc(TF1 &f, float integration_window=30) +{ + //cout << "--------------------" << endl; + cout << "Ampl: " << setw(8) << f.GetParameter(0) << " +/- " << f.GetParError(0) << endl; + cout << "Gain: " << setw(8) << f.GetParameter(1) << " +/- " << f.GetParError(1) << endl; + cout << "Rel.sigma: " << setw(8) << f.GetParameter(2) << " +/- " << f.GetParError(2) << endl; + cout << "Baseline: " << setw(8) << f.GetParameter(4)/integration_window << " +/- " << f.GetParError(4)/integration_window << endl; + cout << "Crosstalk: " << setw(8) << f.GetParameter(3) << " +/- " << f.GetParError(3) << endl; + cout << "Pcrosstalk: " << setw(8) << xtalk(f) << endl; + cout << "Noise: " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl; + cout << "Expo: " << setw(8) << f.GetParameter(6) << " +/- " << f.GetParError(6) << endl; + //cout << "--------------------" << endl; + +} + +// -------------------------------------------------------------------------- + +// The parameters for the function are the filenam from the gain extraction +// and the output filename +int fit_spectra(const char *filename = "20130222_018_018.root", + const char *outfile = "20120802_247_247-fit.root" ) +{ + // It is more correct to fit the integral, but this is super + // slow, especially for 1440 pixel. To get that, one would have + // to analytically integrate and fit this function. + // (Currently the integral is switched off for the 1440 individual + // spectra in all cases) + bool fast = false; // Switch off using integral + + // Map for which pixel shall be plotted and which not + TArrayC usePixel(1440); + memset(usePixel.GetArray(), 1, 1440); + + // List of Pixel that should be ignored in camera view + usePixel[424] = 0; + usePixel[923] = 0; + usePixel[1208] = 0; + usePixel[583] = 0; + usePixel[830] = 0; + usePixel[1399] = 0; + usePixel[113] = 0; + usePixel[115] = 0; + usePixel[354] = 0; + usePixel[423] = 0; + usePixel[1195] = 0; + usePixel[1393] = 0; + + cout << setprecision(3); + + // ====================================================== + // Read data and histograms from file + + TFile file(filename); + if (file.IsZombie()) + { + cout << "Opening file '" << filename << "' failed." << endl; + return 1; + } + + MStatusArray arr; + if (arr.Read()<=0) + { + cout << "Reading of MStatusArray from '" << filename << "' failed." << endl; + return 2; + } + + TH2 *hsignal = (TH2*)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles"); + if (!hsignal) + { + cout << "Histogram Signal not found in '" << filename << "'." << endl; + return 3; + } + TH2 *htime = (TH2*)arr.FindObjectInCanvas("Time", "TH2F", "MHSingles"); + if (!htime) + { + cout << "Histogram Time not found in '" << filename << "'." << endl; + return 3; + } + TProfile2D *hpulse = (TProfile2D*)arr.FindObjectInCanvas("Pulse", "TProfile2D", "MHSingles"); + if (!hpulse) + { + cout << "Histogram Pulse not found in '" << filename << "'." << endl; + return 3; + } + TH2F *hbase = (TH2F*)arr.FindObjectInCanvas("Baseline", "TH2F", "MHBaseline"); + if (!hbase) + { + cout << "Histogram Baseline not found in '" << filename << "'." << endl; + return 3; + } + + MParameterI par("NumEvents"); + if (par.Read()<=0) + { + cout << "NumEvents not found in '" << filename << "'." << endl; + return 4; + } + + MParameterI win("IntegrationWindow"); + if (win.Read()<=0) + { + cout << "IntegrationWindow not found in '" << filename << "'." << endl; + return 5; + } + + Int_t integration_window = win.GetVal(); + + MArrayI ext; + if (ext.Read("ExtractionRange")<=0) + { + cout << "ExtractionRange not found in '" << filename << "'." << endl; + return 6; + + } + + // ====================================================== + + MStatusDisplay *d = new MStatusDisplay; + + // Camera geometry for displays + MGeomCamFACT fact; + + // ------------------ Spectrum Fit --------------- + // Instantiate the display histograms + MHCamera cRate(fact); + MHCamera cGain(fact); + MHCamera cRelSigma(fact); + MHCamera cCrosstalk(fact); + MHCamera cBaseline(fact); + MHCamera cNoise(fact); + MHCamera cChi2(fact); + MHCamera cNormGain(fact); + MHCamera cFitProb(fact); + + // Set name and title for the histograms + cRate.SetNameTitle ("Rate", "Dark count rate"); + cGain.SetNameTitle ("Gain", "Gain distribution"); + cRelSigma.SetNameTitle ("RelSigma", "Rel. Sigma"); + cCrosstalk.SetNameTitle("Crosstalk", "Crosstalk probability"); + cBaseline.SetNameTitle ("Baseline", "Baseline per sample"); + cNoise.SetNameTitle ("Noise", "Noise per sample"); + cChi2.SetNameTitle ("Chi2", "\\Chi^2"); + cNormGain.SetNameTitle ("NormGain", "Normalized gain"); + cFitProb.SetNameTitle ("FitProb", "Root's fit probability"); + + // Instantiate 1D histograms for the distributions + TH1F hRate ("Rate", "Dark count rate", 200, 0, 10); + TH1F hGain ("Gain", "Gain distribution", 100, 0, 400); + TH1F hRelSigma ("RelSigma", "Rel. Sigma", 160, 0, 0.40); + TH1F hCrosstalk("Crosstalk", "Crosstalk probability", 120, 0, 0.30); + TH1F hBaseline ("Baseline", "Baseline per sample", 75, -7.5, 7.5); + TH1F hNoise ("Noise", "Noise per sample", 60, 0, 30); + TH1F hChi2 ("Chi2", "\\Chi^2", 200, 0, 4); + TH1F hNormGain ("NormGain", "Normalized gain", 51, 0.5, 1.5); + TH1F hFitProb ("FitProb", "FitProb distribution", 100, 0, 1); + + // Histigram for the sum of all spectrums + TH1F hSum("Sum1", "Signal spectrum of all pixels", + hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); + hSum.SetXTitle("pulse integral [mV sample]"); + hSum.SetYTitle("Counts"); + hSum.SetStats(false); + hSum.Sumw2(); + + // Histogram for the sum of all pixels excluding the ones with faulty fits + TH1F hSumClear("SumC", "Signal spectrum of all pixels", + hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); + hSumClear.SetXTitle("pulse integral [mV sample]"); + hSumClear.SetYTitle("Counts"); + hSumClear.SetStats(false); + hSumClear.SetLineColor(kBlue); + hSumClear.Sumw2(); + + // Arrival time spectrum of the extracted pulses + TH1F hTime("Time", "Arrival time spectrum", htime->GetNbinsY(), htime->GetYaxis()->GetXmin(), htime->GetYaxis()->GetXmax()); + hTime.SetXTitle("pulse arrival time [sample]"); + hTime.SetYTitle("Counts"); + hTime.SetStats(false); + + // average pulse shape of the extracted pulses + TH1F hPulse("Puls", "Average pulse", hpulse->GetNbinsY(), hpulse->GetYaxis()->GetXmin(), hpulse->GetYaxis()->GetXmax()); + hPulse.SetXTitle("pulse arrival time [sample]"); + hPulse.SetYTitle("Counts"); + hPulse.SetStats(false); + + // Spektrum for the normalized individual spectra + TH1F hSumScale("Sum2", "Signal spectrum of all pixels", 120, 0.05, 12.05); + hSumScale.SetXTitle("pulse integral [pe]"); + hSumScale.SetYTitle("Rate"); + hSumScale.SetStats(false); + hSumScale.Sumw2(); + + // define fit function for Amplitudespectrum + TF1 func("spektrum", fcn_g, 0, hSum.GetXaxis()->GetXmax(), 7); + func.SetNpx(2000); + func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "Expo"); + func.SetLineColor(kRed); + + //--------------------- fill sum spectrum -------------------------------- + + d->SetStatusLine("Calculating sum spectrum", 0); + + // Begin of Loop over Pixels + for (Int_t pixel = 0; pixel < hsignal->GetNbinsX(); pixel++) + { + //jump to next pixel if the current is marked as faulty + if (usePixel[pixel]==0) + continue; + + TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); + hSum.Add(hist); + delete hist; + } + + //----------------- get starting values ------------------------------- + + hSum.GetXaxis()->SetRangeUser(150, hSum.GetXaxis()->GetXmax()); + + const Int_t maxbin = hSum.GetMaximumBin(); + const Double_t maxpos = hSum.GetBinCenter(maxbin); + const Double_t binwidth = hSum.GetBinWidth(maxbin); + const Double_t ampl = hSum.GetBinContent(maxbin); + + double fwhmSum = 0; + + //Calculate full width half Maximum + for (int i=1; iGetXmax(); + + // ------------------- fit -------------------------------- + + //Fit and draw spectrum + func.SetParLimits(0, 0, 2*ampl); // Amplitude + func.SetParLimits(1, 0, 2*maxpos); // Gain + func.SetParLimits(2, 0, 1); // Sigma + func.SetParLimits(3, 0, 1); // Crosstalk + func.SetParLimits(5, 0, 150); // Noise + + func.SetParameter(0, ampl); // Amplitude + func.SetParameter(1, maxpos); // Gain + func.SetParameter(2, 0.1); // Sigma + func.SetParameter(3, 0.15); // Crosstalk + func.SetParameter(4, 0*integration_window); // Baseline + func.SetParameter(5, 1.*sqrt(integration_window)); // Noise + func.SetParameter(6, 0.4); // Expo + + func.SetRange(fitmin, fitmax); + hSum.Fit(&func, fast?"N0QSR":"IN0QSR"); + + Double_t res_par[7]; + func.GetParameters(res_par); + + // ------------------ display result ------------------------------- + + cout << "--------------------" << endl; + cout << "AmplEst: " << ampl << endl; + cout << "GainEst: " << maxpos << endl; + cout << "SigmaEst: " << sigma_est << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + + gROOT->SetSelectedPad(0); + TCanvas &c11 = d->AddTab("SumHist"); + c11.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + hSum.GetXaxis()->SetRange(); + hSum.DrawCopy("hist"); + func.DrawCopy("same"); + + // =================================================================== + // Gain Calculation for each Pixel + // =================================================================== + + // counter for number of processed pixel + int count_ok = 0; + + // Begin of Loop over Pixels + for (Int_t pixel=0; pixelGetNbinsX(); pixel++) + { + // User information + d->SetStatusLine(Form("Fitting pixel #%d", pixel), 0); + d->SetProgressBarPosition((pixel+1.)/hsignal->GetNbinsX(), 1); + + // Skip pixels known to be faulty + if (usePixel[pixel]==0) + continue; + + //Projectipon of a certain Pixel out of Ampl.Spectrum + TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); + hist->SetDirectory(0); + + if (hist->GetEntries()<100) + { + cout << pixel << " ...histogram empty." << endl; + usePixel[pixel] = 0; + delete hist; + continue; + } + + //Rebin Projection + hist->Rebin(2); + + // Fit range + hist->GetXaxis()->SetRangeUser(150, hist->GetXaxis()->GetXmax()); + + // Determine start values + const Int_t maxBin = hist->GetMaximumBin(); + const Double_t maxPos = hist->GetBinCenter(maxBin); + + const Double_t gain = res_par[1]; + const Double_t GainRMS = res_par[2]; + + const double fit_min = maxPos-GainRMS*gain*2.5; + const double fit_max = fitmax;//maxPos+gain*(maxOrder-0.5); + + TArrayD cpy_par(7, res_par); + + cpy_par[0] = hist->GetBinContent(maxBin); + cpy_par[1] = maxPos-res_par[4]; // correct position for avg baseline + + func.SetParameters(cpy_par.GetArray()); + func.SetParLimits(0, 0, 2*cpy_par[0]); + func.SetParLimits(1, 0, 2*cpy_par[1]); + + // For individual spectra, the average fit yields 1 anyway + func.FixParameter(6, 1); // Expo + + // ----------- Fit Pixels spectrum --------------- + + const TFitResultPtr rc = hist->Fit(&func, /*fast?*/"LLN0QS"/*:"LLIN0QS"*/, "", fit_min, fit_max); + + // ----------- Calculate quality parameter --------------- + + Int_t b1 = hist->GetXaxis()->FindFixBin(fit_min); + Int_t b2 = hist->GetXaxis()->FindFixBin(fit_max); + + Double_t chi2 = 0; + Int_t cnt = 0; + for (int i=b1; i<=b2; i++) + { + if (hist->GetBinContent(i)<1.5 || func.Eval(hist->GetBinCenter(i))<1.5) + continue; + + const Double_t y = func.Integral(hist->GetBinLowEdge(i), hist->GetBinLowEdge(i+1)); + const Double_t v = hist->GetBinContent(i)*hist->GetBinWidth(i); + + const Double_t chi = (v-y)/v; + + chi2 += chi*chi; + cnt ++; + } + + chi2 = cnt==0 ? 0 : sqrt(chi2/cnt); + + // ----------------- Fit result -------------------- + + const double fit_prob = rc->Prob(); + + const float fRate = integral(func, *hist)/(ext[pixel]*0.5); + const float fGain = func.GetParameter(1); + const float fGainRMS = func.GetParameter(2); + const float fCrosstlk = xtalk(func); + const float fOffset = func.GetParameter(4); + const float fNoise = func.GetParameter(5)/sqrt(integration_window); + + // Fill histograms with result values + cRate.SetBinContent( pixel+1, fRate); + cGain.SetBinContent( pixel+1, fGain); + cRelSigma.SetBinContent( pixel+1, fGainRMS); + cCrosstalk.SetBinContent( pixel+1, fCrosstlk); + cBaseline.SetBinContent( pixel+1, fOffset/integration_window); + cNoise.SetBinContent( pixel+1, fNoise); + cChi2.SetBinContent( pixel+1, chi2); + cNormGain.SetBinContent( pixel+1, fGain/gain); + cFitProb.SetBinContent( pixel+1, fit_prob); + + // ====================================================== + + // Try to determine faulty fits + + bool ok = int(rc)==0; + + // mark pixels suspicious with failed fit + if (!ok) + cout << pixel << " ...fit failed!" << endl; + + // mark pixels suspicious with negative GainRMS + if (fabs(fGain/gain-1)>0.3) + { + cout << pixel << " ...gain deviates more than 30% from sum-gain." << endl; + ok = 0; + } + + if (fabs(fOffset/integration_window)>3) + { + cout << pixel << " ...baseline deviates." << endl; + ok = 0; + } + + // cancel out pixel where the fit was not succsessfull + usePixel[pixel] = ok; + + // Plot pixel 0 and all faulty fits + if (pixel==0 || !ok) + { + TCanvas &c = d->AddTab(Form("Pix%d", pixel)); + c.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + hist->SetStats(false); + hist->SetXTitle("Extracted signal"); + hist->SetYTitle("Counts"); + hist->SetName(Form("Pix%d", pixel)); + hist->GetXaxis()->SetRange(); + hist->DrawCopy("hist")->SetDirectory(0); + func.DrawCopy("SAME")->SetRange(fit_min, fit_max); + + cout << "--------------------" << endl; + cout << "Pixel: " << pixel << endl; + cout << "fit prob: " << fit_prob << endl; + cout << "AmplEst: " << cpy_par[0] << endl; + cout << "GainEst: " << cpy_par[1] << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + } + + if (!ok) + { + delete hist; + continue; + } + + // Fill Parameters into histograms + hRate.Fill( fRate); + hGain.Fill( fGain); + hRelSigma.Fill( fGainRMS); + hCrosstalk.Fill( fCrosstlk); + hBaseline.Fill( fOffset/integration_window); + hNoise.Fill( fNoise); + hChi2.Fill( chi2); + hNormGain.Fill( fGain/gain); + hFitProb.Fill( fit_prob); + + // Fill sum spectrum + for (int b=1; b<=hist->GetNbinsX(); b++) + hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); + delete hist; + + // Because of the rebinning... + hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); + hSumClear.Add(hist); + delete hist; + + hist = htime->ProjectionY("proj", pixel+1, pixel+1); + hTime.Add(hist); + delete hist; + + hist = hpulse->ProjectionY("proj", pixel+1, pixel+1); + hPulse.Add(hist); + delete hist; + + count_ok++; + } + + //------------------fill histograms----------------------- + // Display only pixels used and with valid fits + + cRate.SetUsed(usePixel); + cGain.SetUsed(usePixel); + cRelSigma.SetUsed(usePixel); + cCrosstalk.SetUsed(usePixel); + cBaseline.SetUsed(usePixel); + cNoise.SetUsed(usePixel); + cChi2.SetUsed(usePixel); + cNormGain.SetUsed(usePixel); + cFitProb.SetUsed(usePixel); + + // -------------------------------------------------------- + // Display data + + TCanvas *canv = &d->AddTab("Cams"); + canv->Divide(4,2); + + canv->cd(1); + cRate.DrawCopy(); + + canv->cd(2); + cGain.DrawCopy(); + + canv->cd(3); + cRelSigma.DrawCopy(); + + canv->cd(4); + cCrosstalk.DrawCopy(); + + canv->cd(5); + cBaseline.DrawCopy(); + + canv->cd(6); + cNoise.DrawCopy(); + + canv->cd(7); + cFitProb.DrawCopy(); + + canv->cd(8); + cChi2.DrawCopy(); + + // -------------------------------------------------------- + + gStyle->SetOptFit(1); + + canv = &d->AddTab("Hists"); + canv->Divide(4,2); + + TH1 *hh = 0; + + canv->cd(1); + hh = hRate.DrawCopy(); + + canv->cd(2); + hh = hGain.DrawCopy(); + hh->Fit("gaus"); + + canv->cd(3); + hh = hRelSigma.DrawCopy(); + hh->Fit("gaus"); + + canv->cd(4); + hh = hCrosstalk.DrawCopy(); + hh->Fit("gaus"); + + canv->cd(5); + hh = hBaseline.DrawCopy(); + hh->Fit("gaus"); + + canv->cd(6); + hh = hNoise.DrawCopy(); + hh->Fit("gaus"); + + canv->cd(7); + gPad->SetLogy(); + hFitProb.DrawCopy(); + + canv->cd(8); + hChi2.DrawCopy(); + + // -------------------------------------------------------- + + canv = &d->AddTab("NormGain"); + canv->Divide(2,1); + + canv->cd(1); + cNormGain.SetMinimum(0.8); + cNormGain.SetMaximum(1.2); + cNormGain.DrawCopy(); + + canv->cd(2); + gPad->SetLogy(); + hh = hNormGain.DrawCopy(); + hh->Fit("gaus"); + + //------------------ Draw gain corrected sum specetrum ------------------- + gROOT->SetSelectedPad(0); + c11.cd(); + hSumClear.DrawCopy("hist same"); + + //-------------------- fit gain corrected sum spectrum ------------------- + + gROOT->SetSelectedPad(0); + TCanvas &c11b = d->AddTab("CleanHist"); + c11b.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + const Int_t maxbin1 = hSumClear.GetMaximumBin(); + const Double_t ampl1 = hSumClear.GetBinContent(maxbin1); + + func.SetParameters(res_par); + func.SetParLimits(0, 0, 2*ampl1); + func.SetParameter(0, ampl1); + func.ReleaseParameter(6); + + hSumClear.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumClear.DrawCopy(); + func.DrawCopy("same"); + + cout << "--------------------" << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + + //-------------------- fit gain corrected sum spectrum ------------------- + + gROOT->SetSelectedPad(0); + TCanvas &c12 = d->AddTab("GainHist"); + c12.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + const Int_t maxbin2 = hSumScale.GetMaximumBin(); + const Double_t ampl2 = hSumScale.GetBinContent(maxbin2); + + //Set fit parameters + Double_t par2[7] = + { + ampl2, 1, 0.1, res_par[3], 0, res_par[5]/res_par[1], res_par[6] + }; + + func.SetParameters(par2); + func.SetParLimits(0, 0, 2*ampl2); + func.FixParameter(1, 1); + func.FixParameter(4, 0); + + func.SetRange(0.62, 9); + hSumScale.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumScale.DrawCopy(); + func.DrawCopy("same"); + + cout << "--------------------" << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + + //--------fit gausses to peaks of gain corrected sum specetrum ----------- + + d->AddTab("ArrTime"); + gPad->SetGrid(); + hTime.DrawCopy(); + + // ----------------------------------------------------------------- + + d->AddTab("Pulse"); + gPad->SetGrid(); + hPulse.DrawCopy(); + + // ================================================================ + + cout << "saving results to rootfile" << endl; + d->SaveAs(outfile); + cout << "..success!" << endl; + + TFile f(outfile, "UPDATE"); + par.Write(); + ext.Write("ExtractionRange"); + + return 0; +} + Index: /trunk/Mars/fact/analysis/gain/runanalysis.sh =================================================================== --- /trunk/Mars/fact/analysis/gain/runanalysis.sh (revision 17033) +++ /trunk/Mars/fact/analysis/gain/runanalysis.sh (revision 17033) @@ -0,0 +1,54 @@ +#!/bin/bash + +script=$PWD/runsinglepe.sh +threshold=5 +maxdist=14 +macro=$PWD/extract_singles.C +outpath=$PWD/files + +mkdir -p $outpath + +where="(((fNumPedestalTrigger=3000 OR fNumPedestalTrigger=5000) AND fRunTypeKEY=2) OR (fNumPedestalTrigger=10000 AND fRunTypeKEY=3) OR fRunTypeKEY=17)" +wheredrs="(drs.fROI=single.fROI AND fDrsStep=2)" + +query="SELECT Concat(fNight, '_', fRunId, '_', (SELECT fRunId FROM RunInfo \`drs\` WHERE drs.fNight=single.fNight AND "$wheredrs" AND (single.fRunId-drs.fRunId)<10 limit 0,1)) AS num " +query=$query" FROM RunInfo \`single\` WHERE "$where + +runpairs=( `mysql -u factread --password=r3adfac! --host=lp-fact factdata -s -e "$query"` ) + +for runpair in ${runpairs[@]} +do + year=`echo $runpair | cut -c1-4` + month=`echo $runpair | cut -c5-6` + day=`echo $runpair | cut -c7-8` + night=`echo $runpair | cut -d_ -f1` + run=`echo $runpair | cut -d_ -f2` + drs=`echo $runpair | cut -d_ -f3` + + runnum=`printf %03d $run` + drsnum=`printf %03d $drs` + + path="/fact/raw/"$year"/"$month"/"$day + runfile=$path"/"$night"_" + drsfile=$path"/"$night"_"$drsnum.drs.fits.gz + log=$outpath"/"$night"_"$runnum + + outfile=$outpath/`basename $runfile`${runnum}_${runnum}.root + + if [ -e $outfile ] + then + continue + fi + + echo "Submitting "${night}"_"${runnum}"["${drsnum}"]" + + cd ~/Mars + + cmd="${ROOTSYS}/root -q -b ${macro}+\($maxdist,$threshold,\\\"$runfile\\\",$run,$run,\\\"$drsfile\\\",\\\"${outpath}\\\"\)" + + echo ${cmd} | qsub -b yes -q test -cwd -e ${log}".err" -o ${log}".log" + + cd - +done + +