Index: trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C (revision 4416) +++ trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C (revision 4416) @@ -0,0 +1,359 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Luisa Sabrina Stark Schneebeli, 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! This macro is partially derived from findstars.C +\* ======================================================================== */ + + +Bool_t HandleInput() +{ + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + while (1) + { + // + // While reading the input process gui events asynchronously + // + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + + if (input=="q\n") + return kFALSE; + + if (input=="\n") + return kTRUE; + }; + + return kFALSE; +} + + +void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0) +{ + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + + MGeomCamMagic geomcam; + MCameraDC dccam; + MStarLocalCam starcam; + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&starcam); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddFile(directory+filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + + MGeomApply geomapl; + TString continuoslightfile = + "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; + MCalibrateDC dccal; + dccal.SetFileName(continuoslightfile); + + + const Int_t numblind = 12; + const Short_t x[numblind] = { 8, 27, + 507, 508, 509, 510, 511, + 543, 559, 560, 561, 567}; + const TArrayS blindpixels(numblind,(Short_t*)x); + Float_t ringinterest = 100; //[mm] + Float_t tailcut = 3.5; + UInt_t integratedevents = 5; + + MFindStars findstars; + // findstars.SetBlindPixels(blindpixels); + findstars.SetRingInterest(ringinterest); + findstars.SetDCTailCut(tailcut); + findstars.SetNumIntegratedEvents(integratedevents); + findstars.SetMinuitPrintOutLevel(0); + + findstars.SetSourceRaH(11); // Mrk421 + findstars.SetSourceRaM(4); + findstars.SetSourceRaS(26.0); + findstars.SetSourceDecD(38); + findstars.SetSourceDecM(12); + findstars.SetSourceDecS(36.0); + + + // prints + MPrint pdc("MCameraDC"); + MPrint pstar("MStarLocalCam"); + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&dccal); + // tlist.AddToList(&pdc, "Currents"); + tlist.AddToList(&findstars, "Currents"); + // tlist.AddToList(&pstar, "Currents"); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + const Int_t numpar = 3; + Float_t starpos[10*numpar]; + Float_t starposcat[10*numpar]; + + if (numEvents > 0) + { + if (!evtloop.Eventloop(numEvents)) + return; + } + else + { + if (!evtloop.PreProcess()) + return; + + MHCamera display(geomcam); + display.SetPrettyPalette(); + + TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600); + + TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2); + h0->GetXaxis()->SetTitle("(deg)"); + h0->GetYaxis()->SetTitle("(deg)"); + h0->Draw(); + + display.Draw(); + gPad->cd(1); + //display.Draw("same"); + starcam.Draw(); + + UShort_t year; + Byte_t month, day; + + MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents"); + + TList *starcatalog = new TList; + TList *stardc = new TList; + TVector3 *res = new TVector3; + + Int_t maxdist = 100; + UInt_t numevents=0; + Int_t numEvt=0; + const Int_t maxnumEvt = 10000; + Int_t col = 2; + const Double_t kDeg2mm = 189/0.6; + + Double_t resx[maxnumEvt]; + Double_t resy[maxnumEvt]; + Double_t time[maxnumEvt]; + + +// +// Find in the file name the date, run and project name + + Size_t pos = filename.Last('/'); + TString iRun = TString(filename(pos+24,5)); + TString iYear = TString(filename(pos+4,4)); + TString iMonth = TString(filename(pos+9,2)); + TString iDay = TString(filename(pos+12,2)); + + TString iHour = TString(filename(pos+15,2)); + TString iMin = TString(filename(pos+18,2)); + TString iSec = TString(filename(pos+21,2)); + + Size_t poslast = filename.Last('.'); + Size_t posfirst = poslast-1; + while (filename[posfirst] != '_') + posfirst--; + + TString iSource = TString(filename(posfirst+1,poslast-posfirst-1)); + char str[100]; + char str2[100]; + sprintf(str,"Date %s %s %s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + + UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec; + UShort_t evtmsec; + + ofstream fposout1; + // fposout1.open (str2, ofstream::app); + fposout1.open (str2, ofstream::trunc); + + while (tlist.Process()) + { + if (numevents >= integratedevents) + { + evttime->GetTime(evthour, evtmin, evtsec,evtmsec); + evttime->Print(); + evthour = evttime->Hour(); + evtmin = evttime->Min(); + evtsec = evttime->Sec(); + + // display.SetCamContent(dccam); + display.SetCamContent(findstars.GetDisplay()); + gPad->Modified(); + gPad->Update(); + Int_t numStars = starcam.GetNumStars(); + starcam.Print("bkmaxposrotsizechi"); + + res = findstars.GetCorrSrcPos(); + cout<<"Corrected source position = "<X()<<" , "<Y()<X(); + resy[numEvt] = res->Y(); + time[numEvt] = numEvt; + numEvt ++; + + // write the corrected source position and the star positions to file + if (!fposout1) + cout << "Could not open output file to write the center of the fit." << endl; + else{ + for ( Int_t i = 0; i < numStars; i++) + { + if (i == 0){ + fposout1 << str<<" "<X()<< " " << res->Y(); + } + fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY(); + } + } + fposout1 << endl; + numevents = 0; + + // Remove the comments if you want to go through the file + // event-by-event: + //if (!HandleInput()) + // break; + } + numevents++; + } + fposout1.close(); + + // Print source position + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptStat(00000000); + + TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200); + c1->Divide(1,3); + + TMath math; + + Double_t minresx, maxresx; + minresx = resx[math.LocMin(numEvt,resx)]; + maxresx = resx[math.LocMax(numEvt,resx)]; + + Double_t minresy, maxresy; + minresy = resy[math.LocMin(numEvt,resy)]; + maxresy = resy[math.LocMax(numEvt,resy)]; + + Double_t minres, maxres; + minres = math.Min(minresx,minresy); + maxres = math.Max(maxresx,maxresy); + + Double_t diff; + diff = maxres - minres; + diff = 0.1*diff; + minres = minres - diff; + maxres = maxres + diff; + + Double_t mintime, maxtime; + mintime = time[math.LocMin(numEvt,time)]; + maxtime = time[math.LocMax(numEvt,time)]; + + c1->cd(1); + TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres); + h1->GetXaxis()->SetTitle("(evt.nr)"); + h1->GetYaxis()->SetTitle("(mm)"); + h1->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + c1->cd(2); + TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx); + h2->GetXaxis()->SetTitle("(evt.nr)"); + h2->GetYaxis()->SetTitle("(mm)"); + h2->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.Draw(); + + c1->cd(3); + TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy); + h3->GetXaxis()->SetTitle("(evt.nr)"); + h3->GetYaxis()->SetTitle("(mm)"); + h3->Draw(); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + evtloop.PostProcess(); + } + + tlist.PrintStatistics(); + +} + +