Changeset 6190 for trunk/MagicSoft


Ignore:
Timestamp:
02/01/05 22:54:10 (20 years ago)
Author:
mazin
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mmpi
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc

    r5490 r6190  
    115115#include "MSkyPlot.h"
    116116
    117 #include <vector>
    118117#include <TF1.h>
    119118#include <TH2.h>
    … …  
    135134#include "MHillasExt.h"
    136135#include "MNewImagePar.h"
     136#include "MHadronness.h"
    137137#include "MTime.h"
    138138#include "MObservatory.h"
    … …  
    151151#include "MLogManip.h"
    152152
     153#include <TOrdCollection.h>
     154
    153155ClassImp(MSkyPlot);
    154156
    … …  
    160162//
    161163MSkyPlot::MSkyPlot(const char *name, const char *title)
    162     : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1)
    163 //    fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
    164 {
    165 
    166 *fLog << warn << "entering default constructor in MSkyPlot" << endl;
     164  : fGeomCam(NULL),   
     165    fTime(NULL),     
     166    fPointPos(NULL), 
     167    fRepDrive(NULL), 
     168    fSrcPosCam(NULL),
     169    fPntPosCam(NULL),
     170    fObservatory(NULL),
     171    fHillas(NULL),   
     172    fHillasExt(NULL),
     173    fHillasSrc(NULL),
     174    fNewImagePar(NULL),
     175    fMm2Deg(-1)
     176{
     177
     178  *fLog << warn << "entering default constructor in MSkyPlot" << endl;
    167179    //
    168180    //   set the name and title of this object
    … …  
    176188    fHistNexcess.SetDirectory(NULL);
    177189    fHistOn.SetDirectory(NULL);
     190    fHistSignifGaus.SetDirectory(NULL);
     191
     192    fHistSignif.UseCurrentStyle();
     193    fHistNexcess.UseCurrentStyle();
     194    fHistOn.UseCurrentStyle();
     195    fHistSignifGaus.UseCurrentStyle();
    178196
    179197    fHistSignif.SetName("SkyPlotSignif");
    … …  
    210228    fGridFineBin = 0.01;   // degrees
    211229
     230// elapsed time:
     231    fElaTime = 0.;
    212232
    213233    // some filter cuts:
    … …  
    217237    fMaxDist = 1.4;          // dist max cut (ever possible)
    218238    fMinDist = 0.1;          // dist max cut (ever possible)
     239    fHadrCut = 0.2;          // hadronness cut
    219240
    220241    fNumBinsAlpha   = 36;        // number of bins for alpha histograms
    … …  
    282303    }
    283304
    284     kSaveAlphaPlots=kTRUE;
    285     kSaveSkyPlots=kTRUE;
    286     kSaveNexPlot=kTRUE;
     305    fSaveAlphaPlots=kTRUE;
     306    fSaveSkyPlots=kTRUE;
     307    fSaveNexPlot=kTRUE;
     308    fUseRF=kFALSE;
    287309    fAlphaHName = "alphaplot.root";
    288310    fSkyHName   = "skyplot.root";
    289311    fNexHName   = "Nexcess.gif";
    290 }
     312
     313    fHistAlpha = new TOrdCollection();
     314    fHistAlpha->SetOwner();   
     315
     316}
     317
     318MSkyPlot::~MSkyPlot()
     319{
     320
     321  if (fHistAlpha)
     322    delete fHistAlpha;
     323}
     324
     325// --------------------------------------------------------------------------
     326//
     327// Get i-th hist
     328//
     329TH1D *MSkyPlot::GetAlphaPlot(Int_t i)
     330{
     331  if (GetSize() == 0)
     332    return NULL;
     333
     334  return static_cast<TH1D*>(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i));
     335}
     336
     337// --------------------------------------------------------------------------
     338//
     339// Get i-th hist
     340//
     341const TH1D *MSkyPlot::GetAlphaPlot(Int_t i) const
     342{
     343  if (GetSize() == 0)
     344    return NULL;
     345
     346  return static_cast<TH1D*>(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i));
     347}
     348
    291349
    292350void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
    … …  
    412470}
    413471
    414 
    415 
    416472// --------------------------------------------------------------------------
    417473//
    … …  
    423479Int_t MSkyPlot::PreProcess(MParList *plist)
    424480{
    425 *fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl;
     481
     482  *fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl;
     483 
    426484    fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
    427485    if (!fGeomCam)
    … …  
    438496
    439497
    440     fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
     498    fSrcPosCam = (MSrcPosCam*)plist->FindCreateObj(AddSerialNumber("MSrcPosCam"));
    441499    if (!fSrcPosCam)
    442500        *fLog << warn << "MSrcPosCam not found... no sky map." << endl;
    … …  
    477535        *fLog << err << "MNewImagePar not found... no sky map." << endl;
    478536
    479 
     537    if(fUseRF)
     538    {
     539       fHadron = (MHadronness*)plist->FindObject(AddSerialNumber("MHadronness"));
     540       if (!fHadron)
     541           *fLog << err << "MHadronness not found although specified... no sky map." << endl;
     542
     543       *fLog << inf << "Hadronness cut set to : " << fHadrCut << endl;
     544    }
    480545
    481546    // FIXME: Because the pointing position could change we must check
    … …  
    495560    fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid,
    496561                        fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    497    // fHistSignif.SetDirectory(NULL);
    498562    fHistSignif.SetFillStyle(4000);
    499     fHistSignif.UseCurrentStyle();
    500563
    501564   // fHistNexcess.SetName("SPNex2ndOrder");
    … …  
    503566    fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid,
    504567                        fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    505    // fHistNexcess.SetDirectory(NULL);
    506568    fHistNexcess.SetFillStyle(4000);
    507     fHistNexcess.UseCurrentStyle();
    508569
    509570   // fHistOn.SetName("SPOnCounted");
    … …  
    511572    fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid,
    512573                        fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
    513    // fHistOn.SetDirectory(NULL);
    514574    fHistOn.SetFillStyle(4000);
    515     fHistOn.UseCurrentStyle();
    516575
    517576    //fHistSignifGaus.SetName("SignifDistrib");
    … …  
    519578    fHistSignifGaus.SetBins(100, -10., 10.);
    520579
    521     // prepare alpha plots
     580    // create alpha histograms
     581    for (Int_t i=0; i< fNumalphahist; i++)
     582      {
    522583        // temp histogram for an alpha plot
    523     TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
    524     temp->SetDirectory(NULL);
    525     fHistAlphaBinWidth = temp->GetBinWidth(1);
    526     // create alpha histograms
    527     for (Int_t i=0; i< fNumalphahist; i++) {
    528         fHistAlpha.push_back(*temp);
    529         fHistAlpha[i].SetDirectory(NULL);
    530     }
    531 //cout  << " fHistAlpha[10].GetBinContent(5) " << fHistAlpha[10].GetBinContent(5) << endl;
    532 
    533     delete temp;
     584        TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
     585        temp->SetDirectory(NULL);
     586        fHistAlpha->AddAt(temp,i);
     587      }
     588
     589    fHistAlphaBinWidth = GetAlphaPlot()->GetBinWidth(1);
     590//cout  << " (*fHistAlpha)[10].GetBinContent(5) " << (*fHistAlpha)[10].GetBinContent(5) << endl;
     591
    534592    return kTRUE;
    535593}
    … …  
    541599Int_t MSkyPlot::Process()
    542600{
    543 // check whether MPointingPos comtains something:
     601
     602//      *fLog << err << "MPointingPos ENTERING the process" << endl;
     603//      *fLog << err << "MPointingPos:: fUseRF " << (int)fUseRF << endl;
     604  // check whether MPointingPos comtains something:
    544605  if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3)
    545   {
    546      *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
    547      return kTRUE;
    548   }
    549 
    550     // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
     606    {
     607      *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
     608      return kTRUE;
     609    }
     610 
     611  // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
    551612  if (fSetCenter==kTRUE)
    552   {
    553      if (fRepDrive)
    554      {
    555          fRa0  = fRepDrive->GetRa();  // [h]
    556          fDec0 = fRepDrive->GetDec(); // [deg]
    557          if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE;  // temp!, should be changed
    558      }
    559      else
    560      {
    561         fRa0 = fPointPos->GetRa();
    562         fDec0 = fPointPos->GetDec();
    563      }
    564      *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl;
    565      fSetCenter=kFALSE;
    566   }
    567 
    568 
    569 // some filter cuts:
     613    {
     614      if (fRepDrive)
     615        {
     616          fRa0  = fRepDrive->GetRa();  // [h]
     617          fDec0 = fRepDrive->GetDec(); // [deg]
     618          if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE;  // temp!, should be changed
     619        }
     620      else
     621        {
     622          fRa0 = fPointPos->GetRa();
     623          fDec0 = fPointPos->GetDec();
     624        }
     625      *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl;
     626      fSetCenter=kFALSE;
     627    }
     628
     629  // some filter cuts:
    570630  if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax)
    571   {
    572 
    573     Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
    574     Double_t dx, dy, beta, tanbeta, alphapar, distpar;
    575     Double_t logsize, lgsize, lgsize2, dist2;
    576     const Double_t log3000 = log(3000.);
    577     const Float_t fDistOffset = 0.9;
    578 
    579     // Get camera rotation angle
    580     Double_t rho = 0;
    581     if (fTime && fObservatory && fPointPos)
    582     {
    583         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
    584 //*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() <<
    585 //                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
    586 
    587     // => coord. system: xtilde, ytilde with center in RA_0, DEC_0
    588 
    589 // Get Rot. Angle:
    590     cosgam = TMath::Cos(rho);
    591     singam = TMath::Sin(rho);
    592 // Get x_0, y_0 for RA_0,DEC_0 of the current event
    593     }
    594     else
    595     {
    596 //      rho = fPointPos->RotationAngle(*fObservatory);
    597       Double_t theta, phi, sing, cosg;
    598       theta = fPointPos->GetZd()*TMath::Pi()/180.;
    599       phi = fPointPos->GetAz()*TMath::Pi()/180.;
    600 //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
    601       fObservatory->RotationAngle(theta, phi, sing, cosg);       
    602       cosgam = cosg;
    603       singam = sing;
    604     }
    605    // if (fPointPos)
    606    //    rho = fPointPos->RotationAngle(*fObservatory);
    607 
    608 /*
    609 //TEMP
    610 //        theta = mcevt->GetTelescopeTheta();
    611     Double_t theta, phi, sing, cosg;
    612     theta = fPointPos->GetZd()*TMath::Pi()/180.;
    613     phi = fPointPos->GetAz()*TMath::Pi()/180.;
    614 //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
    615     fObservatory->RotationAngle(theta, phi, sing, cosg);
     631    {
    616632     
    617 // conclusion: diffference in rho = 7 deg
    618   *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
    619         << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
    620 
    621     Double_t a1 =  0.876627;
    622     Double_t a3 = -0.481171;
    623     Double_t thetaTel=theta, phiTel=phi;
    624 
    625     Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
    626                               ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
    627                               ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
    628     Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
    629     Double_t sinal =    a1 * sin(phiTel) * denom;
    630 
    631   *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
    632         << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
    633 
    634 
    635 
    636 // END TEMP
    637 */
    638 
    639 // make the center of the plot different from the center of the camera
    640 /*
    641     x_0 = fPntPosCam->GetX()*fMm2Deg;
    642     y_0 = fPntPosCam->GetY()*fMm2Deg;
    643 */     
    644     x_0 = 0.; 
    645     y_0 = 0.;
    646 
    647     Int_t index = 0;  // index for alpha histograms
    648     // loop over xtilde
    649     for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
    650     {
    651         ysource = fMinYGrid + gridy * fBinStepGrid;
    652         // loop over ytilde
    653         for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
     633      Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
     634      Double_t dx, dy, beta, tanbeta, alphapar, distpar;
     635      Double_t logsize, lgsize, lgsize2, dist2, hadr;
     636      const Double_t log3000 = log(3000.);
     637      const Float_t fDistOffset = 0.9;
     638     
     639      //Get Hadronness if available:
     640      if(fUseRF)
    654641        {
    655 
    656            xsource = fMinXGrid + gridx * fBinStepGrid;
    657 
    658  /*     derotation    : rotate  into camera coordinates */
    659  /*     formel: (x_cam)    (cos(gam)  -sin(gam))   (xtilde)   (x_0)
    660                 (     )  = (                   ) * (      ) + (   )
    661                 (y_cam)    (sin(gam)   sin(gam))   (ytilde)   (y_0)
    662 */
    663            xsourcam = cosgam * xsource - singam * ysource + x_0;
    664            ysourcam = singam * xsource + cosgam * ysource + y_0;
    665 
    666 
    667  /*    end derotatiom    */
    668 //           xsourcam = xsource;
    669 //           ysourcam = ysource;
    670 
    671        /* calculate ALPHA und DIST according to the source position */
    672            dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
    673            dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
    674            tanbeta = dy / dx ;
    675            beta = TMath::ATan(tanbeta);
    676            alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
    677            distpar = sqrt( dy*dy + dx*dx );
    678            if(alphapar >  90.) alphapar -= 180.;
    679            if(alphapar < -90.) alphapar += 180.;
    680            alphapar = TMath::Abs(alphapar);
    681 
    682 // apply cuts
    683            logsize = log(fHillas->GetSize());
    684            lgsize = logsize-log3000;
    685            lgsize2 = lgsize*lgsize;
    686            dist2 = distpar*distpar - fDistOffset*fDistOffset;
    687 
    688            if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
    689                 (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
    690            if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
    691                 (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
    692            if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
    693                       distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
    694                       distpar < fMaxDist && distpar > fMinDist)
    695            {
    696 // gamma candidates!
    697 //*fLog <<  "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl;
    698 //*fLog <<  "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl;
    699                 fHistAlpha[index].Fill(alphapar);
    700            }
    701            index++;
     642                hadr=fHadron->GetHadronness();
     643//              cout << " will use RF " << hadr << endl;
    702644        }
    703      }
    704   }
    705 
     645      // Get camera rotation angle
     646      Double_t rho = 0;
     647      if (fTime && fObservatory && fPointPos)
     648        {
     649          rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     650          //*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() <<
     651          //                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
     652         
     653          // => coord. system: xtilde, ytilde with center in RA_0, DEC_0
     654         
     655          // Get Rot. Angle:
     656          cosgam = TMath::Cos(rho);
     657          singam = TMath::Sin(rho);
     658          // Get x_0, y_0 for RA_0,DEC_0 of the current event
     659        }
     660      else
     661        {
     662          //    rho = fPointPos->RotationAngle(*fObservatory);
     663          Double_t theta, phi, sing, cosg;
     664          theta = fPointPos->GetZd()*TMath::Pi()/180.;
     665          phi = fPointPos->GetAz()*TMath::Pi()/180.;
     666          printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
     667          fObservatory->RotationAngle(theta, phi, sing, cosg);       
     668          cosgam = cosg;
     669          singam = sing;
     670        }
     671      // if (fPointPos)
     672      //    rho = fPointPos->RotationAngle(*fObservatory);
     673     
     674      /*
     675        //TEMP
     676        //        theta = mcevt->GetTelescopeTheta();
     677        Double_t theta, phi, sing, cosg;
     678        theta = fPointPos->GetZd()*TMath::Pi()/180.;
     679        phi = fPointPos->GetAz()*TMath::Pi()/180.;
     680        //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
     681        fObservatory->RotationAngle(theta, phi, sing, cosg);
     682       
     683        // conclusion: diffference in rho = 7 deg
     684        //  *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
     685        //        << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
     686       
     687        Double_t a1 =  0.876627;
     688        Double_t a3 = -0.481171;
     689        Double_t thetaTel=theta, phiTel=phi;
     690       
     691        Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
     692        ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
     693        ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
     694        Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
     695        Double_t sinal =    a1 * sin(phiTel) * denom;
     696       
     697        //  *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
     698        //        << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
     699       
     700        // END TEMP
     701      */
     702           
     703      // make the center of the plot different from the center of the camera
     704      /*
     705        x_0 = fPntPosCam->GetX()*fMm2Deg;
     706        y_0 = fPntPosCam->GetY()*fMm2Deg;
     707      */     
     708      x_0 = 0.; 
     709      y_0 = 0.;
     710     
     711      Int_t index = 0;  // index for alpha histograms
     712      // loop over xtilde
     713      for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
     714        {
     715          ysource = fMinYGrid + gridy * fBinStepGrid;
     716          // loop over ytilde
     717          for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
     718            {
     719             
     720              xsource = fMinXGrid + gridx * fBinStepGrid;
     721             
     722              /*     derotation    : rotate  into camera coordinates */
     723              /*     formel: (x_cam)      (cos(gam)  -sin(gam))   (xtilde)   (x_0)
     724                     (     )  = - (                   ) * (      ) + (   )
     725                (y_cam)      (sin(gam)   cos(gam))   (ytilde)   (y_0)
     726              */
     727              //xsourcam = - (cosgam * xsource - singam * ysource) + x_0;
     728              //ysourcam = - (singam * xsource + cosgam * ysource) + y_0;
     729             
     730             
     731              /*    end derotatiom    */
     732                         xsourcam = xsource;
     733                         ysourcam = ysource;
     734
     735              /* calculate ALPHA und DIST according to the source position */
     736              dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
     737              dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
     738              tanbeta = dy / dx ;
     739              beta = TMath::ATan(tanbeta);
     740              alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
     741              distpar = sqrt( dy*dy + dx*dx );
     742              if(alphapar >  90.) alphapar -= 180.;
     743              if(alphapar < -90.) alphapar += 180.;
     744              alphapar = TMath::Abs(alphapar);
     745             
     746              if(fUseRF)
     747              {
     748               
     749//              cout << " will use RF " << hadr << endl;
     750                if(hadr<fHadrCut  &&  distpar < fMaxDist && distpar > fMinDist)
     751                {
     752                    TH1D *hist = GetAlphaPlot(index);
     753                    hist->Fill(alphapar);
     754                }
     755              }
     756              else
     757              {
     758                // apply cuts
     759                logsize = log(fHillas->GetSize());
     760                lgsize = logsize-log3000;
     761                lgsize2 = lgsize*lgsize;
     762                dist2 = distpar*distpar - fDistOffset*fDistOffset;
     763             
     764                if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
     765                   (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
     766                        if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
     767                        (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
     768                          if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
     769                          distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
     770                          distpar < fMaxDist && distpar > fMinDist)
     771                        {
     772                      // gamma candidates!
     773                      //*fLog <<  "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl;
     774                      //*fLog <<  "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl;
     775                        TH1D *hist = GetAlphaPlot(index);
     776                        hist->Fill(alphapar);
     777                        }
     778                }
     779                index++;
     780            }
     781        }
     782    }
    706783    return kTRUE;
    707784}
    … …  
    724801    numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.;
    725802
    726     vector <TH1D>::iterator alpha_iterator;
    727803    int index2 = 0;  // index of the TH2F histograms
    728804
    729     alpha_iterator = fHistAlpha.begin();
     805    TOrdCollectionIter Next(fHistAlpha);
     806    TH1D *alpha_iterator = NULL;
    730807
    731808    fHistNexcess.Reset();
    … …  
    734811    fHistSignifGaus.Reset();
    735812   
    736     while( alpha_iterator != fHistAlpha.end() ) {
     813    while ( (alpha_iterator = (TH1D*)Next())) {
    737814         // find the bin in the 2dim histogram
    738          nrow    = index2/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
    739          ncolumn = index2%fHistOn.GetNbinsX()+1;   // column of the histogram (x)
     815         nrow    = index2/fHistOn.GetNbinsX() + 1;              // row of the histogram (y)
     816         ncolumn = index2%fHistOn.GetNbinsX()+1;                // column of the histogram (x)
    740817         //ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
    741818
    … …  
    747824
    748825         // fit parabola to a background region
    749          (*alpha_iterator).Fit(fitbgpar,"EQRN");  // NWR OK?????????????????????????
     826         alpha_iterator->Fit(fitbgpar,"EQRN");  // NWR OK?????????????????????????
    750827         // get Chi2
    751828         chisquarefit = fitbgpar->GetChisquare();
    752829         if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01);
    753830         else  *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl;
    754 
     831         
    755832         // estimate Noff from the fit:
    756833         Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * TMath::Power(fAlphaONMax,3.) +
    757                     (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
    758 
     834                     (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
     835         
    759836         Nex = Non - Noff_fit;
    760837
    761838         fHistNexcess.SetBinContent(ncolumn, nrow, Nex);        // fill in the fHistNexcess
    762 
     839         
    763840         if (Noff_fit<0.) Sign = 0.;
    764 //         else Sign = LiMa17(Non,Noff_fit,1.);
     841         //         else Sign = LiMa17(Non,Noff_fit,1.);
    765842         else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.);
    766 
     843         
    767844         fHistSignif.SetBinContent(ncolumn, nrow, Sign);        // fill in the fHistSignif
    768845         fHistSignifGaus.Fill(Sign);
    769 
    770          alpha_iterator++;
     846         
    771847         index2++;
    772848    }
    773 
    774 // fit with gaus
     849   
     850    // fit with gaus
    775851    fHistSignifGaus.Fit("gaus","N");
    776 
    777 
    778 //temp
    779 /*
    780    Double_t decl, hang;
    781    MStarCamTrans mstarc(*fGeomCam, *fObservatory);
    782    mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
    783 
    784 *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
    785 *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
    786 */
    787 //endtemp
    788 
    789 
    790 // save alpha plots:
    791 //  TString stri1 = "alphaplots.root";
    792   if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName);
    793 
    794 // save sky plots:
    795 //  TString stri2 = "skyplots.root";
    796   if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName);
    797 
    798 // save nex plot:
    799   if(kSaveNexPlot==kTRUE) SaveNexPlot(fNexHName);
    800 
    801   fHistAlpha.clear();
    802   return kTRUE;
     852   
     853   
     854    //temp
     855    /*
     856      Double_t decl, hang;
     857      MStarCamTrans mstarc(*fGeomCam, *fObservatory);
     858      mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
     859     
     860      *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
     861      *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
     862      */
     863    //endtemp
     864   
     865   
     866    // save alpha plots:
     867    //  TString stri1 = "alphaplots.root";
     868    if(fSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName);
     869   
     870    // save sky plots:
     871    //  TString stri2 = "skyplots.root";
     872    if(fSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName);
     873   
     874    // save nex plot:
     875    if(fSaveNexPlot==kTRUE) SaveNexPlot(fNexHName);
     876   
     877    fHistAlpha->Clear();
     878
     879    delete fitbgpar;
     880   
     881    return kTRUE;
    803882}
    804883
    … …  
    843922    gStyle->SetPadLeftMargin(0.13);
    844923
     924    Char_t timet[100];
     925
    845926    TH2D tmp = fHistNexcess;
    846927    tmp.SetMaximum(470);
    … …  
    861942    gPad->SetTheta(20);
    862943    fHistNexcess.Draw("lego2");
     944    TLatex tu(0.5,0.8,"elapsed time:");
     945    tu.Draw();
     946    sprintf(timet,"%.1f min",fElaTime);
     947    TLatex tut(0.5,0.7,timet);
     948    tut.Draw();
     949
    863950    can.Print(nexplotname); 
    864951//    can.Print("test.root");
    … …  
    867954void MSkyPlot::SaveSkyPlots(TString skyplotfilename)
    868955{
    869    TFile rootfile(skyplotfilename, "RECREATE",
    870                    "sky plots in some variations");
    871      fHistSignif.Write();
    872      fHistNexcess.Write();
    873      fHistOn.Write();
    874      fHistSignif.Write();
    875 
    876 // from Wolfgang:
    877    //--------------------------------------------------------------------
    878     // Dec-Hour lines
    879     Double_t rho, sinrho, cosrho;
    880     Double_t theta, phi, sing, cosg;
    881 // do I need it?
    882     if (fTime && fObservatory && fPointPos)
    883     {
    884         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
    885         sinrho=TMath::Sin(rho);
    886         cosrho=TMath::Cos(rho);
    887     }
    888 
    889        theta = fPointPos->GetZd()*TMath::Pi()/180.;
    890        phi = fPointPos->GetAz()*TMath::Pi()/180.;
    891 //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
    892        fObservatory->RotationAngle(theta, phi, sing, cosg);
    893 
    894 *fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
    895 *fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
    896        sinrho=sing;
    897        cosrho=cosg;
    898 
    899     Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
    900     Double_t gridbinning = fGridBinning;
    901     Double_t gridfinebin = fGridFineBin;
    902     Int_t    numgridlines = (Int_t)(4.0/gridbinning);
    903     Double_t aberr = 1.07;
    904     Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
    905 
    906     Double_t declin, hangle;  // [deg], [h]
    907     MStarCamTrans mstarc(*fGeomCam, *fObservatory);
    908     if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
    909     else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD!
    910 
    911 //    std::vector <TGraph> graphvec;
    912     TLatex *pix;
    913 
    914     Char_t tit[100];
    915     Double_t xtxt;
    916     Double_t ytxt;
    917 
    918     Double_t xlo = -534.0 * mmtodeg;
    919     Double_t xup =  534.0 * mmtodeg;
    920 
    921     Double_t ylo = -534.0 * mmtodeg;
    922     Double_t yup =  534.0 * mmtodeg;
    923 
    924     Double_t xx, yy;
    925 
    926 
    927     // direction for the center of the camera
    928     Double_t dec0 = declin;
    929     Double_t h0   = hangle*360./24.; //deg
    930     Double_t RaHOffset = fRepDrive->GetRa() - h0;
     956
     957  TFile rootfile(skyplotfilename, "RECREATE",
     958                 "sky plots in some variations");
     959  fHistSignif.Write();
     960  fHistNexcess.Write();
     961  fHistOn.Write();
     962  fHistSignif.Write();
     963   
     964  // from Wolfgang:
     965  //--------------------------------------------------------------------
     966  // Dec-Hour lines
     967  Double_t rho, sinrho, cosrho;
     968  Double_t theta, phi, sing, cosg;
     969  // do I need it?
     970  if (fTime && fObservatory && fPointPos)
     971    {
     972      rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     973      sinrho=TMath::Sin(rho);
     974      cosrho=TMath::Cos(rho);
     975    }
     976 
     977  theta = fPointPos->GetZd()*TMath::Pi()/180.;
     978  phi = fPointPos->GetAz()*TMath::Pi()/180.;
     979  //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
     980  fObservatory->RotationAngle(theta, phi, sing, cosg);
     981 
     982  *fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
     983  *fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
     984  sinrho=sing;
     985  cosrho=cosg;
     986
     987  Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
     988  Double_t gridbinning = fGridBinning;
     989  Double_t gridfinebin = fGridFineBin;
     990  Int_t    numgridlines = (Int_t)(4.0/gridbinning);
     991  Double_t aberr = 1.07;
     992  Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
     993 
     994  Double_t declin, hangle;  // [deg], [h]
     995  MStarCamTrans mstarc(*fGeomCam, *fObservatory);
     996  if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
     997  else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD!
     998 
     999  TLatex *pix;
     1000 
     1001  Char_t tit[100];
     1002  Double_t xtxt;
     1003  Double_t ytxt;
     1004 
     1005  Double_t xlo = -534.0 * mmtodeg;
     1006  Double_t xup =  534.0 * mmtodeg;
     1007 
     1008  Double_t ylo = -534.0 * mmtodeg;
     1009  Double_t yup =  534.0 * mmtodeg;
     1010 
     1011  Double_t xx, yy;
     1012 
     1013 
     1014  // direction for the center of the camera
     1015  Double_t dec0 = declin;
     1016  Double_t h0   = hangle*360./24.; //deg
     1017  //    Double_t RaHOffset = fRepDrive->GetRa() - h0;
    9311018    //trans.LocToCel(theta0, phi0, dec0, h0);
    932 
    933     gStyle->SetOptFit(0);
    934     gStyle->SetOptStat(0);
    935     gStyle->SetPalette(1);
    936     TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
    937     c1->Divide(2,2);
    938     c1->cd(1);
    939     fHistSignif.Draw("colz");
    940     c1->cd(2);
    941     fHistNexcess.Draw("colz");
    942     c1->cd(3);
    943     fHistOn.Draw("colz");
    944     c1->cd(4);
    945     gPad->SetLogy();
    946     fHistSignifGaus.Draw();
    947 
    948     //-----   lines for fixed dec   ------------------------------------
    949 
    950     const Int_t Ndec = numgridlines;
    951     Double_t dec[Ndec];
    952     Double_t ddec = gridbinning;
    953 
    954     Int_t nh = (Int_t)(4.0/gridfinebin);
    955     const Int_t Nh   = nh+1;
    956     Double_t h[Nh];
    957     Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
    958     if ( dh > 360.0/((Double_t)(Nh-1)) )
    959       dh = 360.0/((Double_t)(Nh-1));
    960 
    961 // start to copy
    962     for (Int_t j=0; j<Ndec; j++)
     1019 
     1020  gStyle->SetOptFit(0);
     1021  gStyle->SetOptStat(0);
     1022  gStyle->SetPalette(1);
     1023  TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
     1024  c1->Divide(2,2);
     1025  c1->cd(1);
     1026  fHistSignif.Draw("colz");
     1027  c1->cd(2);
     1028  fHistNexcess.Draw("colz");
     1029  c1->cd(3);
     1030  fHistOn.Draw("colz");
     1031  c1->cd(4);
     1032  gPad->SetLogy();
     1033  fHistSignifGaus.Draw();
     1034 
     1035  //-----   lines for fixed dec   ------------------------------------
     1036 
     1037  const Int_t Ndec = numgridlines;
     1038  Double_t dec[Ndec];
     1039  Double_t ddec = gridbinning;
     1040 
     1041  Int_t nh = (Int_t)(4.0/gridfinebin);
     1042  const Int_t Nh   = nh+1;
     1043  Double_t h[Nh];
     1044  Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
     1045  if ( dh > 360.0/((Double_t)(Nh-1)) )
     1046    dh = 360.0/((Double_t)(Nh-1));
     1047 
     1048  // start to copy
     1049  for (Int_t j=0; j<Ndec; j++)
    9631050    {
    9641051      dec[j] = dec0 + ((Double_t)j)*ddec
    965                         - ((Double_t)(Ndec/2)-1.0)*ddec;
    966     }
    967 
    968     for (Int_t k=0; k<Nh; k++)
     1052        - ((Double_t)(Ndec/2)-1.0)*ddec;
     1053    }
     1054 
     1055  for (Int_t k=0; k<Nh; k++)
    9691056    {
    9701057      h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
    9711058    }
    972 
    973     Double_t xh[Nh];
    974     Double_t yh[Nh];
    975 
    976     for (Int_t j=0; j<Ndec; j++)
     1059 
     1060  Double_t xh[Nh];
     1061  Double_t yh[Nh];
     1062 
     1063  for (Int_t j=0; j<Ndec; j++)
    9771064    {
    9781065      if (fabs(dec[j]) > 90.0) continue;
    979 
     1066     
    9801067      for (Int_t k=0; k<Nh; k++)
    981       {
    982         Double_t hh0 = h0   *24.0/360.0;
    983         Double_t hx = h[k]*24.0/360.0;
    984         mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
    985                            xx, yy);
    986         xx = xx * mmtodeg * aberr;
    987         yy = yy * mmtodeg * aberr;
    988 //        xh[k] = xx * mmtodeg * aberr;
    989 //        yh[k] = yy * mmtodeg * aberr;
    990         xh[k] = cosg * xx + sing * yy;
    991         yh[k] =-sing * xx + cosg * yy;
    992 //        xh[k] = cosrho * xx + sinrho * yy;
    993 //        yh[k] =-sinrho * xx + cosrho * yy;
    994 
    995 
    996         //gLog << "dec0, h0 = " << dec0 << ",  " << h0
    997         //     << " : dec, h, xh, yh = " << dec[j] << ",  "
    998         //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
    999       }
    1000 
    1001 //      c1->cd(2);
     1068        {
     1069          Double_t hh0 = h0   *24.0/360.0;
     1070          Double_t hx = h[k]*24.0/360.0;
     1071          mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
     1072                              xx, yy);
     1073          xx = xx * mmtodeg * aberr;
     1074          yy = yy * mmtodeg * aberr;
     1075          //        xh[k] = xx * mmtodeg * aberr;
     1076          //        yh[k] = yy * mmtodeg * aberr;
     1077          xh[k] = cosg * xx + sing * yy;
     1078          yh[k] =-sing * xx + cosg * yy;
     1079          //        xh[k] = cosrho * xx + sinrho * yy;
     1080          //        yh[k] =-sinrho * xx + cosrho * yy;
     1081         
     1082         
     1083          //gLog << "dec0, h0 = " << dec0 << ",  " << h0
     1084          //     << " : dec, h, xh, yh = " << dec[j] << ",  "
     1085          //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
     1086        }
     1087     
     1088      //      c1->cd(2);
    10021089      TGraph * graph1 = new TGraph(Nh,xh,yh);
    10031090      //graph1->SetLineColor(j+1);
    … …  
    10171104      graph1->Draw("L");
    10181105      //delete graph1;
    1019 //      graphvec.push_back(*graph1);
    1020 //      graphvec[j].Draw("L");
    1021 
     1106      //      graphvec.push_back(*graph1);
     1107      //      graphvec[j].Draw("L");
     1108     
    10221109      sprintf(tit,"Dec = %6.2f", dec[j]);
    10231110      xtxt = xlo + (xup-xlo)*0.1;
    … …  
    10281115      //pix->Draw("");
    10291116      //delete pix;
    1030 
    1031     }
    1032 //stop copy
    1033     //-----   lines for fixed hour angle   ------------------------------------
    1034 
    1035     Int_t ndec1 = (Int_t)(4.0/gridfinebin);
    1036     const Int_t Ndec1 = ndec1;
    1037     Double_t dec1[Ndec1];
    1038     Double_t ddec1 = gridfinebin;
    1039 
    1040     const Int_t Nh1   = numgridlines;
    1041     Double_t h1[Nh1];
    1042     Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
    1043     if ( dh1 > 360.0/((Double_t)Nh1) )
    1044       dh1 = 360.0/((Double_t)Nh1);
    1045 
    1046 // start copy
    1047     for (Int_t j=0; j<Ndec1; j++)
     1117     
     1118    }
     1119  //stop copy
     1120  //-----   lines for fixed hour angle   ------------------------------------
     1121 
     1122  Int_t ndec1 = (Int_t)(4.0/gridfinebin);
     1123  const Int_t Ndec1 = ndec1;
     1124  Double_t dec1[Ndec1];
     1125  Double_t ddec1 = gridfinebin;
     1126 
     1127  const Int_t Nh1   = numgridlines;
     1128  Double_t h1[Nh1];
     1129  Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
     1130  if ( dh1 > 360.0/((Double_t)Nh1) )
     1131    dh1 = 360.0/((Double_t)Nh1);
     1132 
     1133  // start copy
     1134  for (Int_t j=0; j<Ndec1; j++)
    10481135    {
    10491136      dec1[j] = dec0 + ((Double_t)j)*ddec1
    1050                         - ((Double_t)(Ndec1/2)-1.0)*ddec1;
    1051     }
    1052 
    1053     for (Int_t k=0; k<Nh1; k++)
     1137        - ((Double_t)(Ndec1/2)-1.0)*ddec1;
     1138    }
     1139 
     1140  for (Int_t k=0; k<Nh1; k++)
    10541141    {
    10551142      h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
    10561143    }
    1057 
    1058 
    1059     Double_t xd[Ndec1];
    1060     Double_t yd[Ndec1];
    1061 
    1062     for (Int_t k=0; k<Nh1; k++)
     1144 
     1145  Double_t xd[Ndec1];
     1146  Double_t yd[Ndec1];
     1147 
     1148  for (Int_t k=0; k<Nh1; k++)
    10631149    {
    10641150      Int_t count = 0;
    10651151      for (Int_t j=0; j<Ndec1; j++)
    1066       {
    1067         if (fabs(dec1[j]) > 90.0) continue;
    1068 
    1069         Double_t hh0 = h0   *24.0/360.0;
    1070         Double_t hhx = h1[k]*24.0/360.0;
    1071         mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
     1152        {
     1153          if (fabs(dec1[j]) > 90.0) continue;
     1154         
     1155          Double_t hh0 = h0   *24.0/360.0;
     1156          Double_t hhx = h1[k]*24.0/360.0;
     1157          mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
    10721158                           xx, yy);
    1073 //        xd[count] = xx * mmtodeg * aberr;
    1074 //        yd[count] = yy * mmtodeg * aberr;
    1075 
    1076         xx = xx * mmtodeg * aberr;
    1077         yy = yy * mmtodeg * aberr;
    1078         xd[count] = cosg * xx + sing * yy;
    1079         yd[count] =-sing * xx + cosg * yy;
    1080 
    1081         //gLog << "dec0, h0 = " << dec0 << ",  " << h0
    1082         //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
    1083         //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;
    1084 
    1085         count++;
    1086       }
    1087 
    1088 //      c1->cd(2);
     1159          //        xd[count] = xx * mmtodeg * aberr;
     1160          //        yd[count] = yy * mmtodeg * aberr;
     1161         
     1162          xx = xx * mmtodeg * aberr;
     1163          yy = yy * mmtodeg * aberr;
     1164          xd[count] = cosg * xx + sing * yy;
     1165          yd[count] =-sing * xx + cosg * yy;
     1166
     1167          //gLog << "dec0, h0 = " << dec0 << ",  " << h0
     1168          //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
     1169          //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;
     1170         
     1171          count++;
     1172        }
     1173     
     1174      //      c1->cd(2);
    10891175      TGraph * graph1 = new TGraph(count,xd,yd);
    10901176      //graph1->SetLineColor(k+1);
    … …  
    11111197      //pix->Draw("");
    11121198      //delete pix;
    1113 
    1114     }
    1115 
    1116 //    c1->cd(2);
    1117     sprintf(tit,"Dec0 = %6.2f [deg]   Ra0 = %6.2f [h]", dec0, fRa0);
    1118     xtxt = xlo + (xup-xlo)*0.05 + 0.80;
    1119     ytxt = ylo + (yup-ylo)*0.75;
    1120     pix = new TLatex(xtxt, ytxt, tit);
    1121     pix->SetTextColor(1);
    1122     pix->SetTextSize(.05);
    1123     c1->cd(1);
    1124     pix->Draw("");
    1125     c1->cd(2);
    1126     pix->Draw("");
    1127     c1->cd(3);
    1128     pix->Draw("");
    1129     //delete pix;
    1130 
    1131     c1->Write();
    1132     // we suppose that the {skyplotfilename} ends with .root
    1133     Int_t sizeofout = skyplotfilename.Sizeof();
    1134     TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps";
    1135     c1->Print(outps);  // temporary!!!!!
    1136 
    1137     TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
    1138     c2->Divide(1,2);
    1139     c2->SetBorderMode(0);
    1140     c2->cd(1);
    1141     fHistSignif.Draw("colz");
    1142     c2->cd(2);
    1143     fHistNexcess.Draw("colz");
    1144     c2->Write();
    1145 
    1146    rootfile.Close();
    1147    delete c1;
    1148    delete c2;
    1149    delete pix;
    1150 
     1199     
     1200    }
     1201 
     1202  //    c1->cd(2);
     1203  sprintf(tit,"Dec0 = %6.2f [deg]   Ra0 = %6.2f [h]", dec0, fRa0);
     1204  xtxt = xlo + (xup-xlo)*0.05 + 0.80;
     1205  ytxt = ylo + (yup-ylo)*0.75;
     1206  pix = new TLatex(xtxt, ytxt, tit);
     1207  pix->SetTextColor(1);
     1208  pix->SetTextSize(.05);
     1209  c1->cd(1);
     1210  pix->Draw("");
     1211  c1->cd(2);
     1212  pix->Draw("");
     1213  c1->cd(3);
     1214  pix->Draw("");
     1215  //delete pix;
     1216 
     1217  c1->Write();
     1218  // we suppose that the {skyplotfilename} ends with .root
     1219  Int_t sizeofout = skyplotfilename.Sizeof();
     1220  TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps";
     1221  c1->Print(outps);  // temporary!!!!!
     1222
     1223  TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
     1224  c2->Divide(1,2);
     1225  c2->SetBorderMode(0);
     1226  c2->cd(1);
     1227  fHistSignif.Draw("colz");
     1228  c2->cd(2);
     1229  fHistNexcess.Draw("colz");
     1230  c2->Write();
     1231 
     1232  rootfile.Close();
     1233  delete c1;
     1234  delete c2;
     1235  delete pix;
    11511236
    11521237}
    … …  
    11571242                   "all the alpha plots");
    11581243
    1159     vector <TH1D>::iterator alpha_iterator;
    11601244    int index = 0;  // index of the TH2F histograms
    11611245    Char_t strtitle[100];
    … …  
    11641248    Int_t nrow, ncolumn, non;
    11651249
    1166 
    1167     alpha_iterator = fHistAlpha.begin();
    1168 
    1169     while( alpha_iterator != fHistAlpha.end() ) {
     1250    TH1D *alpha_iterator = NULL;
     1251    TOrdCollectionIter Next(fHistAlpha);
     1252
     1253    while( (alpha_iterator = (TH1D*)Next())) {
    11701254
    11711255          nrow    = index/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
    … …  
    11821266                        xpos, ypos, signif, nex, non);
    11831267           
    1184           (*alpha_iterator).SetName(strname);
    1185           (*alpha_iterator).SetTitle(strtitle);
    1186           (*alpha_iterator).Write();
    1187           alpha_iterator++;
     1268          alpha_iterator->SetName(strname);
     1269          alpha_iterator->SetTitle(strtitle);
     1270          alpha_iterator->Write();
    11881271          index++;
    11891272    }
    … …  
    11921275}
    11931276
    1194 // --------------------------------------------------------------------------
    1195 //
    1196 // Draw the histogram
    1197 //
    1198 void MSkyPlot::Draw(Option_t *opt)
    1199 {
    1200 /*
    1201     TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    1202     pad->SetBorderMode(0);
    1203 
    1204     AppendPad("");
    1205 
    1206     pad->Divide(1, 2, 0, 0.03);
    1207 
    1208     TObject *catalog = GetCatalog();
    1209 
    1210     // Initialize upper part
    1211     pad->cd(1);
    1212     gPad->SetBorderMode(0);
    1213     gPad->Divide(3, 1);
    1214 */
    1215 /*
    1216     // PAD #1
    1217     pad->GetPad(1)->cd(1);
    1218     gPad->SetBorderMode(0);
    1219     fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax);
    1220     TH1 *h3 = fHist.Project3D("yx_on");
    1221     fHist.GetZaxis()->SetRange(0,9999);
    1222     h3->SetDirectory(NULL);
    1223     h3->SetXTitle(fHist.GetXaxis()->GetTitle());
    1224     h3->SetYTitle(fHist.GetYaxis()->GetTitle());
    1225     h3->Draw("colz");
    1226     h3->SetBit(kCanDelete);
    1227     catalog->Draw("mirror same");
    1228 
    1229     // PAD #2
    1230     pad->GetPad(1)->cd(2);
    1231     gPad->SetBorderMode(0);
    1232     fHist.GetZaxis()->SetRange(0,0);
    1233     TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
    1234     fHist.GetZaxis()->SetRange(0,9999);
    1235     h4->SetTitle("Significance");
    1236     h4->SetDirectory(NULL);
    1237     h4->SetXTitle(fHist.GetXaxis()->GetTitle());
    1238     h4->SetYTitle(fHist.GetYaxis()->GetTitle());
    1239     h4->Draw("colz");
    1240     h4->SetBit(kCanDelete);
    1241     catalog->Draw("mirror same");
    1242 
    1243     // PAD #3
    1244     pad->GetPad(1)->cd(3);
    1245     gPad->SetBorderMode(0);
    1246     fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    1247     TH1 *h2 = fHist.Project3D("yx_off");
    1248     h2->SetDirectory(NULL);
    1249     h2->SetXTitle(fHist.GetXaxis()->GetTitle());
    1250     h2->SetYTitle(fHist.GetYaxis()->GetTitle());
    1251     h2->Draw("colz");
    1252     h2->SetBit(kCanDelete);
    1253     catalog->Draw("mirror same");
    1254 
    1255 
    1256     // Initialize lower part
    1257     pad->cd(2);
    1258     gPad->SetBorderMode(0);
    1259     gPad->Divide(3, 1);
    1260 
    1261     // PAD #4
    1262     pad->GetPad(2)->cd(1);
    1263     gPad->SetBorderMode(0);
    1264     TH1 *h1 = fHist.ProjectionZ("Alpha");
    1265     h1->SetDirectory(NULL);
    1266     h1->SetTitle("Distribution of \\alpha");
    1267     h1->SetXTitle(fHist.GetZaxis()->GetTitle());
    1268     h1->SetYTitle("Counts");
    1269     h1->Draw(opt);
    1270     h1->SetBit(kCanDelete);
    1271 
    1272     // PAD #5
    1273     pad->GetPad(2)->cd(2);
    1274     gPad->SetBorderMode(0);
    1275     TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
    1276     h5->Add(h2, -1);
    1277     h5->SetTitle("Difference of on- and off-distribution");
    1278     h5->SetDirectory(NULL);
    1279     h5->Draw("colz");
    1280     h5->SetBit(kCanDelete);
    1281     catalog->Draw("mirror same");
    1282 
    1283     // PAD #6
    1284     pad->GetPad(2)->cd(3);
    1285     gPad->SetBorderMode(0);
    1286     TH1 *h0 = fHist.Project3D("yx_all");
    1287     h0->SetDirectory(NULL);
    1288     h0->SetXTitle(fHist.GetXaxis()->GetTitle());
    1289     h0->SetYTitle(fHist.GetYaxis()->GetTitle());
    1290     h0->Draw("colz");
    1291     h0->SetBit(kCanDelete);
    1292     catalog->Draw("mirror same");
    1293 */
    1294 }
    1295 
    1296 // --------------------------------------------------------------------------
    1297 //
    1298 // Everything which is in the main pad belongs to this class!
    1299 //
    1300 Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py)
    1301 {
    1302     return 0;
    1303 }
    1304 
    1305 // --------------------------------------------------------------------------
    1306 //
    1307 // Set all sub-pads to Modified()
    1308 //
    1309 /*
    1310 void MSkyPlot::Modified()
    1311 {
    1312     if (!gPad)
    1313         return;
    1314 
    1315     TVirtualPad *padsave = gPad;
    1316     padsave->Modified();
    1317     padsave->GetPad(1)->cd(1);
    1318     gPad->Modified();
    1319     padsave->GetPad(1)->cd(3);
    1320     gPad->Modified();
    1321     padsave->GetPad(2)->cd(1);
    1322     gPad->Modified();
    1323     padsave->GetPad(2)->cd(2);
    1324     gPad->Modified();
    1325     padsave->GetPad(2)->cd(3);
    1326     gPad->Modified();
    1327     gPad->cd();
    1328 }
    1329 */
     1277
    13301278// --------------------------------------------------------------------------
    13311279//

  • trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h

    r5490 r6190  
    66#endif
    77
    8 #include <vector>
    9 
    108#ifndef ROOT_TH2
    119#include <TH2.h>
    … …  
    1412#ifndef ROOT_TH1
    1513#include <TH1.h>
     14#endif
     15
     16#ifndef ROOT_TOrdCollection
     17#include <TOrdCollection.h>
    1618#endif
    1719
    … …  
    2931class MHillasSrc;
    3032class MNewImagePar;
     33class MHadronness;
    3134class MGeomCam;
    32 
     35class TOrdCollection;
    3336class MSkyPlot : public MTask
    3437{
    3538private:
    36     MGeomCam      *fGeomCam;        //! container to take the event time from
     39
     40    MGeomCam      *fGeomCam;     //! container to take the event time from
    3741    MTime         *fTime;        //! container to take the event time from
    3842    MPointingPos  *fPointPos;    //! container to take pointing position from
    39     MReportDrive  *fRepDrive;     
    40     MSrcPosCam    *fSrcPosCam;    //! container with x and y of the source
    41     MSrcPosCam    *fPntPosCam;    //! container with x and y of the position MReportDrive.GetRa, MReportDrive.GetDec
     43    MReportDrive  *fRepDrive;    //!     
     44    MSrcPosCam    *fSrcPosCam;   //! container with x and y of the source
     45    MSrcPosCam    *fPntPosCam;   //! container with x and y of the position MReportDrive.GetRa, MReportDrive.GetDec
    4246    MObservatory  *fObservatory; //! container to take observatory location from
    43     MHillas       *fHillas;     
    44     MHillasExt    *fHillasExt;
    45     MHillasSrc    *fHillasSrc;
    46     MNewImagePar  *fNewImagePar;
     47    MHillas       *fHillas;      //!
     48    MHillasExt    *fHillasExt;   //!
     49    MHillasSrc    *fHillasSrc;   //!
     50    MNewImagePar  *fNewImagePar; //!
     51    MHadronness   *fHadron;      //!
    4752
    48     Float_t        fMm2Deg;             // conversion factor for display in degrees
    49     Double_t       fGridBinning;   // degrees
    50     Double_t       fGridFineBin;   // degrees
     53    TOrdCollection *fHistAlpha;  // vector of histograms for alpha
     54
     55    Float_t        fMm2Deg;      // conversion factor for display in degrees
     56    Double_t       fGridBinning; // degrees
     57    Double_t       fGridFineBin; // degrees
    5158
    5259//    Float_t fAlphaCut;           // Alpha cut
    … …  
    5966//    Float_t fMaxLD;              // Maximum distance in percent of dist
    6067
    61     std::vector <TH1D> fHistAlpha;           // vector of histograms for alpha
    6268    Int_t fNumalphahist;                // number of histograms for alpha
    6369    Int_t fNumBinsAlpha;
    … …  
    6571    Float_t fAlphaLeftEdge;
    6672    Float_t fAlphaRightEdge;
    67     Float_t fAlphaONMax;                        //  [deg] , upper cut for alpha ON region in the alpha plot
    68                       // [deg], ON region in the alpha plot, maybe 5 deg is better
    69                       // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
     73    Float_t fAlphaONMax;        //  [deg] , upper cut for alpha ON region in the alpha plot, [deg], ON region in the alpha plot, maybe 5 deg is better,  NOTE: up to now only values of 5, 10, 15, 20 degrees are possible
    7074    Float_t fAlphaBgLow;           // lower limit for bg region in the ON alpha plot
    7175    Float_t fAlphaBgUp;            // upper limit for bg region in the ON alpha plot
    7276   
    73     TH2D     fHistSignif;                // sky plot of significance vs. x and y
    74     TH2D     fHistNexcess;               // sky plot of number of excess events vs. x and y
    75     TH2D     fHistOn;                    // sky plot of events below fAlphaONMax vs. x and y
    76     TH1D     fHistSignifGaus;            // distribution of significance
    77     Bool_t   fSetCenter;                 // used to set the center of these histograms once
    78     Double_t fRa0;
    79     Double_t fDec0;
    80     Bool_t kSaveAlphaPlots;
    81     Bool_t kSaveSkyPlots;
    82     Bool_t kSaveNexPlot;
    83 
    84     Float_t fMinXGrid;                  //  [deg] , left edge of the skyplot
    85     Float_t fMaxXGrid;                  //  [deg] , right edge of the skyplot
    86     Float_t fMinYGrid;                  //  [deg] , upper edge of the skyplot
    87     Float_t fMaxYGrid;                  //  [deg] , lower edge of the skyplot
    88     Float_t fBinStepGrid;                       //  [deg] , grid size
    89     Int_t fNumStepsX;                   // number of bins in x
    90     Int_t fNumStepsY;                   // number of bins in y
     77    TH2D     fHistSignif;          // sky plot of significance vs. x and y
     78    TH2D     fHistNexcess;         // sky plot of number of excess events vs. x and y
     79    TH2D     fHistOn;              // sky plot of events below fAlphaONMax vs. x and y
     80    TH1D     fHistSignifGaus;      // distribution of significance
     81    Bool_t   fSetCenter;           // used to set the center of these histograms once
     82    Bool_t   fUseRF;               // use RF hadronness cut instead of supercuts
     83    Double_t fRa0;                 //   
     84    Double_t fDec0;                //
     85    Bool_t   fSaveAlphaPlots;      //
     86    Bool_t   fSaveSkyPlots;        //
     87    Bool_t   fSaveNexPlot;         //
     88                                   
     89    Float_t fMinXGrid;             //  [deg] , left edge of the skyplot
     90    Float_t fMaxXGrid;             //  [deg] , right edge of the skyplot
     91    Float_t fMinYGrid;             //  [deg] , upper edge of the skyplot
     92    Float_t fMaxYGrid;             //  [deg] , lower edge of the skyplot
     93    Float_t fBinStepGrid;          //  [deg] , grid size
     94    Int_t fNumStepsX;              // number of bins in x
     95    Int_t fNumStepsY;              // number of bins in y
    9196
    9297
    … …  
    97102    Float_t fMaxDist;          // dist max cut (ever possible)
    98103    Float_t fMinDist;          // dist min cut (ever possible)
     104    Float_t fHadrCut;          // hadronness cut
    99105
    100106    // coefficients for the cuts:
    … …  
    107113    TString fSkyHName;            // name for histogram with sky plots
    108114    TString fNexHName;            // name for canvas with Nex plot
    109 
    110     Int_t DistancetoPrimitive(Int_t px, Int_t py);
     115    Float_t fElaTime;             // elapsed time [min]
    111116
    112117    TObject *GetCatalog();
    … …  
    119124public:
    120125    MSkyPlot(const char *name=NULL, const char *title=NULL);
     126    ~MSkyPlot();
    121127
    122 //    TH1 *GetHistByName(const TString name) { return &fHist; }
    123     TH2D *GetHistSignif() { return &fHistSignif; }
    124     TH2D *GetHistNexcess() { return &fHistNexcess; }
    125     TH2D *GetHistOn() { return &fHistOn; }
     128    Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2);
     129
     130    TH2D *GetHistSignif    () { return &fHistSignif;  }
     131    TH2D *GetHistNexcess   () { return &fHistNexcess; }
     132    TH2D *GetHistOn        () { return &fHistOn;      }
    126133    TH1D *GetHistSignifGaus() { return &fHistSignifGaus; }
     134
     135    Int_t GetSize()  const { return fHistAlpha->GetSize(); }
     136
     137    TH1D *GetAlphaPlot( const Int_t i=-1);
     138    const TH1D *GetAlphaPlot( const Int_t=-1) const;
     139
     140    void ReadCuts(const TString parSCinit);
     141
     142    void SaveAlphaPlots(const TString stri2);
     143    void SaveNexPlot(const TString stri3);
     144    void SaveSkyPlots(TString stri);
     145
     146    void SetAlphaCut(Float_t alpha);
     147    void SetAlphaBGLimits(Float_t alphalow, Float_t alphalup);
    127148
    128149    void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance
    … …  
    130151    void SetSizeMin(Float_t size) { fSizeMin = size; } // Absolute minimum Size
    131152    void SetSizeMax(Float_t size) { fSizeMax = size; } // Absolute maximum Size
    132     void ReadCuts(const TString parSCinit);
    133153    void SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step);
    134     Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2);
     154    void SetHadrCut(Float_t b)    { fHadrCut = b;    }  // hadronness cut
    135155
    136156    void SetOutputSkyName(TString outname2)     { fSkyHName = outname2; }
    137     void SaveSkyPlots(TString stri);
    138     void SetNotSaveSkyPlots()                   { kSaveSkyPlots = kFALSE; }
    139 
     157    void SetNotSaveSkyPlots()                   { fSaveSkyPlots = kFALSE; }
    140158
    141159    void SetOutputAlphaName(TString outname1)   { fAlphaHName = outname1; }
    142     void SaveAlphaPlots(const TString stri2);
    143     void SetNotSaveAlphaPlots()                 { kSaveAlphaPlots = kFALSE; }
     160    void SetNotSaveAlphaPlots()                 { fSaveAlphaPlots = kFALSE; }
    144161
    145162    void SetOutputNexName(TString outname3)     { fNexHName = outname3; }
    146     void SaveNexPlot(const TString stri3);
    147     void SetNotSaveNexPlot()                    { kSaveNexPlot = kFALSE; }
     163    void SetElapsedTime(Float_t g)              { fElaTime = g; }
     164    void SetNotSaveNexPlot()                    { fSaveNexPlot = kFALSE; }
    148165
     166    void SetUseRF()                             { fUseRF = kTRUE; }
    149167
    150     void SetAlphaCut(Float_t alpha);
    151     void SetAlphaBGLimits(Float_t alphalow, Float_t alphalup);
    152 
    153     //std::vector <TH1D>::iterator  GetAlphaPlots()      { std::vector <TH1D>::iterator iter; return  iter = fHistAlpha.begin()  ; }
    154     std::vector <TH1D> GetAlphaPlots()      { return  fHistAlpha  ; }
    155  
    156     void Draw(Option_t *option="");
    157 
    158     ClassDef(MSkyPlot, 0) //2D-histogram in alpha, x and y
     168    ClassDef(MSkyPlot, 1) //2D-histogram in alpha, x and y
    159169};
    160170
Note: See TracChangeset for help on using the changeset viewer.