Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 1705) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 1706) @@ -21,7 +21,7 @@ // // $RCSfile: camera.cxx,v $ -// $Revision: 1.50 $ +// $Revision: 1.51 $ // $Author: blanch $ -// $Date: 2003-01-07 16:33:31 $ +// $Date: 2003-01-14 13:40:17 $ // //////////////////////////////////////////////////////////////////////// @@ -443,4 +443,7 @@ int addElecNoise; //@< Will we add ElecNoise? + float riseDiskThres; + float secureDiskThres; + int simulateNSB; //@< Will we simulate NSB? int nphe2NSB; //@< From how many phe will we simulate NSB? @@ -472,5 +475,8 @@ int anaPixels; - + + int nstoredevents = 0; + int flagstoring = 0; + float fCorrection; //@< Factor to apply to pixel values (def. 1.) @@ -609,4 +615,6 @@ Individual_Thres_Pixel = get_indi_thres_pixel(); + get_secure_threhold(&riseDiskThres,&secureDiskThres); + get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm); @@ -811,4 +819,11 @@ Trigger.CheckThreshold(&qThreshold[0]); + // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel + // is active: secureDiskThres*10000+riseDiskThres + + if(riseDiskThres>0.0) + qThreshold[0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres; + else + qThreshold[0]=0.0; // Initialise McTrig information class if we want to save trigger informtion @@ -963,5 +978,5 @@ RunHeader->SetMagicNumber(kMagicNumber); - RunHeader->SetFormatVersion(3); + RunHeader->SetFormatVersion(4); RunHeader->SetSoftVersion((UShort_t) (VERSION*10)); RunHeader->SetRunType(256); @@ -1288,11 +1303,7 @@ if (reflector_file_version<6){ fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile ); - thetashw = mcevth.get_theta(); - phishw = mcevth.get_phi(); } else{ fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile ); - thetashw = mcevth_2.get_theta(); - phishw = mcevth_2.get_phi(); } @@ -1341,48 +1352,78 @@ m1 = sin(thetashw)*sin(phishw); n1 = cos(thetashw); - - // read the deviation of the telescope with respect to the shower - - if (reflector_file_version<6) + + if (reflector_file_version<6){ + + mcevth.get_core(&coreX, &coreY); + + // read the deviation of the telescope with respect to the shower mcevth.get_deviations ( &thetaCT, &phiCT ); - else - mcevth_2.get_deviations ( &thetaCT, &phiCT ); - - if ( (thetaCT == 0.) && (phiCT == 0.) ) { + + if ( (thetaCT == 0.) && (phiCT == 0.) ) { - // CT was looking to the source (both lines are parallel) - // therefore, we calculate the impact parameter as the distance - // between the CT axis and the core position - - if (reflector_file_version<6) - mcevth.get_core(&coreX, &coreY); - else - mcevth_2.get_core(&coreX, &coreY); - - impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. ); + // CT was looking to the source (both lines are parallel) + // therefore, we calculate the impact parameter as the distance + // between the CT axis and the core position + + impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. ); + thetaCT += thetashw; + phiCT += phishw; + + } else { + + // the shower comes off-axis + + // obtain with this the final direction of the CT + + thetaCT += thetashw; + phiCT += phishw; + + // calculate vector for telescope + + l2 = sin(thetaCT)*cos(phiCT); + m2 = sin(thetaCT)*sin(phiCT); + n2 = cos(thetaCT); + + num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY); + den = (SQR(l1*m2 - l2*m1) + + SQR(m1*n2 - m2*n1) + + SQR(n1*l2 - n2*l1)); + den = sqrt(den); + + impactD = fabs(num)/den; + + } + } + else{ + mcevth_2.get_core(&coreX, &coreY); + thetaCT=mcevth_2.get_theta_CT(); + phiCT=mcevth_2.get_phi_CT(); + if ( (thetaCT == thetashw) && (phiCT == phishw) ) { - } else { - - // the shower comes off-axis - - // obtain with this the final direction of the CT - - thetaCT += thetashw; - phiCT += phishw; - - // calculate vector for telescope - - l2 = sin(thetaCT)*cos(phiCT); - m2 = sin(thetaCT)*sin(phiCT); - n2 = cos(thetaCT); - - num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY); - den = (SQR(l1*m2 - l2*m1) + - SQR(m1*n2 - m2*n1) + - SQR(n1*l2 - n2*l1)); - den = sqrt(den); - - impactD = fabs(num)/den; - + // CT was looking to the source (both lines are parallel) + // therefore, we calculate the impact parameter as the distance + // between the CT axis and the core position + + impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. ); + + } else { + + // the shower comes off-axis + + // calculate vector for telescope + + l2 = sin(thetaCT)*cos(phiCT); + m2 = sin(thetaCT)*sin(phiCT); + n2 = cos(thetaCT); + + num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY); + den = (SQR(l1*m2 - l2*m1) + + SQR(m1*n2 - m2*n1) + + SQR(n1*l2 - n2*l1)); + den = sqrt(den); + + impactD = fabs(num)/den; + + } } @@ -1465,5 +1506,4 @@ rho); - } @@ -1551,11 +1591,20 @@ // Set to zero the flag to know if some conditon has triggered btrigger=0; + flagstoring = 0; // Loop over trigger threshold int iconcount; for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){ - for (i=0;i=qThreshold[i])? (Float_t)(fthrescount):qThreshold[i]; + + // Rise the discrimnator threshold to avoid huge rates + + if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe) + for(int ii=0;iiriseDiskThres) + fpixelthres[ii]=secureDiskThres; + } Trigger.SetThreshold(fpixelthres); @@ -1665,4 +1714,8 @@ // + if(!flagstoring) + nstoredevents++; + flagstoring = 1; + if (Write_McEvt) { Float_t ftime, ltime; @@ -1681,6 +1734,6 @@ mcevth.get_coreY(), impactD, - mcevth_2.get_theta_CT(), - mcevth_2.get_phi_CT(), + phiCT, + thetaCT, ftime, ltime, @@ -1701,5 +1754,5 @@ else{ Float_t Nmax, t0, tmax, a, b, c, chi2; - mcevth.get_times(&ftime, <ime); + mcevth_2.get_times(&ftime, <ime); chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c); McEvt->Fill( mcevth_2.get_evt_number(), @@ -1715,6 +1768,6 @@ mcevth_2.get_coreY(), impactD, + mcevth_2.get_phi_CT(), mcevth_2.get_theta_CT(), - mcevth_2.get_phi_CT(), ftime, ltime, @@ -1761,4 +1814,12 @@ for (int i=0;i0.0 && simulateNSB && nphe2NSB<=inumphe) + for(int ii=0;iiriseDiskThres) + fpixelthres[ii]=secureDiskThres; + } + Trigger.SetThreshold(fpixelthres); @@ -1802,4 +1863,6 @@ // fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii)); + + nstoredevents++; if (Write_McTrig){ @@ -1851,6 +1914,6 @@ mcevth.get_coreY(), impactD, - mcevth_2.get_theta_CT(), - mcevth_2.get_phi_CT(), + phiCT, + thetaCT, ftime, ltime, @@ -1871,5 +1934,5 @@ else{ Float_t Nmax, t0, tmax, a, b, c, chi2; - mcevth.get_times(&ftime, <ime); + mcevth_2.get_times(&ftime, <ime); chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c); McEvt->Fill( mcevth_2.get_evt_number(), @@ -1885,6 +1948,6 @@ mcevth_2.get_coreY(), impactD, + mcevth_2.get_phi_CT(), mcevth_2.get_theta_CT(), - mcevth_2.get_phi_CT(), ftime, ltime, @@ -2022,17 +2085,33 @@ get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS); - get_teles_axis(&telestheta,&telesphi); + if (reflector_file_version<6){ + get_teles_axis(&telestheta,&telesphi); + mcevth.get_theta_range(&shthetamin, &shthetamax); + mcevth.get_phi_range(&shphimin,&shphimax); + mcevth.get_theta_range(&shthetamin, &shthetamax); + mcevth.get_phi_range(&shphimin,&shphimax); + corsika=UInt_t(mcevth.get_VersionPGM()*1000); + for (int i=0; i< 10;i++) + heights[i]=mcevth.get_HeightLev (i); + rnum=mcevth.get_RunNumber(); + } + else{ + telestheta=mcevth_2.get_theta_CT(); + telesphi=mcevth_2.get_phi_CT(); + mcevth_2.get_theta_range(&shthetamin, &shthetamax); + mcevth_2.get_phi_range(&shphimin,&shphimax); + corsika=UInt_t(mcevth_2.get_VersionPGM()*1000); + for (int i=0; i< 10;i++) + heights[i]=mcevth_2.get_HeightLev (i); + rnum=mcevth_2.get_RunNumber(); + } + get_source_off(&sofftheta,&soffphi); - mcevth.get_theta_range(&shthetamin, &shthetamax); - mcevth.get_phi_range(&shphimin,&shphimax); - corsika=UInt_t(mcevth.get_VersionPGM()*1000); if(!Trigger_Loop) icontrigger=0; time (<ime); ftime = ((Float_t)ltime)/1000; - for (int i=0; i< 10;i++) - heights[i]=mcevth.get_HeightLev (i); - rnum=mcevth.get_RunNumber(); - - McRunHeader->Fill(rnum, + + if (reflector_file_version<6) + McRunHeader->Fill(rnum, (UInt_t) 0, mcevth.get_DateRun(), @@ -2047,5 +2126,5 @@ CAMERA_PIXELS, (UInt_t)ntshow, - (UInt_t)ntrigger, + (UInt_t)nstoredevents, 0, sfRaH, @@ -2070,7 +2149,52 @@ heights, corsika, - (UInt_t)(REFL_VERSION_A*100), + (UInt_t)(reflector_file_version*100), (UInt_t)(VERSION*100), 0); + else + McRunHeader->Fill(rnum, + (UInt_t) 0, + mcevth_2.get_DateRun(), + ftime, + icontrigger, + !Write_All_Images, + Write_McEvt, + Write_McTrig, + Write_McFADC, + Write_RawEvt, + addElecNoise, + CAMERA_PIXELS, + (UInt_t)ntshow, + (UInt_t)nstoredevents, + 0, + sfRaH, + sfRaM, + sfRaS, + sfDeD, + sfDeM, + sfDeS, + meanNSB, + telestheta, + telesphi, + sofftheta, + soffphi, + shthetamax, + shthetamin, + shphimax, + shphimin, + mcevth_2.get_CWaveLower(), + mcevth_2.get_CWaveUpper(), + mcevth_2.get_slope(), + 1, + heights, + corsika, + (UInt_t)(reflector_file_version*100), + (UInt_t)(VERSION*100), + 0); + // Fill some missing values for MRawRunHeader + + RunHeader->SetRunNumber(rnum); + RunHeader->SetNumEvents(nstoredevents); + // Store qe for each PMT in output file @@ -4028,4 +4152,10 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.50 2003/01/07 16:33:31 blanch +// Star Field Rotation has been implemented by Raul Orduna. Now there is a +// rotation for each shower. It is done by a non enter pixel rotation assuming +// a circular symetry of the camera. It is not exact but it is accurate enough and +// much faster. +// // Revision 1.49 2002/12/13 10:04:07 blanch // *** empty log message *** @@ -4179,4 +4309,10 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.50 2003/01/07 16:33:31 blanch +// Star Field Rotation has been implemented by Raul Orduna. Now there is a +// rotation for each shower. It is done by a non enter pixel rotation assuming +// a circular symetry of the camera. It is not exact but it is accurate enough and +// much faster. +// // Revision 1.49 2002/12/13 10:04:07 blanch // *** empty log message ***