Changeset 19339


Ignore:
Timestamp:
Oct 30, 2018, 5:31:06 PM (3 months ago)
Author:
tbretz
Message:
Added the 8-float format as FillEvtCorsikaThin, calculate absolute wavelength (right now, I don't know what to do with the information from the negative WL, print the origin only if it is defined, added some docu from the corsika docu, removed MMCS specific stuff by a pre-compiler directive -- do we still have MMCS data?
Location:
trunk/Mars/msim
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/msim/MPhotonData.cc

    r18677 r19339  
    190190// contained in the input stream.
    191191//
    192 Int_t MPhotonData::FillRfl(Float_t f[8])
     192Int_t MPhotonData::FillRfl(const Float_t f[8])
    193193{
    194194    // Check coordinate system!!!!
     
    220220// system intpo our own.
    221221//
    222 Int_t MPhotonData::FillCorsika(Float_t f[7], Int_t i)
    223 {
    224     const UInt_t n = TMath::Nint(f[0]);
    225 
    226     if (n==0)
    227         // FIXME: Do we need to decode the rest anyway?
     222Int_t MPhotonData::FillCorsika(const Float_t f[7], Int_t i)
     223{
     224    // From the Corsika manual:
     225    //
     226    // f[0] : n Number of Cherenkov photons in bunch
     227    //          (For THIN option multiplied with thinning weight)
     228    // f[1] : x [cm]
     229    // f[2] : y [cm]
     230    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
     231    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
     232    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
     233    // f[6] : h [ch] bunch production height (except MCERFI=3)
     234    // f[7] : w weight of bunch [only with THIN option]
     235
     236    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
     237    // f[1] = XCER
     238    // f[2] = YXCER
     239    // f[3] = UEMIS
     240    // f[4] = VEMIS
     241    // f[5] = CARTIM
     242    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
     243    // #if __THIN__
     244    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
     245    // #endif
     246    //
     247    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
     248    //
     249    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
     250    //
     251    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
     252    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
     253    //
     254
     255    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
     256
     257    // Empty bunch (this happend at the end of events when
     258    // the remaining block is filled with zeroes)
     259    if (f[0]==0)
    228260        return kCONTINUE;
    229261
     262    if (f[0]!=1)
     263        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
     264
     265#ifdef __MMCS__
    230266    // Check reuse
    231267    if (i >=0)
    232       {
    233        const Int_t reuse = (n/1000)%100; // Force this to be 1!
    234        if (reuse!=i)
    235            return kCONTINUE;
    236       }
     268    {
     269        const Int_t reuse = (n/1000)%100; // Force this to be 1!
     270        if (reuse!=i)
     271        {
     272            cout << "REUSE " << reuse << " " << i << " " << n << endl;
     273            return kCONTINUE;
     274        }
     275    }
    237276
    238277    // This seems to be special to mmcs
    239278    fWavelength = n%1000;
    240279    fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
     280#else
     281    fWavelength = 0;
     282    fPrimary    = MMcEvtBasic::kUNDEFINED;
     283#endif
    241284
    242285    // x=north, y=west
     
    263306}
    264307
     308Int_t MPhotonData::FillCorsikaThin(const Float_t f[8], Int_t i)
     309{
     310    // From the Corsika manual:
     311    //
     312    // f[0] : n Number of Cherenkov photons in bunch
     313    //          (For THIN option multiplied with thinning weight)
     314    // f[1] : x [cm]
     315    // f[2] : y [cm]
     316    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
     317    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
     318    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
     319    // f[6] : h [ch] bunch production height (except MCERFI=3)
     320    // f[7] : w weight of bunch [only with THIN option]
     321
     322    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
     323    // f[1] = XCER
     324    // f[2] = YXCER
     325    // f[3] = UEMIS
     326    // f[4] = VEMIS
     327    // f[5] = CARTIM
     328    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
     329    // #if __THIN__
     330    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
     331    // #endif
     332    //
     333    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
     334    //
     335    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
     336    //
     337    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
     338    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
     339    //
     340
     341    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
     342
     343    // Empty bunch (this happend at the end of events when
     344    // the remaining block is filled with zeroes)
     345    if (f[0]==0)
     346        return kCONTINUE;
     347
     348    if (f[0]!=1)
     349        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
     350
     351
     352    // x=north, y=west
     353    //fPosX = f[1];  // [cm]
     354    //fPosY = f[2];  // [cm]
     355    //fCosU = f[3];  // cos to x
     356    //fCosV = f[4];  // cos to y
     357    // x=west, y=south
     358    fPosX =  f[2];  // [cm]
     359    fPosY = -f[1];  // [cm]
     360
     361    fCosU =  f[4];  // cos to x
     362    fCosV = -f[3];  // cos to y
     363
     364    fTime =  f[5];  // [ns]
     365
     366    fProductionHeight = f[6]; // [cm]
     367
     368    fWavelength = TMath::Abs(f[7]);
     369
     370    // Now reset all data members which are not in the stream
     371    fPrimary = MMcEvtBasic::kUNDEFINED;
     372    fTag     = -1;
     373    fWeight  =  1;
     374
     375    return kTRUE;
     376}
     377
    265378// --------------------------------------------------------------------------
    266379//
     
    272385// system into our own.
    273386//
    274 Int_t MPhotonData::FillEventIO(Short_t f[8])
     387Int_t MPhotonData::FillEventIO(const Short_t f[8])
    275388{
    276389    // From 5.5 compact_bunch:
     
    286399    fTime             =  f[4]/10.;              // a relative arival time [ns]
    287400    fProductionHeight =  pow(10, f[5]/1000.);   // altitude of emission a.s.l. [cm]
    288     fWavelength       =  f[7];                  // wavelength [nm]: 0 undetermined, <0 already in p.e.
     401    fWavelength       =  TMath::Abs(f[7]);      // wavelength [nm]: 0 undetermined, <0 already in p.e.
    289402
    290403    // Now reset all data members which are not in the stream
     
    305418// system into our own.
    306419//
    307 Int_t MPhotonData::FillEventIO(Float_t f[8])
     420Int_t MPhotonData::FillEventIO(const Float_t f[8])
    308421{
    309422    // photons in this bunch
     
    366479{
    367480    gLog << inf << endl;
    368 //    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
    369     gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
     481    //    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
     482    if (fPrimary!=MMcEvtBasic::kUNDEFINED)
     483        gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
    370484    gLog << "Wavelength:       " << fWavelength << "nm" << endl;
    371485    gLog << "Pos X/Y  Cos U/V: " << fPosX << "/" << fPosY << "   " << fCosU << "/" << fCosV << endl;
  • trunk/Mars/msim/MPhotonData.h

    r18549 r19339  
    134134    //Int_t ReadRflEvt(istream &fin);
    135135
    136     Int_t FillCorsika(Float_t f[7], Int_t i);
    137     Int_t FillEventIO(Short_t f[8]);
    138     Int_t FillEventIO(Float_t f[8]);
    139     Int_t FillRfl(Float_t f[8]);
     136    Int_t FillCorsika(const Float_t f[7], Int_t i);
     137    Int_t FillCorsikaThin(const Float_t f[8], Int_t i);
     138    Int_t FillEventIO(const Short_t f[8]);
     139    Int_t FillEventIO(const Float_t f[8]);
     140    Int_t FillRfl(const Float_t f[8]);
    140141
    141142    ClassDef(MPhotonData, 2) //Container to store a cherenkov photon bunch from a CORSUKA file
Note: See TracChangeset for help on using the changeset viewer.