Changeset 9937
- Timestamp:
- 09/22/10 15:25:55 (14 years ago)
- Location:
- trunk/Mars
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/Changelog
r9936 r9937 42 42 * msim/MHPhotonEvent.cc: 43 43 - implemented types 6-8 44 45 * mcorsika/MCorsikaEvtHeader.[h,cc]: 46 - added fNumReuse 47 - added member functions to increase and reset the number of reusages 48 - increased class version number accordingly 49 50 * mcorsika/MCorsikaRead.[h,cc]: 51 - count the number of showers (events) times its reusage 52 for fNumTotalEvents 53 - adapted ReadEvent to re-usage of showers (still needs to be tested for 54 EventIO) 55 - Fixed return type (should be Int_t) of ReadEvent 56 - take the number of reusages in PostProcess into account 57 58 * mcorsika/MCorsikaRunHeader.[h,cc]: 59 - added total number of reusages fNumReuse 60 - increased class version number accordingly 61 62 * mmc/MMcEvt.hxx: 63 - added setter for the reuse counter 64 65 * msim/MPhotonData.[h,cc]: 66 - removed fNumPhotons (we allow 1-photon bunches only anyway and MMCS 67 doesn't even distribute this number) 68 - increased class version number accordingly 69 70 * msim/MSimMMCS.cc: 71 - correctly propagate the number of events (number of showers times reusage) 72 - propagate counter for reusage to MMcEvt 73 74 * msim/MPhotonEvent.[h,cc]: 75 - added functions to calculate mean x and mean y 76 - propagate re-usage counter through ReadEvent 44 77 45 78 -
trunk/Mars/NEWS
r9878 r9937 67 67 mapping of the regular patches, resmc/fact-trigger-all.txt the mapping 68 68 of all patches. 69 70 * Added reading of re-used corsika showers (only supported if MMCS is used) 69 71 70 72 ;star: -
trunk/Mars/mcorsika/MCorsikaEvtHeader.cc
r9616 r9937 32 32 // 33 33 // 34 // Class Version 3: 35 // ---------------- 36 // - UInt_t fNumReuse 37 // 38 // 34 39 ///////////////////////////////////////////////////////////////////////////// 35 40 #include "MCorsikaEvtHeader.h" … … 55 60 // 56 61 MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title) 57 : f X(0), fY(0)62 : fNumReuse((UInt_t)-1), fX(0), fY(0) 58 63 { 59 64 fName = name ? name : "MCorsikaEvtHeader"; … … 104 109 { 105 110 *fLog << all; 106 *fLog << "Event Number: " << dec << fEvtNumber << endl;111 *fLog << "Event Number: " << dec << fEvtNumber << " (reused=" << fNumReuse << ")" << endl; 107 112 // *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; 108 113 *fLog << "Energy: " << fTotalEnergy << "GeV" << endl; … … 122 127 // return FALSE if there is no header anymore, else TRUE 123 128 // 124 Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat * 129 Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat *fInFormat) 125 130 { 126 127 131 if (!fInFormat->SeekNextBlock("EVTH", 1202)) 128 132 return kFALSE; … … 151 155 fMomentumZ = -f[8]; // Does this minus make sense?! 152 156 153 // FIXME: Correct for direction of magnetic field!154 157 fZd = f[9]; 155 158 fAz = TMath::Pi()-f[10]; 156 159 157 const Int_t n = TMath::Nint(f[96]); 158 if (n!=1) 159 { 160 *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; 161 *fLog << warn << "WARNING - This error is replaced by a warning." << endl; 162 } 160 // f[96] // Number of reuse of event [max=20] 163 161 164 fX = f[117]; //fX = f[97]; 165 fY = -f[97]; //fY = f[117]; 162 // FIXME: Check fNumReuse<20 163 fX = f[117 + fNumReuse]; 164 fY = -f[ 97 + fNumReuse]; 165 166 fWeightedNumPhotons = 0; 166 167 167 168 return !fInFormat->Eof(); … … 171 172 // this member function is for reading the event end block 172 173 173 Bool_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)174 Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) 174 175 { 175 176 176 if (!fInFormat->SeekNextBlock("EVTE", 1209)) 177 177 return kFALSE; … … 192 192 } 193 193 194 // FIXME: What's the meaning of this if in case of reusing the event this number 195 // exists only once? 194 196 fWeightedNumPhotons = f[1]; 195 197 196 198 return !fInFormat->Eof(); 197 199 } 198 -
trunk/Mars/mcorsika/MCorsikaEvtHeader.h
r9616 r9937 20 20 private: 21 21 UInt_t fEvtNumber; // Event number 22 UInt_t fNumReuse; // Counter of the reuse of the same shower 22 23 // UInt_t fParticleID; // Particle ID (see MMcEvtBasic or CORSIKA manual) 23 24 Float_t fTotalEnergy; // [GeV] … … 46 47 47 48 UInt_t GetEvtNumber() const { return fEvtNumber; } 49 UInt_t GetNumReuse() const { return fNumReuse; } 48 50 // UInt_t GetParticleID() const { return fParticleID; } 49 51 … … 62 64 Double_t GetImpact() const; 63 65 64 Int_t ReadEvt(MCorsikaFormat * fInFormat); // read in event header block65 Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); // read in event end block66 void IncNumReuse() { fNumReuse++; } 67 void ResetNumReuse() { fNumReuse=(UInt_t)-1; } 66 68 67 ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER 69 Int_t ReadEvt(MCorsikaFormat *informat); // read in event header block 70 Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block 71 72 ClassDef(MCorsikaEvtHeader, 3) // Parameter Conatiner for raw EVENT HEADER 68 73 }; 69 74 -
trunk/Mars/mcorsika/MCorsikaRead.cc
r9616 r9937 319 319 } 320 320 321 fNumTotalEvents += fRunHeader->GetNumEvents() ;321 fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse(); 322 322 continue; 323 323 } … … 393 393 // Read a single event from the stream 394 394 // 395 Bool_t MCorsikaRead::ReadEvent() 396 { 397 // 398 // Read in the next EVENT HEADER (see specification), 399 // if there is no next event anymore stop eventloop 400 // 401 Int_t rc = fEvtHeader->ReadEvt(fInFormat); //read event header block 402 if (!rc) 403 return kFALSE; 404 405 rc = fEvent->ReadCorsikaEvt(fInFormat); 406 407 /* 408 // Check if we are allowed to stop reading in the middle of the data 409 if (rc==kFALSE && !fForce) 410 { 411 *fLog << err; 412 *fLog << "ERROR - End of file in the middle of the data stream!" << endl; 413 *fLog << " Set the force-option to force processing." << endl; 414 return kERROR; 415 } 416 */ 417 418 return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fInFormat) : rc; 395 Int_t MCorsikaRead::ReadEvent() 396 { 397 // Store the position to re-read a single event several times 398 // FIXME: Does this work with EventIO, too? 399 if (fEvtHeader->GetNumReuse()>=fRunHeader->GetNumReuse()-1) 400 { 401 fInFormat->StorePos(); 402 fEvtHeader->ResetNumReuse(); 403 } 404 else 405 fInFormat->ResetPos(); 406 407 fEvtHeader->IncNumReuse(); 408 409 // Read the event header 410 const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); 411 if (rc1!=kTRUE) 412 return rc1; 413 414 // Check if reusage number should be reset and increase it 415 // Note, that the trick here is that after reset it is set to -1 416 417 // Read the photons corresponding to the i-th core location 418 const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1); 419 if (rc2!=kTRUE) 420 return rc2; 421 422 // read event end 423 return fEvtHeader->ReadEvtEnd(fInFormat); 419 424 } 420 425 … … 434 439 { 435 440 // Read a single event from file 436 const Bool_t rc = ReadEvent();441 const Int_t rc = ReadEvent(); 437 442 438 443 // kFALSE means: end of file (try next one) … … 465 470 Int_t MCorsikaRead::PostProcess() 466 471 { 472 const UInt_t n = fNumEvents*fRunHeader->GetNumReuse(); 473 467 474 // 468 475 // Sanity check for the number of events 469 476 // 470 if ( fNumEvents==GetNumExecutions()-1 || GetNumExecutions()==0)477 if (n==GetNumExecutions()-1 || GetNumExecutions()==0) 471 478 return kTRUE; 472 479 … … 474 481 *fLog << "Warning - number of read events (" << GetNumExecutions()-1; 475 482 *fLog << ") doesn't match number in run header(s) ("; 476 *fLog << f NumEvents<< ")." << endl;483 *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl; 477 484 478 485 return kTRUE; -
trunk/Mars/mcorsika/MCorsikaRead.h
r9616 r9937 29 29 UInt_t fNumTotalEvents; //! total number of events in all files 30 30 31 ifstream *fIn;//! input stream (file to read from)32 MCorsikaFormat *fInFormat; //! access to input corsika data31 ifstream *fIn; //! input stream (file to read from) 32 MCorsikaFormat *fInFormat; //! access to input corsika data 33 33 34 34 MParList *fParList; //! tasklist to call ReInit from … … 42 42 Int_t OpenNextFile(Bool_t print=kTRUE); 43 43 Bool_t CalcNumTotalEvents(); 44 Bool_tReadEvent();44 Int_t ReadEvent(); 45 45 46 46 Int_t PreProcess(MParList *pList); -
trunk/Mars/mcorsika/MCorsikaRunHeader.cc
r9616 r9937 42 42 // + Float_t fAtmosphericCoeffC[5] 43 43 // + UInt_t fCerenkovFlag 44 // 45 // Class Version 3: 46 // ---------------- 47 // + UInt_t fNumReuse 44 48 // 45 49 //////////////////////////////////////////////////////////////////////////// … … 164 168 fInFormat->UnreadLastHeader(); 165 169 166 const Int_t n= TMath::Nint(g[96]); // Number i of uses of each cherenkov event167 if ( n!=1)170 fNumReuse = TMath::Nint(g[96]); // Number i of uses of each cherenkov event 171 if (fNumReuse!=1) 168 172 { 169 173 *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; … … 269 273 *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; 270 274 if (fNumEvents>0) 271 *fLog << "Num Events: " << fNumEvents << endl;275 *fLog << "Num Events: " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl; 272 276 *fLog << "Obs Level: "; 273 277 for (Byte_t i=0; i<fNumObsLevel; i++) -
trunk/Mars/mcorsika/MCorsikaRunHeader.h
r9616 r9937 33 33 34 34 UInt_t fRunNumber; // Run number 35 UInt_t fNumReuse; // Number of how many times the same shower is used 35 36 UInt_t fParticleID; // Particle ID (see MMcEvtBasic or CORSIKA manual) 36 37 UInt_t fNumEvents; // Number of events … … 76 77 UInt_t GetParticleID() const { return fParticleID; } 77 78 UInt_t GetNumEvents() const { return fNumEvents; } 79 UInt_t GetNumReuse() const { return fNumReuse; } 78 80 79 81 const MTime &GetRunStart() const { return fRunStart; } … … 128 130 void Print(Option_t *t=NULL) const; 129 131 130 ClassDef(MCorsikaRunHeader, 2) // storage container for general info132 ClassDef(MCorsikaRunHeader, 3) // storage container for general info 131 133 }; 132 134 #endif -
trunk/Mars/mmc/MMcEvt.hxx
r9595 r9937 103 103 104 104 void SetEvtNumber(UInt_t n) { fEvtNumber=n; } 105 void SetEventReuse(UInt_t n) { fEventReuse=n; } 105 106 void SetPhotElfromShower(UInt_t n) { fPhotElfromShower=n; } 106 107 -
trunk/Mars/msim/MPhotonData.cc
r9820 r9937 16 16 ! 17 17 ! 18 ! Author(s): Thomas Bretz, 12/2000 <mailto:t bretz@astro.uni-wuerzburg.de>18 ! Author(s): Thomas Bretz, 12/2000 <mailto:thomas.bretz@epfl.ch> 19 19 ! Author(s): Qi Zhe, 06/2007 <mailto:qizhe@astro.uni-wuerzburg.de> 20 20 ! 21 ! Copyright: CheObs Software Development, 2000-20 0921 ! Copyright: CheObs Software Development, 2000-2010 22 22 ! 23 23 ! … … 34 34 // Version 1: 35 35 // ---------- 36 // - First implementation 36 // * First implementation 37 // 38 // Version 2: 39 // ---------- 40 // - fNumPhotons 37 41 // 38 42 ///////////////////////////////////////////////////////////////////////////// … … 57 61 MPhotonData::MPhotonData(/*const char *name, const char *title*/) 58 62 : fPosX(0), fPosY(0), fCosU(0), fCosV(0), fTime(0), fWavelength(0), 59 fNumPhotons(1),fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED),63 /*fNumPhotons(1),*/ fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED), 60 64 fTag(-1), fWeight(1) 61 65 { … … 82 86 MPhotonData &d = static_cast<MPhotonData&>(obj); 83 87 84 d.fNumPhotons = fNumPhotons;88 // d.fNumPhotons = fNumPhotons; 85 89 d.fPosX = fPosX; 86 90 d.fPosY = fPosY; … … 185 189 186 190 fPrimary = MMcEvtBasic::kUNDEFINED; 187 fNumPhotons = 1;191 // fNumPhotons = 1; 188 192 fTag = -1; 189 193 fWeight = 1; … … 201 205 // system intpo our own. 202 206 // 203 Int_t MPhotonData::FillCorsika(Float_t f[7] )207 Int_t MPhotonData::FillCorsika(Float_t f[7], Int_t i) 204 208 { 205 209 const UInt_t n = TMath::Nint(f[0]); 206 210 if (n==0) 211 // FIXME: Do we need to decode the rest anyway? 212 return kCONTINUE; 213 214 // Check reuse 215 const Int_t reuse = (n/1000)%100; // Force this to be 1! 216 if (reuse!=i) 207 217 return kCONTINUE; 208 218 … … 210 220 fWavelength = n%1000; 211 221 fPrimary = MMcEvtBasic::ParticleId_t(n/100000); 212 fNumPhotons = (n/1000)%100; // Force this to be 1!213 214 if (fNumPhotons!=1)215 {216 // FIXME: Could be done in MPhotonEvent::ReadCorsikaEvent217 gLog << err << "ERROR - MPhotonData::FillCorsika: fNumPhotons not 1, but " << fNumPhotons << endl;218 gLog << " This is not yet supported." << endl;219 return kERROR;220 }221 222 222 223 // x=north, y=west … … 254 255 Int_t MPhotonData::FillEventIO(Float_t f[8]) 255 256 { 257 if (TMath::Nint(f[6])!=1) 258 { 259 gLog << err << "ERROR - Bunch sizes != 1 are not supported." << endl; 260 return kFALSE; 261 } 262 256 263 fPosX = f[1]; // xpos relative to telescope [cm] 257 264 fPosY = -f[0]; // ypos relative to telescope [cm] … … 260 267 fTime = f[4]; // a relative arival time [ns] 261 268 fProductionHeight = f[5]; // altitude of emission [cm] 262 fNumPhotons = TMath::Nint(f[6]); // photons in this bunch269 // fNumPhotons = TMath::Nint(f[6]); // photons in this bunch 263 270 fWavelength = TMath::Nint(f[7]); // so far always zeor = unspec. [nm] 264 265 271 266 272 // Now reset all data members which are not in the stream … … 271 277 return kTRUE; 272 278 } 273 279 /* 274 280 // -------------------------------------------------------------------------- 275 281 // … … 299 305 return rc==kTRUE ? !fin.eof() : rc; 300 306 } 307 */ 301 308 302 309 // -------------------------------------------------------------------------- … … 308 315 { 309 316 gLog << inf << endl; 310 gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl; 317 // gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl; 318 gLog << "Origin: " << MMcEvtBasic::GetParticleName(fPrimary) << endl; 311 319 gLog << "Wavelength: " << fWavelength << "nm" << endl; 312 320 gLog << "Pos X/Y Cos U/V: " << fPosX << "/" << fPosY << " " << fCosU << "/" << fCosV << endl; -
trunk/Mars/msim/MPhotonData.h
r9616 r9937 37 37 38 38 Float_t fTime; // [ns] Time since first interaction or entrance into atmosphere 39 // 17M40 39 UShort_t fWavelength; // [nm] Wavelength 41 // 19M 42 UInt_t fNumPhotons; // Number of cherenkov photons ins bunch 40 // UInt_t fNumPhotons; // Number of cherenkov photons ins bunch 43 41 Float_t fProductionHeight; // [cm] Height of bunch production 44 42 MMcEvtBasic::ParticleId_t fPrimary; // Type of emitting particle 45 // 22M46 // gzip47 // 25M48 // raw49 // 32M50 43 51 44 Int_t fTag; //! A tag for external use … … 129 122 130 123 // I/O 131 Int_t ReadCorsikaEvt(istream &fin);132 Int_t ReadRflEvt(istream &fin);124 //Int_t ReadCorsikaEvt(istream &fin); 125 //Int_t ReadRflEvt(istream &fin); 133 126 134 Int_t FillCorsika(Float_t f[7] );127 Int_t FillCorsika(Float_t f[7], Int_t i); 135 128 Int_t FillEventIO(Float_t f[7]); 136 129 Int_t FillRfl(Float_t f[8]); 137 130 138 ClassDef(MPhotonData, 1) //Container to store a cherenkov photon bunch from a CORSUKA file131 ClassDef(MPhotonData, 2) //Container to store a cherenkov photon bunch from a CORSUKA file 139 132 }; 140 133 -
trunk/Mars/msim/MPhotonEvent.cc
r9616 r9937 449 449 } 450 450 451 Double_t MPhotonEvent::GetMeanX() const 452 { 453 const UInt_t n = GetNumPhotons(); 454 455 Double_t mean = 0; 456 for (UInt_t i=0; i<n; i++) 457 mean += operator[](i).GetPosX(); 458 459 return mean / n; 460 } 461 462 Double_t MPhotonEvent::GetMeanY() const 463 { 464 const UInt_t n = GetNumPhotons(); 465 466 Double_t mean = 0; 467 for (UInt_t i=0; i<n; i++) 468 mean += operator[](i).GetPosY(); 469 470 return mean / n; 471 } 472 451 473 // -------------------------------------------------------------------------- 452 474 // 453 475 // Read the Event section from the file 454 476 // 455 Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat * fInFormat)477 Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i) 456 478 { 457 479 Int_t n = 0; … … 509 531 { 510 532 511 const Int_t rc = Add(n).FillCorsika(buffer );533 const Int_t rc = Add(n).FillCorsika(buffer, i); 512 534 switch (rc) 513 535 { … … 577 599 } 578 600 579 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin )601 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i) 580 602 { 581 603 Int_t n = 0; … … 635 657 // Get/Add the n-th entry from the array and 636 658 // fill it with the current 7 floats 637 const Int_t rc = Add(n).FillCorsika(ptr );659 const Int_t rc = Add(n).FillCorsika(ptr, i); 638 660 ptr += 7; 639 661 … … 707 729 708 730 // -------------------------------------------------------------------------- 709 / /710 Int_t MPhotonEvent::ReadRflEvt(std::istream &fin )731 /* 732 Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i) 711 733 { 712 734 Int_t n = 0; … … 735 757 // Now we read a single cherenkov bunch 736 758 //const Int_t rc = static_cast<MPhotonData*>(o)->ReadRflEvt(fin); 737 const Int_t rc = Add(n).ReadRflEvt(fin );759 const Int_t rc = Add(n).ReadRflEvt(fin, i); 738 760 739 761 // Evaluate result from reading event … … 754 776 SetReadyToSave(); 755 777 756 // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;778 // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; 757 779 return kTRUE; 758 780 } 759 781 */ 760 782 // -------------------------------------------------------------------------- 761 783 // -
trunk/Mars/msim/MPhotonEvent.h
r9616 r9937 37 37 Double_t GetTimeMedianDev() const; 38 38 39 Double_t GetMeanX() const; 40 Double_t GetMeanY() const; 41 39 42 TClonesArray &GetArray() { return fData; } 40 43 const TClonesArray &GetArray() const { return fData; } … … 53 56 54 57 // I/O 55 Int_t ReadCorsikaEvt(MCorsikaFormat * fInFormat);56 Int_t ReadCorsikaEvt(istream &fin );57 Int_t ReadRflEvt(istream &fin);58 Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i); 59 Int_t ReadCorsikaEvt(istream &fin, Int_t i); 60 //Int_t ReadRflEvt(istream &fin, Int_t i); 58 61 59 62 // TObject -
trunk/Mars/msim/MSimMMCS.cc
r9595 r9937 131 131 132 132 // FIXME: Is there a way to write them as LAST entry in the file? 133 fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents() );133 fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()); 134 134 fMcRunHeader->SetCorsikaVersion(TMath::Nint(fRunHeader->GetProgramVersion()*100)); 135 135 … … 211 211 212 212 fMcEvt->SetEvtNumber(fEvtHeader->GetEvtNumber()); 213 fMcEvt->SetEventReuse(fEvtHeader->GetNumReuse()); 213 214 fMcEvt->SetPhotElfromShower(0); 214 215
Note:
See TracChangeset
for help on using the changeset viewer.