Changeset 9362 for trunk/MagicSoft


Ignore:
Timestamp:
02/22/09 23:36:01 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r9361 r9362  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20 2009/02/23 Thomas Bretz
     21
     22   * mcorsika/MCorsikaRunHeader.cc:
     23     - switched off the impact max workaround - it doesn't work
     24     - initialize fImpact Max
     25     - print also number of used ATMEXT
     26
     27   * mhbase/MH.cc:
     28     - improved setting of binning from the parameter list in the
     29       n-D case
     30
     31   * mhbase/MH3.cc:
     32     - imporved output
     33
     34   * mhflux/MHCollectionArea.cc:
     35     - also collect the maximum impact while running for a nicer
     36       behaviour of the plots
     37
     38   * mhflux/MHThreshold.cc:
     39     - allow setting of a dedicated Threshold binning
     40
     41   * mimage/MHHillasExt.cc:
     42     - converted slope binning to deg
     43
     44   * mjobs/MJSimulation.cc:
     45     - write MMcCorsikaRunHeader
     46     - changed binnings
     47     - added new binning
     48     - now display the signal unscaled
     49     - show threshold and collection area only for data runs
     50
     51   * msim/MSimMMCS.cc:
     52     - make setting of directions depending on view-cone option
     53
     54   * msim/MSimPointingPos.[h,cc]:
     55     - added class description
     56     - removed obsolete fPointingCorsika
     57     - improved output
     58     - added option for homogenous distribution
     59
     60
     61
    2062 2009/02/20 Thomas Bretz
    2163

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc

    r9359 r9362  
    6666//
    6767MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title)
    68     : fNumObsLevel(0), fZdMin(0), fZdMax(-1), fAzMin(0), fAzMax(0),
    69     fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
     68    : fNumObsLevel(0), fImpactMax(-1), fZdMin(0), fZdMax(-1),
     69    fAzMin(0), fAzMax(0),  fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
    7070{
    7171    fName  = name  ? name  : "MCorsikaRunHeader";
    … …  
    183183    // g[93] angle between array x-direction and magnetic north
    184184
     185    fImpactMax = -1;
     186/*
    185187    // This is a trick to use CERARY for storage of the
    186188    // maximum simulated impact
    187     fImpactMax = -1;
    188189    if (TMath::Nint(g[84])==1 && TMath::Nint(g[85])==1 &&
    189190        TMath::Nint(g[88])==1 && TMath::Nint(g[89])==1 &&
    190191        g[86]==g[87])
    191192        fImpactMax = g[86];
    192 
     193 */
    193194    fWavelengthMin = g[94];        // Cherenkov bandwidth lower end in nm
    194195    fWavelengthMax = g[95];        // Cherenkov bandwidth upper end in nm
    … …  
    272273
    273274    if (fViewConeOuterAngle>0)
    274         *fLog << "ViewCone:       " << fViewConeInnerAngle << UTF8::kDeg << "-" << fViewConeOuterAngle << UTF8::kDeg << endl;
     275        *fLog << "ViewCone:       " << fViewConeInnerAngle << UTF8::kDeg << " - " << fViewConeOuterAngle << UTF8::kDeg << endl;
    275276
    276277    if (fZdMax>=0)
    … …  
    300301        *fLog << " CEFFIC";
    301302    if (Has(kAtmext))
    302         *fLog << " ATMEXT";
     303        *fLog << " ATMEXT" << GetNumAtmosphericModel();
    303304    if (Has(kRefraction))
    304305        *fLog << " AtmRefraction";

  • trunk/MagicSoft/Mars/mhbase/MH.cc

    r9315 r9362  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.45 2009-02-11 12:08:23 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.46 2009-02-22 23:36:00 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    698698
    699699    if (bins->IsDefault())
     700    {
     701        gLog << inf << "Object 'Binning" << name << "' [MBinning] is default... binning unchanged." << endl;
    700702        return kTRUE;
     703    }
    701704
    702705    SetBinning(h, bins);
    … …  
    707710{
    708711    const MBinning *binsx = (MBinning*)plist.FindObject("Binning"+x);
     712    const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
     713
     714    if (!binsx && !binsy)
     715    {
     716        gLog << inf << "Neither 'Binning" << x << "' nor 'Binning" << y;
     717        gLog << "' [MBinning] found... no binning applied." << endl;
     718        return kFALSE;
     719    }
     720
    709721    if (!binsx)
    710     {
    711         gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... no binning applied." << endl;
     722        gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... binning unchanged." << endl;
     723    if (!binsy)
     724        gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... binning unchanged." << endl;
     725
     726    if (binsx && binsx->IsDefault())
     727    {
     728        gLog << inf << "Object 'Binning" << x << "' [MBinning] is default... binning unchanged." << endl;
     729        binsx = 0;
     730    }
     731
     732    if (binsy && binsy->IsDefault())
     733    {
     734        gLog << inf << "Object 'Binning" << y << "' [MBinning] is default... binning unchanged." << endl;
     735        binsy = 0;
     736    }
     737
     738    MBinning binsxx, binsyy;
     739    binsxx.SetEdges(*h, 'x');
     740    binsyy.SetEdges(*h, 'y');
     741
     742    SetBinning(h, binsx?binsx:&binsxx, binsy?binsy:&binsyy);
     743
     744    return kTRUE;
     745}
     746
     747Bool_t MH::ApplyBinning(const MParList &plist, TString x, TString y, TString z, TH3 *h)
     748{
     749    const MBinning *binsx = (MBinning*)plist.FindObject("Binning"+x);
     750    const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
     751    const MBinning *binsz = (MBinning*)plist.FindObject("Binning"+z);
     752
     753    if (!binsx && !binsy && !binsz)
     754    {
     755        gLog << inf << "Neither 'Binning" << x << "', 'Binning" << y;
     756        gLog << "' nor 'Binning" << z << "' [MBinning] found... ";
     757        gLog << "no binning applied." << endl;
    712758        return kFALSE;
    713759    }
    714     const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
     760
     761    if (!binsx)
     762        gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... binning unchanged." << endl;
    715763    if (!binsy)
    716     {
    717         gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... no binning applied." << endl;
    718         return kFALSE;
    719     }
    720 
    721     if (binsx->IsDefault() && binsy->IsDefault())
    722         return kTRUE;
    723 
    724     // -------------------------
    725     /*
    726      MBinning binsx, binsy, binsz;
    727      binsx.SetEdges(fHist, 'x');
    728      binsy.SetEdges(fHist, 'y');
    729      binsz.SetEdges(fHist, 'z');
    730      */
    731     // -------------------------
    732 
    733     SetBinning(h, binsx, binsy);
    734 
    735     return kTRUE;
    736 }
    737 
    738 Bool_t MH::ApplyBinning(const MParList &plist, TString x, TString y, TString z, TH3 *h)
    739 {
    740     const MBinning *binsx = (MBinning*)plist.FindObject("Binning"+x);
    741     if (!binsx)
    742     {
    743         gLog << inf << "Object 'Binning" << x << "' [MBinning] not found... no binning applied." << endl;
    744         return kFALSE;
    745     }
    746     const MBinning *binsy = (MBinning*)plist.FindObject("Binning"+y);
    747     if (!binsy)
    748     {
    749         gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... no binning applied." << endl;
    750         return kFALSE;
    751     }
    752     const MBinning *binsz = (MBinning*)plist.FindObject("Binning"+z);
     764        gLog << inf << "Object 'Binning" << y << "' [MBinning] not found... binning unchanged." << endl;
    753765    if (!binsz)
    754     {
    755         gLog << inf << "Object 'Binning" << z << "' [MBinning] not found... no binning applied." << endl;
    756         return kFALSE;
    757     }
    758 
    759     if (binsx->IsDefault() && binsy->IsDefault() && binsz->IsDefault())
    760         return kTRUE;
    761 
    762     SetBinning(h, binsx, binsy, binsz);
     766        gLog << inf << "Object 'Binning" << z << "' [MBinning] not found... binning unchanged." << endl;
     767
     768    if (binsx && binsx->IsDefault())
     769    {
     770        gLog << inf << "Object 'Binning" << x << "' [MBinning] is default... binning unchanged." << endl;
     771        binsx = 0;
     772    }
     773
     774    if (binsy && binsy->IsDefault())
     775    {
     776        gLog << inf << "Object 'Binning" << y << "' [MBinning] is default... binning unchanged." << endl;
     777        binsy = 0;
     778    }
     779
     780    if (binsz && binsz->IsDefault())
     781    {
     782        gLog << inf << "Object 'Binning" << z << "' [MBinning] is default... binning unchanged." << endl;
     783        binsy = 0;
     784    }
     785
     786    MBinning binsxx, binsyy, binszz;
     787    binsxx.SetEdges(*h, 'x');
     788    binsyy.SetEdges(*h, 'y');
     789    binszz.SetEdges(*h, 'z');
     790
     791    SetBinning(h, binsx?binsx:&binsxx, binsy?binsy:&binsyy, binsz?binsz:&binszz);
     792
    763793    return kTRUE;
    764794}

  • trunk/MagicSoft/Mars/mhbase/MH3.cc

    r9302 r9362  
    583583                if (!binsx)
    584584                {
    585                     *fLog << err << dbginf << "Neither '" << bx << "' nor '" << fName << "' found... aborting." << endl;
     585                    *fLog << err << dbginf << "Neither '" << bx << "' nor 'Binning" << fName << "' found... aborting." << endl;
    586586                    return kFALSE;
    587587                }

  • trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

    r9359 r9362  
    211211void MHCollectionArea::GetImpactMax()
    212212{
    213     if (fHeader->GetImpactMax()<=fMcAreaRadius*100)
     213    if (fHeader->GetImpactMax()<=fMcAreaRadius*100 || fHeader->GetImpactMax()<0)
    214214        return;
    215215
    … …  
    454454Int_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
    455455{
     456    // This is not perfect because it selects the maximum impact only
     457    // from the selected events. Hoever, it will get overwritten
     458    // in finalize anyway.
     459    const Double_t impact = fMcEvt->GetImpact()*0.01; // cm->m
     460    if (impact>fMcAreaRadius)
     461        fMcAreaRadius = impact;
     462
    456463    const Double_t energy = fMcEvt->GetEnergy();
    457464    const Double_t theta  = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();

  • trunk/MagicSoft/Mars/mhflux/MHThreshold.cc

    r9301 r9362  
    9797        return kFALSE;
    9898    }
    99 /*
    100     MBinning binse;
    101     binse.SetEdges(fHEnergy, 'x');
    10299
    103     MBinning *bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
    104     if (bins)
    105         binse.SetEdges(*bins);
     100    ApplyBinning(*pl, "Threshold", &fHEnergy);
    106101
    107     binse.Apply(fHEnergy);
    108   */
    109102    return kTRUE;
    110103}

  • trunk/MagicSoft/Mars/mimage/MHHillasExt.cc

    r9350 r9362  
    130130
    131131    binsx.SetEdges(100,     0,  1.5);
    132     binsy.SetEdges(100, -0.04, 0.04);
     132    binsy.SetEdges(100,    -9,    9);
    133133    MH::SetBinning(&fHSlopeL, &binsx, &binsy);
    134134}

  • trunk/MagicSoft/Mars/mjobs/MJSimulation.cc

    r9359 r9362  
    215215{
    216216    // Common run headers
     217    write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
    217218    write.AddContainer("MCorsikaRunHeader",   "RunHeaders", kFALSE);
    218219    write.AddContainer("MRawRunHeader",       "RunHeaders");
    … …  
    346347    // -------------------------------------------------------------------
    347348
    348     MBinning binse( 100,     1,   100000, "BinningEnergy", "log");
    349     MBinning binss( 100,     1, 10000000, "BinningSize",   "log");
     349    MBinning binse( 100,     1,   100000, "BinningEnergy",    "log");
     350    MBinning binsth( 35,   0.9,   900000, "BinningThreshold", "log");
     351    MBinning binsee( 35,   0.9,   900000, "BinningEnergyEst", "log");
     352    MBinning binss( 100,     1, 10000000, "BinningSize",      "log");
    350353    MBinning binsi( 100,  -500,      500, "BinningImpact");
    351354    MBinning binsh( 150,     0,       50, "BinningHeight");
    … …  
    353356    MBinning binszd( 70,     0,       70, "BinningZd");
    354357    MBinning binsvc( 45,     0,        9, "BinningViewCone");
    355     MBinning binsew(150,     0,       50, "BinningEvtWidth");
     358    MBinning binsel(150,     0,       50, "BinningTotLength");
     359    MBinning binsew(150,     0,       15, "BinningMedLength");
    356360    MBinning binsd("BinningDistC");
    357361    MBinning binsd0("BinningDist");
    358362    MBinning binstr("BinningTrigPos");
    359363
     364    plist.AddToList(&binsee);
    360365    plist.AddToList(&binse);
     366    plist.AddToList(&binsth);
    361367    plist.AddToList(&binss);
    362368    plist.AddToList(&binsi);
    … …  
    365371    plist.AddToList(&binsaz);
    366372    plist.AddToList(&binsvc);
     373    plist.AddToList(&binsel);
    367374    plist.AddToList(&binsew);
    368375    plist.AddToList(&binstr);
    … …  
    380387
    381388    MH3 mhew("MPhotonStatistics.fLength");
    382     mhew.SetName("EvtWidth");
    383     mhew.SetTitle("Time between first and last photon hitting a detector");
     389    mhew.SetName("TotLength");
     390    mhew.SetTitle("Time between first and last photon hitting a detector;L [ns]");
    384391
    385392    MH3 mhed("MPhotonStatistics.fTimeMedDev");
    386     mhed.SetName("EvtMedDev;EvtWidth");
     393    mhed.SetName("MedLength");
     394    mhed.SetTitle("Median deviation (1\\sigma);L [ns]");
    387395
    388396    MFillH fillh1(&mhn1, "", "FillCorsika");
    … …  
    488496    //  NSB = 0.2*6e-3
    489497
     498    // FIXME: Simulate photons before QE!
     499
    490500    MSimRandomPhotons  simnsb;
    491501//    simnsb.SetFreq(0.005, 0.004);    // 5MHz/cm^2, 4MHz
    … …  
    530540    hcalc.Disable(MHillasCalc::kCalcConc);
    531541
    532     MHCamEvent evt0a(/*10*/0, "Signal",    "Average signal;;S [ph]");
    533     MHCamEvent evt0c(/*10*/0, "MaxSignal", "Maximum signal;;S [ph]");
     542    MHCamEvent evt0a(/*10*/3, "Signal",    "Average signal (absolute);;S [ph]");
     543    MHCamEvent evt0c(/*10*/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
    534544    MHCamEvent evt0d(/*11*/8, "ArrTm",     "Time after first photon;;T [ns]");
    535545    evt0a.SetErrorSpread(kFALSE);
    … …  
    653663        //tasks.AddToList(&fillx2);
    654664        tasks.AddToList(&fillx3);
     665        //tasks.AddToList(&fillx4);
     666        //tasks.AddToList(&fillx5);
     667    }
     668    if (header.IsDataRun())
     669    {
    655670        tasks.AddToList(&fillth);
    656671        tasks.AddToList(&fillca);
    657         //tasks.AddToList(&fillx4);
    658         //tasks.AddToList(&fillx5);
    659672    }
    660673    //-------------------------------------------

  • trunk/MagicSoft/Mars/msim/MSimMMCS.cc

    r9359 r9362  
    174174    // Convert from corsika frame to telescope frame, taking
    175175    // the magnetic field into account: tel = corsika+offset
    176     fMcEvt->SetTheta(fEvtHeader->GetZd());
    177     fMcEvt->SetPhi(fEvtHeader->GetAz()+fRunHeader->GetMagneticFieldAz());
     176    if (fRunHeader->HasViewCone())
     177    {
     178        fMcEvt->SetTheta(fPointingTel->GetZdRad());
     179        fMcEvt->SetPhi(fPointingTel->GetAzRad());
     180    }
     181    else
     182    {
     183        fMcEvt->SetTheta(fEvtHeader->GetZd());
     184        fMcEvt->SetPhi(fEvtHeader->GetAz()+fRunHeader->GetMagneticFieldAz());
     185    }
    178186
    179187    fMcEvt->SetEvtNumber(fEvtHeader->GetEvtNumber());

  • trunk/MagicSoft/Mars/msim/MSimPointingPos.cc

    r9354 r9362  
    2525//////////////////////////////////////////////////////////////////////////////
    2626//
    27 //  MSimPointingPos
     27// MSimPointingPos
     28//
    2829//
    2930// This task is meant to simulate the pointing position (mirror orientation).
    … …  
    3132// on the user request (e.g. Wobble mode).
    3233//
    33 // WARNING: Currently the telescope orientation is just fixed to the
    34 //          direction in the run-header (if a view cone was given) or
    35 //          the direction in the evt-header (if no view cone given)
     34//
     35// If fOffTragetDistance==0 the telescope is oriented depending on the
     36// view cone option. If a view cone was given the orientation is fixed to
     37// the main direction around the view cone was produced. If no view
     38// cone was given the telescope is oriented parallel to the shower axis.
     39//
     40// If no view cone option was given and off-target observations are switched
     41// on by setting fOffTargetDistance!=0 the poitnting position is calculated:
     42//
     43//  1) fOffTargetDistance < 0:
     44//     The pointing position is randomly distributed in a disk of radius
     45//     -fOffTragetDistance. fOffTargetDistance is silently ignored.
     46//
     47//  2) fOffTargetDistance > 0:
     48//     The pointing position is set to a position in the given distance
     49//     away from the shower axis. If fOffTargetPhi>=0 it is fixed at
     50//     this phi value. For phi<0 it is randomly distributed at distances
     51//     fOffTargetDistance. (phi==0 is the direction of positive theta)
     52//
    3653//
    3754//  Input Containers:
    … …  
    4158//  Output Containers:
    4259//   MPointingPos
    43 //   PointingCorsika [MPointingPos]
    4460//
    4561//////////////////////////////////////////////////////////////////////////////
    … …  
    6884//
    6985MSimPointingPos::MSimPointingPos(const char* name, const char *title)
    70     : fRunHeader(0), fEvtHeader(0), fPointingCorsika(0), fPointingLocal(0),
    71     fOffTargetDistance(-1), fOffTargetPhi(-1)
     86    : fRunHeader(0), fEvtHeader(0), fPointing(0),
     87    fOffTargetDistance(0), fOffTargetPhi(-1)
    7288
    7389{
    … …  
    8399Int_t MSimPointingPos::PreProcess(MParList *pList)
    84100{
    85 //    fPointingCorsika = (MPointingPos*)pList->FindCreateObj("MPointingPos", "PointingCorsika");
    86 //    if (!fPointingCorsika)
    87 //        return kFALSE;
    88 
    89     fPointingLocal = (MPointingPos*)pList->FindCreateObj("MPointingPos");
    90     if (!fPointingLocal)
     101    fPointing = (MPointingPos*)pList->FindCreateObj("MPointingPos");
     102    if (!fPointing)
    91103        return kFALSE;
    92104
    … …  
    105117    }
    106118
    107     // FIXED offset
    108     // Diffuse?  ( FOV of camera folded with mirror diameter as Corsika input? )
    109     // Hour angle!
    110     // Angle to magnetic field!
    111 
    112     if (IsOffTargetObservation())
    113     {
    114         *fLog << inf;
    115         *fLog << "Off target observations switched on with" << endl;
     119    if (!IsOffTargetObservation())
     120        return kTRUE;
     121
     122    *fLog << inf;
     123    *fLog << "Off target observations switched on with" << endl;
     124    if (fOffTargetDistance>0)
     125    {
    116126        *fLog <<"   a pointing distance of " << GetOffTargetDistance() << "deg ";
    117127        if (fOffTargetPhi<0)
    … …  
    120130            *fLog << "and phi=" << GetOffTargetPhi() << "deg." << endl;
    121131    }
     132    else
     133        *fLog << "   a homogenous distribution up to a distance of " << -GetOffTargetDistance() << "deg " << endl;
    122134
    123135    return kTRUE;
    124136}
    125137
    126 /*
    127138Bool_t MSimPointingPos::ReInit(MParList *pList)
    128139{
     140    if (fRunHeader->HasViewCone() && IsOffTargetObservation())
     141    {
     142        *fLog << warn;
     143        *fLog << "WARNING - Combining the view cone option with off-target observations doesn't make sense." << endl;
     144        *fLog << "          Option for off-target observations will be ignored." << endl;
     145    }
    129146    // FIXME: Check also the enlightened region on the ground!
    130147    return kTRUE;
    131148}
    132 */
     149
     150void MSimPointingPos::GetDelta(Double_t &dtheta, Double_t &dphi) const
     151{
     152    if (fOffTargetDistance>0)
     153    {
     154        dtheta = fOffTargetDistance;
     155        dphi   = fOffTargetPhi>0 ? fOffTargetPhi : gRandom->Uniform(TMath::TwoPi());
     156    }
     157    else
     158    {
     159        dtheta = TMath::Sqrt(gRandom->Uniform(fOffTargetDistance));
     160        dphi   = gRandom->Uniform(TMath::TwoPi());
     161    }
     162}
    133163
    134164// --------------------------------------------------------------------------
    … …  
    136166Int_t MSimPointingPos::Process()
    137167{
    138 
    139168    // If a view cone is given use the fixed telescope orientation
    140     const Bool_t fixed = fRunHeader->HasViewCone();
     169    const Bool_t viewcone = fRunHeader->HasViewCone();
    141170
    142171    // Local sky coordinates (direction of telescope axis)
    143     /*const*/ Double_t zd = fixed ? fRunHeader->GetZdMin() : fEvtHeader->GetZd()*TMath::RadToDeg();  // x==north
    144     /*const*/ Double_t az = fixed ? fRunHeader->GetAzMin() : fEvtHeader->GetAz()*TMath::RadToDeg();  // y==west
    145 
    146     if (!fixed && IsOffTargetObservation())
    147     {
     172    Double_t zd = viewcone ? fRunHeader->GetZdMin() : fEvtHeader->GetZd()*TMath::RadToDeg();  // x==north
     173    Double_t az = viewcone ? fRunHeader->GetAzMin() : fEvtHeader->GetAz()*TMath::RadToDeg();  // y==west
     174
     175    if (!viewcone)
     176    {
     177        Double_t dtheta, dphi;
     178        GetDelta(dtheta, dphi);
     179
    148180        const Double_t theta = zd*TMath::DegToRad();
    149181        const Double_t phi   = az*TMath::DegToRad();
    150182
    151         /*const*/ TVector3 source;
    152         source.SetMagThetaPhi(1, theta, phi);
    153 
    154         /*const*/ TVector3 point;
    155         point.SetMagThetaPhi(1, theta+fOffTargetDistance, phi);
    156 
    157         const Double_t delta = fOffTargetPhi>0 ? fOffTargetPhi :
    158             gRandom->Uniform(TMath::TwoPi());
    159 
    160         point.Rotate(delta, source);
    161 
    162         zd = point.Theta()*TMath::RadToDeg();
    163         az = point.Phi()*TMath::RadToDeg();
    164     }
     183        TVector3 src, pnt;
     184        src.SetMagThetaPhi(1, theta,        phi);
     185        pnt.SetMagThetaPhi(1, theta+dtheta, phi);
     186
     187        pnt.Rotate(dphi, src);
     188
     189        zd = pnt.Theta()*TMath::RadToDeg();
     190        az = pnt.Phi()  *TMath::RadToDeg();
     191    }
     192
     193    // Transform the corsika coordinate system (north is magnetic north)
     194    // into the telescopes local coordinate system. Note, that all vectors
     195    // are already correctly oriented.
     196    az += fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg();
    165197
    166198    // Setup the pointing position
    167 //    fPointingCorsika->SetLocalPosition(zd, az/*+fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg()*/);
    168     fPointingLocal->SetLocalPosition(zd, az+fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg());
     199    fPointing->SetLocalPosition(zd, az);
    169200
    170201    // Calculate incident angle between magnetic field direction

  • trunk/MagicSoft/Mars/msim/MSimPointingPos.h

    r9336 r9362  
    1414{
    1515private:
    16     MCorsikaRunHeader *fRunHeader;         //! Header storing event information
    17     MCorsikaEvtHeader *fEvtHeader;         //! Header storing event information
    18     MPointingPos      *fPointingCorsika;   //! Output storing telescope poiting position in corsika coordinate system (modulo magnetig field declination)
    19     MPointingPos      *fPointingLocal;     //! Output storing telescope poiting position in local (telescope) coordinate system
     16    MCorsikaRunHeader *fRunHeader;  //! Header storing event information
     17    MCorsikaEvtHeader *fEvtHeader;  //! Header storing event information
     18    MPointingPos      *fPointing;   //! Output storing telescope poiting position in local (telescope) coordinate system
    2019
    21     Double_t fOffTargetDistance;  // [rad] Distance of the observed off-target position from the source
    22     Double_t fOffTargetPhi;       // [rad] Rotation angle of the off-target position (phi==0 means south, phi=90 west) [0;2pi], phi<0 means random
     20    Double_t fOffTargetDistance;    // [rad] Distance of the observed off-target position from the source
     21    Double_t fOffTargetPhi;         // [rad] Rotation angle of the off-target position (phi==0 means south, phi=90 west) [0;2pi], phi<0 means random
     22
     23    // MSimPointingPos
     24    void GetDelta(Double_t &dtheta, Double_t &dphi) const;
    2325
    2426    // MParContainer
    … …  
    2729    // MTask
    2830    Int_t  PreProcess(MParList *pList);
    29     //Bool_t ReInit(MParList *pList);
     31    Bool_t ReInit(MParList *pList);
    3032    Int_t  Process();
    3133
    … …  
    3436
    3537    // Getter
    36     Double_t GetOffTargetDistance() const { return fOffTargetDistance*TMath::RadToDeg(); }
     38    Double_t GetOffTargetDistance() const { return fOffTargetDistance==0 ? 0 : fOffTargetDistance*TMath::RadToDeg(); }
    3739    Double_t GetOffTargetPhi() const { return fOffTargetPhi*TMath::RadToDeg(); }
    3840
    3941    // Setter
    40     void SetOffTargetDistance(Double_t d=-1) { fOffTargetDistance = d*TMath::DegToRad(); }
    41     void SetOffTargetPhi(Double_t p=-1) { fOffTargetPhi = p*TMath::DegToRad(); }
     42    void SetOffTargetDistance(Double_t d=0) { fOffTargetDistance = d==0 ? 0 : d*TMath::DegToRad(); }
     43    void SetOffTargetPhi(Double_t p=0) { fOffTargetPhi = p*TMath::DegToRad(); }
    4244
    4345    // MSimPointingPos
    44     Bool_t IsOffTargetObservation() const { return fOffTargetDistance>0; }
     46    Bool_t IsOffTargetObservation() const { return fOffTargetDistance!=0; }
    4547
    4648    // TObject
Note: See TracChangeset for help on using the changeset viewer.