Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 1075) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 1078) @@ -21,7 +21,7 @@ // // $RCSfile: camera.cxx,v $ -// $Revision: 1.28 $ +// $Revision: 1.29 $ // $Author: blanch $ -// $Date: 2001-10-26 16:31:45 $ +// $Date: 2001-11-14 17:38:23 $ // //////////////////////////////////////////////////////////////////////// @@ -69,4 +69,5 @@ #include "MMcEvt.hxx" #include "MMcTrig.hxx" +#include "MMcRunHeader.hxx" #include "MMcTrigHeader.hxx" #include "MMcFadcHeader.hxx" @@ -431,5 +432,5 @@ int simulateNSB; //@< Will we simulate NSB? - int nphe2NSB; //@< From how many ohe will we simulate NSB? + int nphe2NSB; //@< From how many phe will we simulate NSB? float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel) float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel @@ -455,4 +456,5 @@ EXTWAVEBAND5 }; + float zenfactor=1.0; // correction factor calculated from the extinction float qTailCut; //@< Tail Cut value @@ -473,5 +475,5 @@ float fpixelthres[TRIGGER_PIXELS]; //@< Threshold values - TArrayC *fadcValues; //@< the analog Fadc siganl for pixels + TArrayC *fadcValues; //@< the analog Fadc signal for pixels float plateScale_cm2deg; //@< plate scale (deg/cm) @@ -479,4 +481,22 @@ int still_in_loop = FALSE; + + //@< Variables to fill the McRunHeader + Int_t sfRaH = 5; + Int_t sfRaM = 34; + Int_t sfRaS = 32; + Int_t sfDeD = 22; + Int_t sfDeM = 00; + Int_t sfDeS = 52; + Float_t shthetamax = 0.0; + Float_t shthetamin = 0.0; + Float_t shphimax = 0.0; + Float_t shphimin = 0.0; + Float_t telestheta = 0.0; + Float_t telesphi = 0.0; + Float_t sofftheta = 0.0; + Float_t soffphi = 0.0; + UInt_t corsika = 5200 ; + struct camera cam; // structure holding the camera definition @@ -743,7 +763,7 @@ // Initialise McTrig information class if we want to save trigger informtion - MMcTrig **McTrig; - MMcTrigHeader **HeaderTrig; - MMcFadcHeader **HeaderFadc; + MMcTrig **McTrig = NULL; + MMcTrigHeader **HeaderTrig = NULL; + MMcFadcHeader **HeaderFadc = NULL; if (Write_McTrig){ @@ -796,12 +816,10 @@ // Header Tree - MRawRunHeader **RunHeader; - - RunHeader = new MRawRunHeader * [1]; - RunHeader[0] = new MRawRunHeader(); + MRawRunHeader *RunHeader= new MRawRunHeader(); + MMcRunHeader *McRunHeader = new MMcRunHeader(); // Header branch - MRawEvtHeader **EvtHeader; + MRawEvtHeader **EvtHeader = NULL; if (Write_RawEvt) { @@ -815,5 +833,5 @@ // Data branch - MRawEvtData **EvtData; // Data branch + MRawEvtData **EvtData = NULL; // Data branch if (Write_RawEvt) { @@ -822,4 +840,6 @@ for (i=0;iInit(RunHeader); // We need thr RunHeader to read + // number of pixels } } @@ -833,6 +853,4 @@ TFile outfile_temp ( rootname , "RECREATE" ); - Int_t bsize=128000; Int_t split=1; - // create a Tree for the Header Event TTree HeaderTree("RunHeaders","Headers of Run"); @@ -843,10 +861,13 @@ HeaderTree.Branch("MRawRunHeader","MRawRunHeader", - &RunHeader[0], bsize, split); + &RunHeader); + + HeaderTree.Branch("MMcRunHeader","MMcRunHeader", + &McRunHeader); if(!Trigger_Loop && Write_McTrig){ HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", - &HeaderTrig[0], bsize, split); + &HeaderTrig[0]); } if (Trigger_Loop && Write_McTrig){ @@ -858,5 +879,5 @@ strcat (branchname, "."); HeaderTree.Branch(branchname,"MMcTrigHeader", - &HeaderTrig[i], bsize, split); + &HeaderTrig[i]); } } @@ -865,5 +886,5 @@ HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", - &HeaderFadc[0], bsize, split); + &HeaderFadc[0]); } if (Trigger_Loop && Write_McFADC){ @@ -875,5 +896,5 @@ strcat (branchname, "."); HeaderTree.Branch(branchname,"MMcFadcHeader", - &HeaderFadc[i], bsize, split); + &HeaderFadc[i]); } } @@ -881,10 +902,12 @@ // Fill branches for MRawRunHeader - RunHeader[0]->SetMagicNumber(kMagicNumber); - RunHeader[0]->SetFormatVersion(2); - RunHeader[0]->SetSoftVersion((UShort_t) (VERSION*10)); - RunHeader[0]->SetRunType(256); - RunHeader[0]->SetRunNumber(0); - RunHeader[0]->SetNumSamples(0,FADC_SLICES); + RunHeader->SetMagicNumber(kMagicNumber); + RunHeader->SetFormatVersion(2); + RunHeader->SetSoftVersion((UShort_t) (VERSION*10)); + RunHeader->SetRunType(256); + RunHeader->SetRunNumber(0); + RunHeader->SetNumSamples(0,FADC_SLICES); + RunHeader->SetNumCrates(1); + RunHeader->SetNumPixInCrate(CAMERA_PIXELS); @@ -956,5 +979,5 @@ // Fill the Header Tree with the current leaves of each branch - HeaderTree.Fill() ; + // HeaderTree.Fill() ; @@ -965,5 +988,5 @@ EvtTree.Branch("MMcEvt","MMcEvt", - &McEvt, bsize, split); + &McEvt); } @@ -972,11 +995,11 @@ if (Write_RawEvt){ EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", - &EvtHeader[0], bsize, split); + &EvtHeader[0]); EvtTree.Branch("MRawEvtData","MRawEvtData", - &EvtData[0], bsize, split); + &EvtData[0]); } if (Write_McTrig){ EvtTree.Branch("MMcTrig","MMcTrig", - &McTrig[0], bsize, split); + &McTrig[0]); } } @@ -990,5 +1013,5 @@ strcat (branchname, "."); EvtTree.Branch(branchname,"MMcTrig", - &McTrig[i], bsize, split); + &McTrig[i]); } } @@ -1003,5 +1026,5 @@ strcat (branchname, "."); EvtTree.Branch(branchname,"MRawEvtHeader", - &EvtHeader[i], bsize, split); + &EvtHeader[i]); } for(char branchname[15],i=0;iIsBatch()) { - fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); - // return 1; - } - if(FADC_Scan){ - TApplication theAppFadc("App", &argc, argv); - - if (gROOT->IsBatch()) { - fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); - // return 1; - } + if (gROOT->IsBatch()) { + fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); + // return 1; + } + } + if(FADC_Scan){ + //TApplication theAppFadc("App", &argc, argv); + + if (gROOT->IsBatch()) { + fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); + // return 1; + } } @@ -1060,5 +1082,5 @@ k = produce_nsbrates( starfieldname, &cam, - nsbrate_phepns ); + nsbrate_phepns); if (k != 0){ @@ -1075,8 +1097,13 @@ // calculate nsb rate including diffuse and starlight - for(i=0; iSetFirstLevel (Lev1); @@ -1700,4 +1734,44 @@ } // end if else found start of run } // end big while loop + + //<@ Finally we should fill th McRunHeader + + + get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS); + get_teles_axis(&telestheta,&telesphi); + get_source_off(&sofftheta,&soffphi); + mcevth.get_theta_range(&shthetamin, &shthetamax); + mcevth.get_phi_range(&shphimin,&shphimax); + corsika=get_corsika_ver(); + if(!Trigger_Loop) icontrigger=0; + McRunHeader->Fill(icontrigger, + !Write_All_Images, + Write_McEvt, + Write_McTrig, + Write_McFADC, + CAMERA_PIXELS, + (UInt_t)ntshow, + (UInt_t)ntrigger, + sfRaH, + sfRaM, + sfRaS, + sfDeD, + sfDeM, + sfDeS, + meanNSB, + telestheta, + telesphi, + sofftheta, + soffphi, + shthetamax, + shthetamin, + shphimax, + shphimin, + corsika, + (UInt_t)(REFL_VERSION*100), + (UInt_t)(VERSION*100)); + + // Fill the Header Tree with the current leaves of each branch + HeaderTree.Fill() ; //++ @@ -2621,11 +2695,11 @@ /******** bpoint_is_in_pix() ***************************************/ - +/* int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){ - /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */ - /* the pixel is assumed to be a "closed set" */ - - double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */ - double c, xx, yy; /* auxiliary variable */ + // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) + // the pixel is assumed to be a "closed set" + + double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel + double c, xx, yy; // auxiliary variable b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum]; @@ -2643,9 +2717,50 @@ if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) || ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){ - return(TRUE); /* point is inside */ + return(TRUE); // point is inside } else{ - return(FALSE); /* point is outside */ - } + return(FALSE); // point is outside + } +} +*/ + +#define sqrt13 0.577350269 // = 1./sqrt(3.) +#define sqrt3 1.732050807 // = sqrt(3.) + +int bpoint_is_in_pix(double dx, double dy, struct camera *pcam) +{ + /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */ + /* the pixel is assumed to be a "closed set" */ + + /* + a = length of one of the edges of one pixel, + b = half the width of one pixel + */ + + const int num = pcam->inumpixels; + const double diam = pcam->dpixdiameter_cm; + + const double diam2 = diam/2; + const double diam3 = diam/sqrt3; + + for (int i=0; idpixsizefactor[i]; + + const double b = diam2 * size; + const double a = diam3 * size; + + const double xx = dx - pcam->dxc[i]; + const double yy = dy - pcam->dyc[i]; + + if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) || + ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){ + + return i; // inside i + } + + // return -1; // outside + } + return -1; // outside } @@ -2789,9 +2904,13 @@ // read the photons data - fread ( flag, SIZE_OF_FLAGS, 1, sp ); // loop over the photons - while ( !isA( flag, FLAG_END_OF_EVENT ) ) { + while (1) { + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + if (isA( flag, FLAG_END_OF_EVENT )) + break; memcpy( (char*)&photon, flag, SIZE_OF_FLAGS ); @@ -2808,11 +2927,5 @@ if (t-*tmin_ns>TOTAL_TRIGGER_TIME) { - - // read next Photon - - fread( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop, the photon is lost - continue; + continue; } @@ -2833,45 +2946,16 @@ if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){ - - // cout << " lost first check\n"; - - // read next CPhoton - - fread ( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop, the photon is lost continue; } - ipixnum = -1; - - for(i=0; iinumpixels; i++){ - if( bpoint_is_in_pix( cx, cy, i, cam) ){ - ipixnum = i; - break; - } - } - - if(ipixnum==-1){// the photon is in none of the pixels - - // cout << " lost pixlization\n"; - - // read next CPhoton - - fread ( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop, the photon is lost - continue; - } - - if(ipixnum==0) {// the phton is in the central pixel, which is not used for trigger - // read next CPhoton - - fread ( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop, the photon is lost - continue; - } + ipixnum = bpoint_is_in_pix(cx, cy, cam); + + // -1 = the photon is in none of the pixels + // 0 = the phton is in the central pixel, which is not used for trigger + if (ipixnum==-1 || ipixnum==0) { + continue; + } + //+++ // QE simulation @@ -2885,12 +2969,4 @@ if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){ - - // cout << " lost\n"; - - // read next Photon - - fread ( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop continue; @@ -2911,14 +2987,5 @@ if ( (RandomNumber) > qe ) { - - // cout << " lost\n"; - - // read next Photon - - fread ( flag, SIZE_OF_FLAGS, 1, sp ); - - // go to beginning of loop continue; - } @@ -2938,10 +3005,5 @@ *itotnphe += 1; - - // read next Photon - - fread( flag, SIZE_OF_FLAGS, 1, sp ); - - } // end while, i.e. found end of event + } return(0); @@ -3248,4 +3310,7 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.28 2001/10/26 16:31:45 blanch +// Removing several warnings. +// // Revision 1.27 2001/09/05 10:04:33 blanch // *** empty log message *** @@ -3317,4 +3382,7 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.28 2001/10/26 16:31:45 blanch +// Removing several warnings. +// // Revision 1.27 2001/09/05 10:04:33 blanch // *** empty log message ***