Ignore:
Timestamp:
05/08/15 12:32:23 (10 years ago)
Author:
dbaack
Message:
 
Location:
branches/Corsika7405Compatibility
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

  • branches/Corsika7405Compatibility/Makefile

    r17980 r18179  
    7878          mmovie \
    7979          mtools \
    80           mmuon
     80          mmuon \
     81          mdebug
    8182
    8283LIBNOVA_PATH := `which libnovaconfig`

  • branches/Corsika7405Compatibility/mbase/MQuaternion.h

    r9564 r18179  
    99#ifndef ROOT_TQuaternion
    1010#include <math.h>
    11 #define sqrt ::sqrt
     11//#define sqrt ::sqrt
    1212#include <TQuaternion.h>
    1313#undef sqrt

  • branches/Corsika7405Compatibility/mcorsika/MCorsikaFormat.cc

    r10213 r18179  
    5757    }
    5858
    59     char *buffer = new char[5];
    60     memset(buffer, 0, 5);
    61     fileIn->read(buffer, 4);
    62     fileIn->seekg(-4, ios::cur);
     59    char *buffer = new char[9];
     60    memset(buffer, 0, 9);
     61    fileIn->read(buffer, 8);
     62    fileIn->seekg(-8, ios::cur);
    6363
    6464    if (strcmp(buffer, "RUNH") == 0)
    … …  
    6767        return new MCorsikaFormatRaw(fileIn);
    6868    }
    69 
     69    else if(strcmp(&buffer[4], "RUNH") == 0)
     70    {
     71        fileIn->seekg(4, ios::cur);
     72        delete [] buffer;
     73        return new MCorsikaFormatRaw(fileIn, true);
     74    }
    7075    if (*reinterpret_cast<unsigned int*>(buffer) == kSyncMarker)
    7176    {
     77        //std::cout << "EventIO!!" << std::endl;
    7278        delete [] buffer;
    7379        return new MCorsikaFormatEventIO(fileIn);
    … …  
    8490Bool_t MCorsikaFormat::Read(void *ptr, Int_t i) const
    8591{
     92        //std::cout << "RawRead of " << i << " Bytes" << std::endl;
    8693   fIn->read((char*)ptr, i);
    8794   return !fIn->fail();
    … …  
    114121                                    Int_t & blockLength) const
    115122{
    116 
     123        //std::cout << "||||||| NextBlock " << std::endl;
    117124    int blockHeader;
    118125    fIn->read((char*)&blockHeader, 4);
    … …  
    120127        return kFALSE;
    121128
     129
    122130    blockVersion = 0;
    123131    blockIdentifier = 0;
    … …  
    134142      case 1213486661 : // EVTH
    135143         if (readState != 10)
    136             blockType = 1202;
     144         { blockType = 1202; }
    137145         else
    138146            {
    139             blockType = 1105;
     147                blockType = 1105;
    140148            fIn->seekg(-4, ios::cur);
    141149            blockLength += 4;
    … …  
    159167    for (int i=1; i<22; i++)
    160168    {
    161         fIn->seekg(-i*273*4, ios::end);
     169        if(fFortranRaw)
     170                fIn->seekg(-i*273*4 - 4, ios::end);
     171        else
     172                fIn->seekg(-i*273*4, ios::end);
    162173
    163174        char runh[5]="\0\0\0\0";
    … …  
    196207      fIn->read((char*)blockHeader, 3 * sizeof(int));
    197208
    198       if (fIn->eof())
     209     if (fIn->eof())
    199210         return kFALSE;
    200211
    … …  
    204215      blockIdentifier = blockHeader[1];
    205216      blockLength     = blockHeader[2] & 0x3FFFFFFF;
     217
     218      std::cout << "ReadState " << readState << " | Type:"
     219                << blockType << " Version:" << blockVersion << " Length:" << blockLength << " Id:" << blockIdentifier << std::endl;
     220
    206221      }
    207222   else
    … …  
    210225       fIn->read((char*)blockHeader, 4 * sizeof(int));
    211226
    212        if (fIn->eof())
     227      if (fIn->eof())
    213228           return kFALSE;
    214229
    … …  
    218233       blockIdentifier = blockHeader[2];
    219234       blockLength     = blockHeader[3] & 0x3FFFFFFF;
     235
     236       std::cout << "ReadState " << readState << " | Type:"
     237                 << blockType << " Version:" << blockVersion << " Length:" << blockLength << " Id:" << blockIdentifier << std::endl;
     238
    220239
    221240       if (blockType == 1200  || blockType == 1210 ||

  • branches/Corsika7405Compatibility/mcorsika/MCorsikaFormat.h

    r10213 r18179  
    1717protected:
    1818   std::istream  *fIn;
     19
     20   Bool_t fMMCS = true;
    1921
    2022public:
    … …  
    3840           Bool_t Read(void *ptr, Int_t i) const;
    3941
     42
    4043   static MCorsikaFormat *CorsikaFormatFactory(const char *fileName);
     44
     45   const Bool_t isMMCS() {return fMMCS;}
    4146};
    4247
    … …  
    4550{
    4651private:
     52        Bool_t fFortranRaw = false;
    4753
    4854public:
    4955   MCorsikaFormatRaw(std::istream * in)
    5056        : MCorsikaFormat(in) {}
     57   MCorsikaFormatRaw(std::istream* in, Bool_t fortranRaw)
     58            : MCorsikaFormat(in), fFortranRaw (fortranRaw) {}
    5159
    5260   Bool_t NextBlock(Int_t   readState, Int_t & blockType, Int_t & blockVersion,

  • branches/Corsika7405Compatibility/mcorsika/MCorsikaRead.cc

    r14892 r18179  
    290290                  status = fInFormat->Read(buffer, 272 * sizeof(Float_t));
    291291                  status = fRunHeader->ReadEventHeader(buffer);
     292                  //std::cout << "\x1b[31;1m EVTH Header: \x1b[0m";
     293                  //fRunHeader->Print(nullptr);
    292294                  break;
    293295                  }
    … …  
    329331Int_t MCorsikaRead::ReadTelescopePosition()
    330332{
    331    if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
     333        if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
    332334
    333335   if (fTelescopeIdx>=fNumTelescopes)
    … …  
    474476Int_t MCorsikaRead::Process()
    475477{
     478        //std::cout << "Process \n\n" << std::endl;
    476479   while (1)  // loop for multiple input files
    477480   {
     481           //std::cout << "OuterLoop \n\n" << std::endl;
     482
    478483      if (fInFormat)
    479484      {
    … …  
    492497
    493498         Int_t status = kTRUE;
     499
     500         //std::cout << "||||ReadState: " << fReadState
     501        //               << "  Block: " << fBlockType << "\n\n" << std::endl;
     502
    494503         switch (fBlockType)
    495504            {
    … …  
    500509
    501510            case 1201:       // telescope position
     511            //  std::cout << "TelescopPos ||||||" << std::endl;
    502512               status = ReadTelescopePosition();
    503513               break;
    … …  
    505515            case 1202:  // the event header
    506516               Float_t buffer[272];
     517               static int asd = 0;
     518              // std::cout << "EventHeader |||||| " << ++asd << std::endl;
    507519               if (!fInFormat->Read(buffer, 272 * sizeof(Float_t)))
    508520                  return kFALSE;
    … …  
    511523                  {
    512524                  fRunHeader->ReadEventHeader(buffer);
    513                   fRunHeader->Print();
     525                  //fRunHeader->Print();
    514526                  }
    515527
    516528               status = fEvtHeader->ReadEvt(buffer);
     529
     530               //fEvtHeader->Print();
     531               //std::cout << "TotReuse: " << fEvtHeader->GetTotReuse();
     532
    517533               if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
    518534                  {
    … …  
    527543
    528544            case 1204: // top level block for one array (only for eventio data)
     545            //  std::cout << "EventIO1 |||||||" << std::endl;
    529546               if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
    530547                  {
    … …  
    535552                  // skip this array of telescopes
    536553                  fInFormat->Seek(fBlockLength);
    537 
     554                   std::cout << "Skip Array!" << std::endl;
    538555               break;
    539556
    … …  
    541558            case 1205: // eventio data
    542559               {
     560              // std::cout << "EventIO2 |||||||" << std::endl;
    543561               Int_t telIdx   = fBlockIdentifier % 1000;
    544                if (fBlockVersion == 0                               &&
    545                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
    546                   {
     562           //    std::cout << "TelIDX: " << telIdx << " fTelescopeIdx:" << fTelescopeIdx
     563            //             << " fBlockVersion" << fBlockVersion << std::endl;
     564
     565               if (fBlockVersion == 0  && (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx) )
     566               {
    547567                  status = fEvent->ReadEventIoEvt(fInFormat);
     568
     569                  //std::cout << " ArrayIdx:" << fBlockIdentifier / 1000.0 << std::endl;
    548570
    549571                  Int_t arrayIdx = fBlockIdentifier / 1000;
    … …  
    560582               else
    561583                  // skip this telescope
     584               {
    562585                  fInFormat->Seek(fBlockLength);
     586
     587                   std::cout << "Skip this telescope!" << std::endl;
     588               }
    563589               }
    564590               break;
    565591
    566592            case 1209:  // the event end
     593           //   std::cout << "EventEnd ||||||" << std::endl;
    567594               status = fEvtHeader->ReadEvtEnd(fInFormat);
    568595               
    … …  
    581608
    582609            case 1210:  // the run end
     610                //std::cout << "RunEnd ||||||" << std::endl;
    583611               status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
    584612               fNumEvents += fRunHeader->GetNumEvents();
    … …  
    594622
    595623            case 1105:  // event block of raw format
     624                static int NmrEventBlock = 0;
     625                //std::cout << "EventBlock ||||||" << ++NmrEventBlock << std::endl;
    596626               if (fReadState == 2 || fReadState == 10)
    597627                  {
    … …  
    607637
    608638            case 1109:  // save corsika events
     639                //std::cout << "SaveCorsika ||||||" << std::endl;
    609640               fEvtHeader->InitXY();
    610641               status = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0],
    … …  
    630661            default:
    631662               // unknown block, we skip it
     663               std::cout << "UnknownBlock " << fBlockType << std::endl;
    632664               fReadState = 15;
    633665               fInFormat->Seek(fBlockLength);
    634666
    635667            }
     668
     669         std::cout << "Status: " << status << std::endl;
    636670
    637671         if (status != kTRUE)
    … …  
    640674
    641675         Int_t headerStatus =  ReadNextBlockHeader();
     676
     677         std::cout << "Header Status: " << headerStatus << std::endl;
     678
     679
    642680         if (headerStatus == kFALSE)
    643681            // end of file

  • branches/Corsika7405Compatibility/mcorsika/MCorsikaRunHeader.cc

    r10060 r18179  
    179179    fCerenkovFlag = TMath::Nint(g[75]);
    180180
    181     fZdMin = g[79];                // lower edge of theta in °
    182     fZdMax = g[80];                // upper edge of theta in °
    183     fAzMin = 180-g[81];            // lower edge of phi   in °
    184     fAzMax = 180-g[82];            // upper edge of phi   in °
     181    fZdMin = g[79];                // lower edge of theta in �
     182    fZdMax = g[80];                // upper edge of theta in �
     183    fAzMin = 180-g[81];            // lower edge of phi   in �
     184    fAzMax = 180-g[82];            // upper edge of phi   in �
    185185    // FIXME: Correct for direction of magnetic field!
    186186
    … …  
    207207    fWavelengthMax = g[95];        // Cherenkov bandwidth upper end in nm
    208208
    209     fViewConeInnerAngle = g[151];  // inner angle of view cone (°)
    210     fViewConeOuterAngle = g[152];  // outer angle of view cone (°)
     209    std::cout << "Version: " << fProgramVersion << std::endl;
     210    std::cout << "Wavelength: " << fWavelengthMin << " - " << fWavelengthMax << std::endl;
     211    if((abs(fWavelengthMax - 700)< 0.1) && (abs(fProgramVersion - 7.4005)< 1E-5)) fWavelengthMax = 900;
     212
     213    fViewConeInnerAngle = g[151];  // inner angle of view cone (�)
     214    fViewConeOuterAngle = g[152];  // outer angle of view cone (�)
    211215
    212216    return kTRUE;

  • branches/Corsika7405Compatibility/msim/MPhotonData.cc

    r10060 r18179  
    273273Int_t MPhotonData::FillEventIO(Float_t f[8])
    274274{
    275     // photons in this bunch
     275        // photons in this bunch
    276276    const UInt_t n = TMath::Nint(f[6]);
    277277    if (n==0)
    … …  
    292292    fTag        = -1;
    293293    fWeight     =  1;
     294
    294295
    295296    return n;

  • branches/Corsika7405Compatibility/msim/MPhotonData.h

    r9941 r18179  
    1616#ifndef ROOT_TQuaternion
    1717#include <math.h>
    18 #define sqrt ::sqrt
     18//#define sqrt ::sqrt
    1919#include <TQuaternion.h>
    2020#undef sqrt

  • branches/Corsika7405Compatibility/msim/MPhotonEvent.cc

    r10167 r18179  
    512512   fInFormat->Read(bunchHeader, 3 * sizeof(Int_t));
    513513
     514
     515   //std::cout << "BunchHeader: " << bunchHeader[0] << " "
     516        //         << bunchHeader[1] << " "
     517        //         << bunchHeader[2] << std::endl;
     518
    514519   Int_t n = 0;
    515520
    … …  
    520525
    521526      if (Add(n).FillEventIO(buffer))
    522          n++;
     527          n++;
    523528      }
    524529

  • branches/Corsika7405Compatibility/msimcamera/MSimAPD.cc

    r17643 r18179  
    307307            *fLog << err << "ERROR - MSimAPD: Weight of " << i << "-th photon not 1, but " << ph.GetWeight() << endl;
    308308            ph.Print();
    309             return kERROR;
     309            //return kERROR;
    310310        }
    311311
Note: See TracChangeset for help on using the changeset viewer.