Index: trunk/MagicSoft/Simulation/Detector/Camera/MDiag.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/MDiag.cxx (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/MDiag.cxx (revision 308) @@ -0,0 +1,62 @@ +//!///////////////////////////////////////////////////////////////////// +// +// camera +// +// @file MDiag.cxx +// @title definition of the diagnostic output root tree +// @author DP +// @date Thu Oct 21 +// +//---------------------------------------------------------------------- +// +// $RCSfile: MDiag.cxx,v $ +// $Revision: 1.1.1.1 $ +// $Author: harald $ +// $Date: 1999-11-05 11:59:32 $ +// +//////////////////////////////////////////////////////////////////////// + +#include "MDiag.h" + +ClassImp(MDiagEventobject) + +MDiagEventobject::MDiagEventobject(){ // the constructor function + // initialize with invalid values + n = -99.; + primary = -99.; + energy = -99.; + cored = -99.; + impact = -99.; + xcore = -99.; + ycore = -99.; + theta = -99.; + phi = -99.; + deviations = -99.; + dtheta = -99.; + dphi = -99.; + trigger = -99.; + ncphs = -99.; + nphes = -99.; + maxpassthr_phe = -99.; + nphes2 = -99.; + length = -99.; + width = -99.; + dist = -99.; + xdist = -99.; + azw = -99.; + miss = -99.; + alpha = -99.; + conc2 = -99.; + conc3 = -99.; + conc4 = -99.; + conc5 = -99.; + conc6 = -99.; + conc7 = -99.; + conc8 = -99.; + conc9 = -99.; + conc10 = -99.; + asymx = -99.; + asymy = -99.; + phiasym = -99.; + +} Index: trunk/MagicSoft/Simulation/Detector/Camera/MDiag.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/MDiag.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/MDiag.h (revision 308) @@ -0,0 +1,74 @@ +//!///////////////////////////////////////////////////////////////////// +// +// camera +// +// @file MDiag.h +// @title definition of the diagnostic output root tree +// @author DP +// @date Thu Oct 21 +// +//---------------------------------------------------------------------- +// +// $RCSfile: MDiag.h,v $ +// $Revision: 1.1.1.1 $ +// $Author: harald $ +// $Date: 1999-11-05 11:59:33 $ +// +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include "TROOT.h" +#include "TFile.h" +#include "TRandom.h" +#include "TTree.h" + +class MDiagEventobject : public TObject{ +public: + + Float_t n; + Float_t primary; + Float_t energy; + Float_t cored; + Float_t impact; + Float_t xcore; + Float_t ycore; + Float_t theta; + Float_t phi; + Float_t deviations; + Float_t dtheta; + Float_t dphi; + Float_t trigger; + Float_t ncphs; + Float_t maxpassthr_phe; + Float_t nphes; + Float_t nphes2; + Float_t length; + Float_t width; + Float_t dist; + Float_t xdist; + Float_t azw; + Float_t miss; + Float_t alpha; + Float_t conc2; + Float_t conc3; + Float_t conc4; + Float_t conc5; + Float_t conc6; + Float_t conc7; + Float_t conc8; + Float_t conc9; + Float_t conc10; + Float_t asymx; + Float_t asymy; + Float_t phiasym; + + MDiagEventobject(); + virtual ~MDiagEventobject() {} + + ClassDef(MDiagEventobject,1) + +}; Index: trunk/MagicSoft/Simulation/Detector/Camera/MDiagLinkDef.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/MDiagLinkDef.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/MDiagLinkDef.h (revision 308) @@ -0,0 +1,9 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MDiagEventobject; + +#endif Index: trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.cxx (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.cxx (revision 308) @@ -0,0 +1,588 @@ +/******************************************************** +* MDiagdict.cxx +********************************************************/ +#include "MDiagdict.h" + +#ifdef G__MEMTEST +#undef malloc +#endif + +extern "C" void G__cpp_reset_tagtableMDiagdict(); + +extern "C" void G__set_cpp_environmentMDiagdict() { + G__add_macro("__MAKECINT__"); + G__add_compiledheader("TROOT.h"); + G__add_compiledheader("TMemberInspector.h"); + G__add_compiledheader("MDiag.h"); + G__cpp_reset_tagtableMDiagdict(); +} +extern "C" int G__cpp_dllrevMDiagdict() { return(51111); } + +/********************************************************* +* Member function Interface Method +*********************************************************/ + +/* MDiagEventobject */ +static int G__MDiagEventobject_MDiagEventobject_0_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + MDiagEventobject *p=NULL; + if(G__getaryconstruct()) p=new MDiagEventobject[G__getaryconstruct()]; + else p=new MDiagEventobject; + result7->obj.i = (long)p; + result7->ref = (long)p; + result7->type = 'u'; + result7->tagnum = G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_DeclFileName_2_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,67,(long)((MDiagEventobject*)(G__getstructoffset()))->DeclFileName()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_DeclFileLine_3_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,105,(long)((MDiagEventobject*)(G__getstructoffset()))->DeclFileLine()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_ImplFileName_4_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,67,(long)((MDiagEventobject*)(G__getstructoffset()))->ImplFileName()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_ImplFileLine_5_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,105,(long)((MDiagEventobject*)(G__getstructoffset()))->ImplFileLine()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_Class_Version_6_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,115,(long)((MDiagEventobject*)(G__getstructoffset()))->Class_Version()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_Class_7_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,85,(long)((MDiagEventobject*)(G__getstructoffset()))->Class()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_Dictionary_8_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__setnull(result7); + ((MDiagEventobject*)(G__getstructoffset()))->Dictionary(); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_IsA_9_0(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__letint(result7,85,(long)((MDiagEventobject*)(G__getstructoffset()))->IsA()); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_ShowMembers_0_1(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__setnull(result7); + ((MDiagEventobject*)(G__getstructoffset()))->ShowMembers(*(TMemberInspector*)libp->para[0].ref,(char*)G__int(libp->para[1])); + return(1 || funcname || hash || result7 || libp) ; +} + +static int G__MDiagEventobject_Streamer_1_1(G__value *result7,char *funcname,struct G__param *libp,int hash) { + G__setnull(result7); + ((MDiagEventobject*)(G__getstructoffset()))->Streamer(*(TBuffer*)libp->para[0].ref); + return(1 || funcname || hash || result7 || libp) ; +} + +// automatic copy constructor +static int G__MDiagEventobject_MDiagEventobject_2_1(G__value *result7,char *funcname,struct G__param *libp,int hash) +{ + MDiagEventobject *p; + if(1!=libp->paran) ; + p=new MDiagEventobject(*(MDiagEventobject*)G__int(libp->para[0])); + result7->obj.i = (long)p; + result7->ref = (long)p; + result7->type = 'u'; + result7->tagnum = G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject); + return(1 || funcname || hash || result7 || libp) ; +} + +// automatic destructor +static int G__MDiagEventobject_wAMDiagEventobject_3_1(G__value *result7,char *funcname,struct G__param *libp,int hash) { + if(G__getaryconstruct()) + if(G__PVOID==G__getgvp()) + delete[] (MDiagEventobject *)(G__getstructoffset()); + else + for(int i=G__getaryconstruct()-1;i>=0;i--) + delete (MDiagEventobject *)((G__getstructoffset())+sizeof(MDiagEventobject)*i); + else delete (MDiagEventobject *)(G__getstructoffset()); + G__setnull(result7); + return(1 || funcname || hash || result7 || libp) ; +} + + +/* Setting up global function */ + +/********************************************************* +* Member function Stub +*********************************************************/ + +/* MDiagEventobject */ + +/********************************************************* +* Global function Stub +*********************************************************/ + +/********************************************************* +* Get size of pointer to member function +*********************************************************/ +class G__Sizep2memfuncMDiagdict { + public: + G__Sizep2memfuncMDiagdict() {p=&G__Sizep2memfuncMDiagdict::sizep2memfunc;} + size_t sizep2memfunc() { return(sizeof(p)); } + private: + size_t (G__Sizep2memfuncMDiagdict::*p)(); +}; + +size_t G__get_sizep2memfuncMDiagdict() +{ + G__Sizep2memfuncMDiagdict a; + G__setsizep2memfunc((int)a.sizep2memfunc()); + return((size_t)a.sizep2memfunc()); +} + + +/********************************************************* +* virtual base class offset calculation interface +*********************************************************/ + + /* Setting up class inheritance */ + +/********************************************************* +* Inheritance information setup/ +*********************************************************/ +extern "C" void G__cpp_setup_inheritanceMDiagdict() { + + /* Setting up class inheritance */ + if(0==G__getnumbaseclass(G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject))) { + MDiagEventobject *G__Lderived; + G__Lderived=(MDiagEventobject*)0x1000; + { + TObject *G__Lpbase=(TObject*)G__Lderived; + G__inheritance_setup(G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject),G__get_linked_tagnum(&G__MDiagdictLN_TObject),(long)G__Lpbase-(long)G__Lderived,1,1); + } + } +} + +/********************************************************* +* typedef information setup/ +*********************************************************/ +extern "C" void G__cpp_setup_typetableMDiagdict() { + + /* Setting up typedef entry */ + G__search_typename2("Char_t",99,-1,0, +-1); + G__setnewtype(-1,"Signed Character 1 byte",0); + G__search_typename2("UChar_t",98,-1,0, +-1); + G__setnewtype(-1,"Unsigned Character 1 byte",0); + G__search_typename2("Short_t",115,-1,0, +-1); + G__setnewtype(-1,"Signed Short integer 2 bytes",0); + G__search_typename2("UShort_t",114,-1,0, +-1); + G__setnewtype(-1,"Unsigned Short integer 2 bytes",0); + G__search_typename2("Int_t",105,-1,0, +-1); + G__setnewtype(-1,"Signed integer 4 bytes",0); + G__search_typename2("UInt_t",104,-1,0, +-1); + G__setnewtype(-1,"Unsigned integer 4 bytes",0); + G__search_typename2("Seek_t",105,-1,0, +-1); + G__setnewtype(-1,"File pointer",0); + G__search_typename2("Long_t",108,-1,0, +-1); + G__setnewtype(-1,"Signed long integer 8 bytes",0); + G__search_typename2("ULong_t",107,-1,0, +-1); + G__setnewtype(-1,"Unsigned long integer 8 bytes",0); + G__search_typename2("Float_t",102,-1,0, +-1); + G__setnewtype(-1,"Float 4 bytes",0); + G__search_typename2("Double_t",100,-1,0, +-1); + G__setnewtype(-1,"Float 8 bytes",0); + G__search_typename2("Text_t",99,-1,0, +-1); + G__setnewtype(-1,"General string",0); + G__search_typename2("Bool_t",98,-1,0, +-1); + G__setnewtype(-1,"Boolean (0=false, 1=true)",0); + G__search_typename2("Byte_t",98,-1,0, +-1); + G__setnewtype(-1,"Byte (8 bits)",0); + G__search_typename2("Version_t",115,-1,0, +-1); + G__setnewtype(-1,"Class version identifier",0); + G__search_typename2("Option_t",99,-1,0, +-1); + G__setnewtype(-1,"Option string",0); + G__search_typename2("Ssiz_t",105,-1,0, +-1); + G__setnewtype(-1,"String size",0); + G__search_typename2("Real_t",102,-1,0, +-1); + G__setnewtype(-1,"TVector and TMatrix element type",0); + G__search_typename2("VoidFuncPtr_t",89,-1,0, +-1); + G__setnewtype(-1,"pointer to void function",0); + G__search_typename2("FreeHookFun_t",89,-1,0, +-1); + G__setnewtype(-1,NULL,0); + G__search_typename2("ReAllocFun_t",81,-1,0, +-1); + G__setnewtype(-1,NULL,0); + G__search_typename2("ReAllocCFun_t",81,-1,0, +-1); + G__setnewtype(-1,NULL,0); + G__search_typename2("Axis_t",102,-1,0, +-1); + G__setnewtype(-1,"Axis values type",0); + G__search_typename2("Stat_t",100,-1,0, +-1); + G__setnewtype(-1,"Statistics type",0); + G__search_typename2("Font_t",115,-1,0, +-1); + G__setnewtype(-1,"Font number",0); + G__search_typename2("Style_t",115,-1,0, +-1); + G__setnewtype(-1,"Style number",0); + G__search_typename2("Marker_t",115,-1,0, +-1); + G__setnewtype(-1,"Marker number",0); + G__search_typename2("Width_t",115,-1,0, +-1); + G__setnewtype(-1,"Line width",0); + G__search_typename2("Color_t",115,-1,0, +-1); + G__setnewtype(-1,"Color number",0); + G__search_typename2("SCoord_t",115,-1,0, +-1); + G__setnewtype(-1,"Screen coordinates",0); + G__search_typename2("Coord_t",102,-1,0, +-1); + G__setnewtype(-1,"Pad world coordinates",0); + G__search_typename2("Angle_t",102,-1,0, +-1); + G__setnewtype(-1,"Graphics angle",0); + G__search_typename2("Size_t",102,-1,0, +-1); + G__setnewtype(-1,"Attribute size",0); +} + +/********************************************************* +* Data Member information setup/ +*********************************************************/ + + /* Setting up class,struct,union tag member variable */ + + /* MDiagEventobject */ +static void G__setup_memvarMDiagEventobject(void) { + G__tag_memvar_setup(G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject)); + { MDiagEventobject *p; p=(MDiagEventobject*)0x1000; if (p) { } + G__memvar_setup((void*)((long)(&p->n)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"n=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->primary)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"primary=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->energy)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"energy=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->cored)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"cored=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->impact)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"impact=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->xcore)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"xcore=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->ycore)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"ycore=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->theta)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"theta=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->phi)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"phi=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->deviations)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"deviations=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->dtheta)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"dtheta=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->dphi)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"dphi=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->trigger)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"trigger=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->ncphs)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"ncphs=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->maxpassthr_phe)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"maxpassthr_phe=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->nphes)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"nphes=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->nphes2)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"nphes2=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->length)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"length=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->width)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"width=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->dist)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"dist=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->xdist)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"xdist=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->azw)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"azw=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->miss)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"miss=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->alpha)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"alpha=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc2)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc2=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc3)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc3=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc4)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc4=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc5)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc5=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc6)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc6=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc7)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc7=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc8)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc8=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc9)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc9=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->conc10)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"conc10=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->asymx)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"asymx=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->asymy)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"asymy=",0,(char*)NULL); + G__memvar_setup((void*)((long)(&p->phiasym)-(long)(p)),102,0,0,-1,G__defined_typename("Float_t"),-1,1,"phiasym=",0,(char*)NULL); + G__memvar_setup((void*)NULL,85,0,0,G__get_linked_tagnum(&G__MDiagdictLN_TClass),-1,-2,4,"fgIsA=",0,(char*)NULL); + } + G__tag_memvar_reset(); +} + +extern "C" void G__cpp_setup_memvarMDiagdict() { +} +/*********************************************************** +************************************************************ +************************************************************ +************************************************************ +************************************************************ +************************************************************ +************************************************************ +***********************************************************/ + +/********************************************************* +* Member function information setup for each class +*********************************************************/ +static void G__setup_memfuncMDiagEventobject(void) { + /* MDiagEventobject */ + G__tag_memfunc_setup(G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject)); + G__memfunc_setup("MDiagEventobject",1595,G__MDiagEventobject_MDiagEventobject_0_0,105,G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject),-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("DeclFileName",1145,G__MDiagEventobject_DeclFileName_2_0,67,-1,-1,0,0,1,1,1,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("DeclFileLine",1152,G__MDiagEventobject_DeclFileLine_3_0,105,-1,-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("ImplFileName",1171,G__MDiagEventobject_ImplFileName_4_0,67,-1,-1,0,0,1,1,1,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("ImplFileLine",1178,G__MDiagEventobject_ImplFileLine_5_0,105,-1,-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("Class_Version",1339,G__MDiagEventobject_Class_Version_6_0,115,-1,G__defined_typename("Version_t"),0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("Class",502,G__MDiagEventobject_Class_7_0,85,G__get_linked_tagnum(&G__MDiagdictLN_TClass),-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("Dictionary",1046,G__MDiagEventobject_Dictionary_8_0,121,-1,-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,0); + G__memfunc_setup("IsA",253,G__MDiagEventobject_IsA_9_0,85,G__get_linked_tagnum(&G__MDiagdictLN_TClass),-1,0,0,1,1,8,"",(char*)NULL,(void*)NULL,1); + G__memfunc_setup("ShowMembers",1132,G__MDiagEventobject_ShowMembers_0_1,121,-1,-1,0,2,1,1,0, +"u 'TMemberInspector' - 1 - insp C - - 0 - parent",(char*)NULL,(void*)NULL,1); + G__memfunc_setup("Streamer",835,G__MDiagEventobject_Streamer_1_1,121,-1,-1,0,1,1,1,0,"u 'TBuffer' - 1 - b",(char*)NULL,(void*)NULL,1); + // automatic copy constructor + G__memfunc_setup("MDiagEventobject",1595,G__MDiagEventobject_MDiagEventobject_2_1,(int)('i'),G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject),-1,0,1,1,1,0,"u 'MDiagEventobject' - 1 - -",(char*)NULL,(void*)NULL,0); + // automatic destructor + G__memfunc_setup("~MDiagEventobject",1721,G__MDiagEventobject_wAMDiagEventobject_3_1,(int)('y'),-1,-1,0,0,1,1,0,"",(char*)NULL,(void*)NULL,1); + G__tag_memfunc_reset(); +} + + +/********************************************************* +* Member function information setup +*********************************************************/ +extern "C" void G__cpp_setup_memfuncMDiagdict() { +} + +/********************************************************* +* Global variable information setup for each class +*********************************************************/ +extern "C" void G__cpp_setup_globalMDiagdict() { + + /* Setting up global variables */ + G__resetplocal(); + + + G__resetglobalenv(); +} + +/********************************************************* +* Global function information setup for each class +*********************************************************/ +extern "C" void G__cpp_setup_funcMDiagdict() { + G__lastifuncposition(); + + + G__resetifuncposition(); +} + +/********************************************************* +* Class,struct,union,enum tag information setup +*********************************************************/ +/* Setup class/struct taginfo */ +G__linked_taginfo G__MDiagdictLN_TClass = { "TClass" , 99 , -1 }; +G__linked_taginfo G__MDiagdictLN_TObject = { "TObject" , 99 , -1 }; +G__linked_taginfo G__MDiagdictLN_MDiagEventobject = { "MDiagEventobject" , 99 , -1 }; + +/* Reset class/struct taginfo */ +extern "C" void G__cpp_reset_tagtableMDiagdict() { + G__MDiagdictLN_TClass.tagnum = -1 ; + G__MDiagdictLN_TObject.tagnum = -1 ; + G__MDiagdictLN_MDiagEventobject.tagnum = -1 ; +} + +extern "C" void G__cpp_setup_tagtableMDiagdict() { + + /* Setting up class,struct,union tag entry */ + G__tagtable_setup(G__get_linked_tagnum(&G__MDiagdictLN_MDiagEventobject),sizeof(MDiagEventobject),-1,0,(char*)NULL,G__setup_memvarMDiagEventobject,G__setup_memfuncMDiagEventobject); +} +extern "C" void G__cpp_setupMDiagdict() { + G__check_setup_version(51111,"G__cpp_setupMDiagdict()"); + G__set_cpp_environmentMDiagdict(); + G__cpp_setup_tagtableMDiagdict(); + + G__cpp_setup_inheritanceMDiagdict(); + + G__cpp_setup_typetableMDiagdict(); + + G__cpp_setup_memvarMDiagdict(); + + G__cpp_setup_memfuncMDiagdict(); + G__cpp_setup_globalMDiagdict(); + G__cpp_setup_funcMDiagdict(); + + if(0==G__getsizep2memfunc()) G__get_sizep2memfuncMDiagdict(); + return; +} +class G__cpp_setup_initMDiagdict { + public: + G__cpp_setup_initMDiagdict() { G__add_setup_func("MDiagdict",&G__cpp_setupMDiagdict); } + ~G__cpp_setup_initMDiagdict() { G__remove_setup_func("MDiagdict"); } +}; +G__cpp_setup_initMDiagdict G__cpp_setup_initializerMDiagdict; + +// +// File generated by /cern/root/bin/rootcint at Fri Oct 29 11:52:06 1999. +// Do NOT change. Changes will be lost next time file is generated +// + +#include "TBuffer.h" +#include "TMemberInspector.h" +#include "TError.h" + +//______________________________________________________________________________ +TBuffer &operator>>(TBuffer &buf, MDiagEventobject *&obj) +{ + // Read a pointer to an object of class MDiagEventobject. + + obj = (MDiagEventobject *) buf.ReadObject(MDiagEventobject::Class()); + return buf; +} + +//______________________________________________________________________________ +void MDiagEventobject::Streamer(TBuffer &R__b) +{ + // Stream an object of class MDiagEventobject. + + if (R__b.IsReading()) { + Version_t R__v = R__b.ReadVersion(); if (R__v) { } + TObject::Streamer(R__b); + R__b >> n; + R__b >> primary; + R__b >> energy; + R__b >> cored; + R__b >> impact; + R__b >> xcore; + R__b >> ycore; + R__b >> theta; + R__b >> phi; + R__b >> deviations; + R__b >> dtheta; + R__b >> dphi; + R__b >> trigger; + R__b >> ncphs; + R__b >> maxpassthr_phe; + R__b >> nphes; + R__b >> nphes2; + R__b >> length; + R__b >> width; + R__b >> dist; + R__b >> xdist; + R__b >> azw; + R__b >> miss; + R__b >> alpha; + R__b >> conc2; + R__b >> conc3; + R__b >> conc4; + R__b >> conc5; + R__b >> conc6; + R__b >> conc7; + R__b >> conc8; + R__b >> conc9; + R__b >> conc10; + R__b >> asymx; + R__b >> asymy; + R__b >> phiasym; + } else { + R__b.WriteVersion(MDiagEventobject::IsA()); + TObject::Streamer(R__b); + R__b << n; + R__b << primary; + R__b << energy; + R__b << cored; + R__b << impact; + R__b << xcore; + R__b << ycore; + R__b << theta; + R__b << phi; + R__b << deviations; + R__b << dtheta; + R__b << dphi; + R__b << trigger; + R__b << ncphs; + R__b << maxpassthr_phe; + R__b << nphes; + R__b << nphes2; + R__b << length; + R__b << width; + R__b << dist; + R__b << xdist; + R__b << azw; + R__b << miss; + R__b << alpha; + R__b << conc2; + R__b << conc3; + R__b << conc4; + R__b << conc5; + R__b << conc6; + R__b << conc7; + R__b << conc8; + R__b << conc9; + R__b << conc10; + R__b << asymx; + R__b << asymy; + R__b << phiasym; + } +} + +//______________________________________________________________________________ +void MDiagEventobject::ShowMembers(TMemberInspector &R__insp, char *R__parent) +{ + // Inspect the data members of an object of class MDiagEventobject. + + TClass *R__cl = MDiagEventobject::IsA(); + Int_t R__ncp = strlen(R__parent); + if (R__ncp || R__cl || R__insp.IsA()) { } + R__insp.Inspect(R__cl, R__parent, "n", &n); + R__insp.Inspect(R__cl, R__parent, "primary", &primary); + R__insp.Inspect(R__cl, R__parent, "energy", &energy); + R__insp.Inspect(R__cl, R__parent, "cored", &cored); + R__insp.Inspect(R__cl, R__parent, "impact", &impact); + R__insp.Inspect(R__cl, R__parent, "xcore", &xcore); + R__insp.Inspect(R__cl, R__parent, "ycore", &ycore); + R__insp.Inspect(R__cl, R__parent, "theta", &theta); + R__insp.Inspect(R__cl, R__parent, "phi", &phi); + R__insp.Inspect(R__cl, R__parent, "deviations", &deviations); + R__insp.Inspect(R__cl, R__parent, "dtheta", &dtheta); + R__insp.Inspect(R__cl, R__parent, "dphi", &dphi); + R__insp.Inspect(R__cl, R__parent, "trigger", &trigger); + R__insp.Inspect(R__cl, R__parent, "ncphs", &ncphs); + R__insp.Inspect(R__cl, R__parent, "maxpassthr_phe", &maxpassthr_phe); + R__insp.Inspect(R__cl, R__parent, "nphes", &nphes); + R__insp.Inspect(R__cl, R__parent, "nphes2", &nphes2); + R__insp.Inspect(R__cl, R__parent, "length", &length); + R__insp.Inspect(R__cl, R__parent, "width", &width); + R__insp.Inspect(R__cl, R__parent, "dist", &dist); + R__insp.Inspect(R__cl, R__parent, "xdist", &xdist); + R__insp.Inspect(R__cl, R__parent, "azw", &azw); + R__insp.Inspect(R__cl, R__parent, "miss", &miss); + R__insp.Inspect(R__cl, R__parent, "alpha", &alpha); + R__insp.Inspect(R__cl, R__parent, "conc2", &conc2); + R__insp.Inspect(R__cl, R__parent, "conc3", &conc3); + R__insp.Inspect(R__cl, R__parent, "conc4", &conc4); + R__insp.Inspect(R__cl, R__parent, "conc5", &conc5); + R__insp.Inspect(R__cl, R__parent, "conc6", &conc6); + R__insp.Inspect(R__cl, R__parent, "conc7", &conc7); + R__insp.Inspect(R__cl, R__parent, "conc8", &conc8); + R__insp.Inspect(R__cl, R__parent, "conc9", &conc9); + R__insp.Inspect(R__cl, R__parent, "conc10", &conc10); + R__insp.Inspect(R__cl, R__parent, "asymx", &asymx); + R__insp.Inspect(R__cl, R__parent, "asymy", &asymy); + R__insp.Inspect(R__cl, R__parent, "phiasym", &phiasym); + TObject::ShowMembers(R__insp, R__parent); +} + Index: trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/MDiagdict.h (revision 308) @@ -0,0 +1,35 @@ +/******************************************************************** +* MDiagdict.h +********************************************************************/ +#ifdef __CINT__ +#error MDiagdict.h/C is only for compilation. Abort cint. +#endif +#include +#include +#include +#include +#include +extern "C" { +#define G__ANSIHEADER +#include "G__ci.h" +extern void G__cpp_setup_tagtableMDiagdict(); +extern void G__cpp_setup_inheritanceMDiagdict(); +extern void G__cpp_setup_typetableMDiagdict(); +extern void G__cpp_setup_memvarMDiagdict(); +extern void G__cpp_setup_globalMDiagdict(); +extern void G__cpp_setup_memfuncMDiagdict(); +extern void G__cpp_setup_funcMDiagdict(); +extern void G__set_cpp_environmentMDiagdict(); +} + + +#include "TROOT.h" +#include "TMemberInspector.h" +#include "MDiag.h" + +#ifndef G__MEMFUNCBODY +#endif + +extern G__linked_taginfo G__MDiagdictLN_TClass; +extern G__linked_taginfo G__MDiagdictLN_TObject; +extern G__linked_taginfo G__MDiagdictLN_MDiagEventobject; Index: trunk/MagicSoft/Simulation/Detector/Camera/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/Makefile (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/Makefile (revision 308) @@ -0,0 +1,427 @@ +################################################################## +# +# makefile +# +# @file makefile +# @title Simulation of the camera and trigger logic +# @author J C Gonz\'alez +# @email gonzalez@mppmu.mpg.de +# @date Fri Mar 12 11:51:11 MET 1999 +# +#_______________________________________________________________ +# +# Created: Fri Mar 12 11:51:11 MET 1999 +# Author: Jose Carlos Gonzalez +# Purpose: Makefile for the compilation of the camera program +# Notes: +# +#--------------------------------------------------------------- +# +# $RCSfile: Makefile,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 1999-11-05 11:59:32 $ +# +################################################################## +# @maintitle + +# @code + +INCLUDEMK = config.mk.${OSTYPE} +include ${INCLUDEMK} + +# @endcode + +# @code + +# common flags +INCLUDES = -I${INCLUDE} \ + -I${INCLUDE_COR} \ + -I${INCLUDE_MC} \ + -I${INCLUDE_EVITA} \ + -I${INCLUDE_TRIGGER} \ + -I${INCLUDE_REFL} \ + -I${INCLUDE_ROOT} \ + -I/usr/include/g++ + +RANLIB = -L${RANLIBDIR} -lranlib + +# what is needed for ROOT + +ROOTLIBS = -L$(ROOTSYS)/lib -lNew -lBase -lCint -lClib \ + -lCont -lFunc -lGraf -lGraf3d -lHist -lHtml \ + -lMatrix -lMeta -lMinuit -lNet -lPostscript \ + -lProof -lTree -lUnix -lZip -lRint +#ROOTLIBS = + +ROOTGLIBS = -lGpad -lGui -lGX11 -lX3d -lX11 +#ROOTGLIBS = + +GLIBS = $(ROOTLIBS) $(ROOTGLIBS) -L/usr/X11R6/lib \ + -lXpm -lX11 -lm -ldl -rdynamic +# special flags + +osf_FORLIBS = -lUfor -lfor -lutil -lots -lm +#linux_FORLIBS = -lf2c -lm /usr/lib/libc.a +linux_FORLIBS = -lm -ldl +#linux_FORLIBS = -lm -ldl -rdynamic +generic_FORLIBS = -lm + +FORLIBS = ${${SYSTEM}_FORLIBS} + +# compilation and linking flags + +CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} + CFLAGS = ${CXXFLAGS} +FFLAGS = ${CXXFLAGS} +LIBS = ${RANLIB} ${ROOTLIBS} ${ROOTGLIBS} ${GLIBS} + +#------------------------------------------------------------------------------ + +#.SILENT: + +.SUFFIXES: .c .cxx .C .c++ .h .hxx .H .h++ .o .so .f + +SRCS = \ + ${INCLUDE_MC}/MCEventHeader.cxx \ + ${INCLUDE_MC}/MCCphoton.cxx \ + ${INCLUDE_TRIGGER}/MTrigger.cxx \ + ${INCLUDE_EVITA}/MRawPixel.cxx \ + ${INCLUDE_EVITA}/MRawEvt.cxx \ + ${INCLUDE_EVITA}/MMcEvt.cxx \ + MCamCint.cxx \ + MDiag.cxx \ + moments.cxx \ + creadparam.cxx \ + camera.cxx + +HEADERS = \ + MCEventHeader.hxx \ + MCCphoton.hxx \ + MTRigger.hxx \ + MRawPixel.h \ + MRawEvt.h \ + MMcEvt.h \ + lagrange.h \ + atm.h \ + MDiag.h \ + moments.h \ + creadparam.h \ + camera.h + +OBJS = \ + ${INCLUDE_MC}/MCEventHeader.o \ + ${INCLUDE_MC}/MCCphoton.o \ + ${INCLUDE_TRIGGER}/MTrigger.o \ + ${INCLUDE_EVITA}/MRawPixel.o \ + ${INCLUDE_EVITA}/MRawEvt.o \ + ${INCLUDE_EVITA}/MMcEvt.o \ + MCamCint.o \ + MDiag.o \ + MDiagdict.o \ + moments.o \ + creadparam.o \ + camera.o + +############################################################ + +all: ${PROGRAM} + +depend: + @makedepend $(SRCS) $(INCLUDES) -fMakefile 2> kk.kk ; cat kk.kk + +doc: camera-doc + +camera-doc: + @echo "Generating documentation for camera . . . " + $(DOCUM) -latex -o camera.tex \ + camera.cxx camera.h \ + creadparam.cxx creadparam.h \ + moments.cxx moments.h + latex "\nonstopmode\input{camera.tex}" && \ + makeindex camera && \ + latex "\nonstopmode\input{camera.tex}" && \ + latex "\nonstopmode\input{camera.tex}" + @echo "Files camera.tex and camera.dvi generated." + +${PROGRAM}: $(OBJS) MDiag.so + @echo "Linking..." + echo `ls -m $(OBJS)|sed 's/,/ +/g' `" + libraries => " $@ + $(CXX) $(CXXFLAGS) $(OBJS) $(LIBS) -o $@ + @echo "done." + +MDiagdict.o: MDiagdict.cxx MDiag.h + $(CXX) $(CXXFLAGS) -I${INCLUDE_ROOT} -c MDiagdict.cxx + +MDiagdict.cxx: MDiag.h MDiagLinkDef.h + ${ROOTSYS}/bin/rootcint -f MDiagdict.cxx -c MDiag.h MDiagLinkDef.h + +MDiag.o: MDiag.cxx MDiag.h + $(CXX) $(CXXFLAGS) -I${INCLUDE_ROOT} -c MDiag.cxx -o MDiag.o + +# the following shared object library is for being loaded into ROOT using +#".L MDiag.so" if the diagnostic output is to be read e.g. with a TTree viewer +MDiag.so: MDiag.o MDiagdict.o + $(CXX) -shared -g -I${INCLUDE_ROOT} MDiag.o MDiagdict.o -o MDiag.so + + +MCamCint.cxx: ${INCLUDE_EVITA}/MRawPixel.h \ + ${INCLUDE_EVITA}/MRawEvt.h \ + ${INCLUDE_EVITA}/MMcEvt.h \ + ${INCLUDE_EVITA}/Mdefine.h + + @echo + @echo "Generating dictionary ..." + @echo + + @$(ROOTSYS)/bin/rootcint -f \ + MCamCint.cxx -c \ + ${INCLUDE_EVITA}/MRawPixel.h \ + ${INCLUDE_EVITA}/MRawEvt.h \ + ${INCLUDE_EVITA}/MMcEvt.h \ +# ${INCLUDE_EVITA}/MCameraDisplay.h \ + ${INCLUDE_EVITA}/Mdefine.h \ + ${INCLUDE_EVITA}/LinkDef.h + + @echo + @echo "Dictionary done" + @echo + +.cxx.o: + @echo "Compiling " $< + $(CXX) $(CXXFLAGS) -c $< -o $@ + +.c.o: + @echo "Compiling " $< + $(CC) $(CFLAGS) -c $< -o $@ + +lclean: + @echo "Cleanning..." + @rm -f *.o core + +clean: + @echo "Cleanning..." + @rm -f $(OBJS) core + @rm -f MCamCint.cxx MCamCint.h + +mrproper: clean + @echo "Mr.Proper in action . . ." + @rm -f $(PROGRAM) + +ctags: + @echo "Creating CTAGS file . . ." + @ctags -txw $(SRCS) $(HEADERS) > CTAGS + +etags: + @echo "Creating TAGS file . . ." + @etags -C $(SRCS) $(HEADERS) + +listsrc: + @ls -m $(SRCS) $(HEADERS) | sed 's/,//g' + +redo: clean all + +cflags: + @echo $(INCLUDES) $(CXXFLAGS) + +# @endcode + +# DO NOT DELETE THIS LINE -- make depend depends on it. + +../include-MC/MCEventHeader.o: ../include-MC/MCEventHeader.hxx +../include-MC/MCEventHeader.o: ../include-GENERAL/Rtypes.h +../include-MC/MCEventHeader.o: /usr/include/g++/iostream.h +../include-MC/MCEventHeader.o: /usr/include/g++/streambuf.h +../include-MC/MCEventHeader.o: /usr/include/libio.h /usr/include/features.h +../include-MC/MCEventHeader.o: /usr/include/sys/cdefs.h +../include-MC/MCEventHeader.o: /usr/include/gnu/stubs.h +../include-MC/MCEventHeader.o: /usr/include/_G_config.h +../include-MC/MCEventHeader.o: /usr/include/gnu/types.h +../include-MC/MCEventHeader.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../include-MC/MCEventHeader.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../include-MC/MCEventHeader.o: /usr/include/g++/iomanip.h +../include-MC/MCEventHeader.o: /usr/include/g++/fstream.h +../include-MC/MCEventHeader.o: /usr/include/stdlib.h /usr/include/sys/types.h +../include-MC/MCEventHeader.o: /usr/include/time.h /usr/include/endian.h +../include-MC/MCEventHeader.o: /usr/include/bytesex.h +../include-MC/MCEventHeader.o: /usr/include/sys/select.h +../include-MC/MCEventHeader.o: /usr/include/selectbits.h +../include-MC/MCEventHeader.o: /usr/include/alloca.h /usr/include/math.h +../include-MC/MCEventHeader.o: /usr/include/huge_val.h +../include-MC/MCEventHeader.o: /usr/include/mathcalls.h +../include-MC/MCEventHeader.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +../include-MC/MCEventHeader.o: ../include-CORSIKA/COREventHeader.hxx +../include-MC/MCCphoton.o: ../include-MC/MCCphoton.hxx +../include-MC/MCCphoton.o: ../include-GENERAL/Rtypes.h +../include-MC/MCCphoton.o: /usr/include/g++/iostream.h +../include-MC/MCCphoton.o: /usr/include/g++/streambuf.h /usr/include/libio.h +../include-MC/MCCphoton.o: /usr/include/features.h /usr/include/sys/cdefs.h +../include-MC/MCCphoton.o: /usr/include/gnu/stubs.h /usr/include/_G_config.h +../include-MC/MCCphoton.o: /usr/include/gnu/types.h +../include-MC/MCCphoton.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../include-MC/MCCphoton.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../include-MC/MCCphoton.o: /usr/include/g++/iomanip.h +../include-MC/MCCphoton.o: /usr/include/g++/fstream.h /usr/include/stdlib.h +../include-MC/MCCphoton.o: /usr/include/sys/types.h /usr/include/time.h +../include-MC/MCCphoton.o: /usr/include/endian.h /usr/include/bytesex.h +../include-MC/MCCphoton.o: /usr/include/sys/select.h +../include-MC/MCCphoton.o: /usr/include/selectbits.h /usr/include/alloca.h +../include-MC/MCCphoton.o: /usr/include/string.h /usr/include/math.h +../include-MC/MCCphoton.o: /usr/include/huge_val.h /usr/include/mathcalls.h +../include-MC/MCCphoton.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +../include-MC/MCCphoton.o: ../include-GENERAL/jcmacros.h +../include-MTrigger/MTrigger.o: ../include-MTrigger/MTrigger.hxx +../include-MTrigger/MTrigger.o: /usr/include/g++/iostream.h +../include-MTrigger/MTrigger.o: /usr/include/g++/streambuf.h +../include-MTrigger/MTrigger.o: /usr/include/libio.h /usr/include/features.h +../include-MTrigger/MTrigger.o: /usr/include/sys/cdefs.h +../include-MTrigger/MTrigger.o: /usr/include/gnu/stubs.h +../include-MTrigger/MTrigger.o: /usr/include/_G_config.h +../include-MTrigger/MTrigger.o: /usr/include/gnu/types.h +../include-MTrigger/MTrigger.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../include-MTrigger/MTrigger.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../include-MTrigger/MTrigger.o: /usr/include/math.h /usr/include/huge_val.h +../include-MTrigger/MTrigger.o: /usr/include/mathcalls.h +../include-MTrigger/MTrigger.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +../include-MTrigger/MTrigger.o: /cern/root/include/TObject.h +../include-MTrigger/MTrigger.o: ../../../include-Classes/Mdefine.h +../../../include-Classes/MRawPixel.o: /usr/include/g++/iostream.h +../../../include-Classes/MRawPixel.o: /usr/include/g++/streambuf.h +../../../include-Classes/MRawPixel.o: /usr/include/libio.h +../../../include-Classes/MRawPixel.o: /usr/include/features.h +../../../include-Classes/MRawPixel.o: /usr/include/sys/cdefs.h +../../../include-Classes/MRawPixel.o: /usr/include/gnu/stubs.h +../../../include-Classes/MRawPixel.o: /usr/include/_G_config.h +../../../include-Classes/MRawPixel.o: /usr/include/gnu/types.h +../../../include-Classes/MRawPixel.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../../../include-Classes/MRawPixel.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../../../include-Classes/MRawPixel.o: /cern/root/include/TClonesArray.h +../../../include-Classes/MRawPixel.o: /cern/root/include/TString.h +../../../include-Classes/MRawPixel.o: /usr/include/string.h +../../../include-Classes/MRawPixel.o: /cern/root/include/TRandom.h +../../../include-Classes/MRawPixel.o: ../../../include-Classes/MRawPixel.h +../../../include-Classes/MRawPixel.o: /cern/root/include/TObject.h +../../../include-Classes/MRawPixel.o: ../../../include-Classes/Mdefine.h +../../../include-Classes/MRawEvt.o: /usr/include/g++/iostream.h +../../../include-Classes/MRawEvt.o: /usr/include/g++/streambuf.h +../../../include-Classes/MRawEvt.o: /usr/include/libio.h +../../../include-Classes/MRawEvt.o: /usr/include/features.h +../../../include-Classes/MRawEvt.o: /usr/include/sys/cdefs.h +../../../include-Classes/MRawEvt.o: /usr/include/gnu/stubs.h +../../../include-Classes/MRawEvt.o: /usr/include/_G_config.h +../../../include-Classes/MRawEvt.o: /usr/include/gnu/types.h +../../../include-Classes/MRawEvt.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../../../include-Classes/MRawEvt.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../../../include-Classes/MRawEvt.o: /cern/root/include/TClonesArray.h +../../../include-Classes/MRawEvt.o: /cern/root/include/TString.h +../../../include-Classes/MRawEvt.o: /usr/include/string.h +../../../include-Classes/MRawEvt.o: /cern/root/include/TRandom.h +../../../include-Classes/MRawEvt.o: ../../../include-Classes/MRawEvt.h +../../../include-Classes/MRawEvt.o: /cern/root/include/TObject.h +../../../include-Classes/MRawEvt.o: ../../../include-Classes/Mdefine.h +../../../include-Classes/MRawEvt.o: ../../../include-Classes/MRawPixel.h +../../../include-Classes/MMcEvt.o: /usr/include/g++/iostream.h +../../../include-Classes/MMcEvt.o: /usr/include/g++/streambuf.h +../../../include-Classes/MMcEvt.o: /usr/include/libio.h +../../../include-Classes/MMcEvt.o: /usr/include/features.h +../../../include-Classes/MMcEvt.o: /usr/include/sys/cdefs.h +../../../include-Classes/MMcEvt.o: /usr/include/gnu/stubs.h +../../../include-Classes/MMcEvt.o: /usr/include/_G_config.h +../../../include-Classes/MMcEvt.o: /usr/include/gnu/types.h +../../../include-Classes/MMcEvt.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +../../../include-Classes/MMcEvt.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +../../../include-Classes/MMcEvt.o: ../../../include-Classes/MMcEvt.h +../../../include-Classes/MMcEvt.o: /usr/include/stdlib.h +../../../include-Classes/MMcEvt.o: /usr/include/sys/types.h +../../../include-Classes/MMcEvt.o: /usr/include/time.h /usr/include/endian.h +../../../include-Classes/MMcEvt.o: /usr/include/bytesex.h +../../../include-Classes/MMcEvt.o: /usr/include/sys/select.h +../../../include-Classes/MMcEvt.o: /usr/include/selectbits.h +../../../include-Classes/MMcEvt.o: /usr/include/alloca.h /usr/include/stdio.h +../../../include-Classes/MMcEvt.o: /usr/include/stdio_lim.h +../../../include-Classes/MMcEvt.o: /usr/include/string.h +../../../include-Classes/MMcEvt.o: /usr/include/unistd.h +../../../include-Classes/MMcEvt.o: /usr/include/posix_opt.h +../../../include-Classes/MMcEvt.o: /usr/include/confname.h +../../../include-Classes/MMcEvt.o: /usr/include/fcntl.h +../../../include-Classes/MMcEvt.o: /usr/include/fcntlbits.h +../../../include-Classes/MMcEvt.o: /cern/root/include/TObject.h +MDiag.o: MDiag.h /usr/include/stdlib.h /usr/include/features.h +MDiag.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +MDiag.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +MDiag.o: /usr/include/sys/types.h /usr/include/gnu/types.h +MDiag.o: /usr/include/time.h /usr/include/endian.h /usr/include/bytesex.h +MDiag.o: /usr/include/sys/select.h /usr/include/selectbits.h +MDiag.o: /usr/include/alloca.h /usr/include/g++/iostream.h +MDiag.o: /usr/include/g++/streambuf.h /usr/include/libio.h +MDiag.o: /usr/include/_G_config.h +MDiag.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +MDiag.o: /usr/include/g++/fstream.h /usr/include/string.h +MDiag.o: /cern/root/include/TROOT.h /cern/root/include/TFile.h +MDiag.o: /cern/root/include/TRandom.h /cern/root/include/TTree.h +moments.o: moments.h /usr/include/g++/iostream.h /usr/include/g++/streambuf.h +moments.o: /usr/include/libio.h /usr/include/features.h +moments.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +moments.o: /usr/include/_G_config.h /usr/include/gnu/types.h +moments.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +moments.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +moments.o: /usr/include/g++/fstream.h /usr/include/stdlib.h +moments.o: /usr/include/sys/types.h /usr/include/time.h /usr/include/endian.h +moments.o: /usr/include/bytesex.h /usr/include/sys/select.h +moments.o: /usr/include/selectbits.h /usr/include/alloca.h +moments.o: /usr/include/stdio.h /usr/include/stdio_lim.h +moments.o: /usr/include/string.h /usr/include/math.h /usr/include/huge_val.h +moments.o: /usr/include/mathcalls.h +moments.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +moments.o: ../include-GENERAL/jcmacros.h ../include-GENERAL/jcdebug.h +moments.o: camera-v.h +creadparam.o: creadparam.h /usr/include/g++/iostream.h +creadparam.o: /usr/include/g++/streambuf.h /usr/include/libio.h +creadparam.o: /usr/include/features.h /usr/include/sys/cdefs.h +creadparam.o: /usr/include/gnu/stubs.h /usr/include/_G_config.h +creadparam.o: /usr/include/gnu/types.h +creadparam.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +creadparam.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +creadparam.o: /usr/include/g++/fstream.h /usr/include/stdlib.h +creadparam.o: /usr/include/sys/types.h /usr/include/time.h +creadparam.o: /usr/include/endian.h /usr/include/bytesex.h +creadparam.o: /usr/include/sys/select.h /usr/include/selectbits.h +creadparam.o: /usr/include/alloca.h /usr/include/stdio.h +creadparam.o: /usr/include/stdio_lim.h /usr/include/string.h +creadparam.o: /usr/include/math.h /usr/include/huge_val.h +creadparam.o: /usr/include/mathcalls.h +creadparam.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +creadparam.o: ../include-GENERAL/jcmacros.h ../include-GENERAL/jcdebug.h +creadparam.o: camera-v.h +camera.o: /cern/root/include/TROOT.h /cern/root/include/TFile.h +camera.o: /cern/root/include/TTree.h /cern/root/include/TBranch.h MDiag.h +camera.o: /usr/include/stdlib.h /usr/include/features.h +camera.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +camera.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stddef.h +camera.o: /usr/include/sys/types.h /usr/include/gnu/types.h +camera.o: /usr/include/time.h /usr/include/endian.h /usr/include/bytesex.h +camera.o: /usr/include/sys/select.h /usr/include/selectbits.h +camera.o: /usr/include/alloca.h /usr/include/g++/iostream.h +camera.o: /usr/include/g++/streambuf.h /usr/include/libio.h +camera.o: /usr/include/_G_config.h +camera.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/stdarg.h +camera.o: /usr/include/g++/fstream.h /usr/include/string.h +camera.o: /cern/root/include/TRandom.h ../include-MTrigger/MTrigger.hxx +camera.o: /usr/include/math.h /usr/include/huge_val.h +camera.o: /usr/include/mathcalls.h +camera.o: /usr/lib/gcc-lib/i486-linux/2.7.2.3/include/float.h +camera.o: /cern/root/include/TObject.h ../../../include-Classes/Mdefine.h +camera.o: ../../../include-Classes/MRawEvt.h +camera.o: /cern/root/include/TClonesArray.h ../../../include-Classes/MMcEvt.h +camera.o: /usr/include/stdio.h /usr/include/stdio_lim.h /usr/include/unistd.h +camera.o: /usr/include/posix_opt.h /usr/include/confname.h +camera.o: /usr/include/fcntl.h /usr/include/fcntlbits.h camera.h +camera.o: /usr/include/dirent.h /usr/include/direntry.h +camera.o: /usr/include/posix1_lim.h /usr/include/local_lim.h +camera.o: /usr/include/linux/limits.h /usr/include/libgen.h camera-v.h +camera.o: ../include-GENERAL/jcmacros.h ../include-GENERAL/jcdebug.h +camera.o: creadparam.h ../Reflector/atm.h ../Reflector/reflector-v.h +camera.o: moments.h ../include-GENERAL/lagrange.h +camera.o: ../include-MC/MCEventHeader.hxx ../include-GENERAL/Rtypes.h +camera.o: /usr/include/g++/iomanip.h ../include-CORSIKA/COREventHeader.hxx +camera.o: ../include-MC/MCCphoton.hxx ../include-GENERAL/ranlib.h Index: trunk/MagicSoft/Simulation/Detector/Camera/camera-v.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera-v.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera-v.h (revision 308) @@ -0,0 +1,12 @@ +#ifndef VERSION + +#define PROGRAM camera +#define VERSION 0.2 + +#define GLUE_prep(x,y) #x" "#y +#define GLUE_postp(x,y) GLUE_prep(x,y) + +const char SIGNATURE[] = GLUE_postp( PROGRAM, VERSION ); + +#endif // ! VERSION + Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 308) @@ -0,0 +1,3131 @@ +//!///////////////////////////////////////////////////////////////////// +// +// camera +// +// @file camera.cxx +// @title Camera simulation +// @subtitle Code for the simulation of the camera phase +// @desc Code for the simulation of the camera of CT1 and MAGIC +// @author J C Gonzalez +// @email gonzalez@mppmu.mpg.de +// @date Thu May 7 16:24:22 1998 +// +//---------------------------------------------------------------------- +// +// Created: Thu May 7 16:24:22 1998 +// Author: Jose Carlos Gonzalez +// Purpose: Program for reflector simulation +// Notes: See files README for details +// +//---------------------------------------------------------------------- +// +// $RCSfile: camera.cxx,v $ +// $Revision: 1.1.1.1 $ +// $Author: harald $ +// $Date: 1999-11-05 11:59:31 $ +// +//////////////////////////////////////////////////////////////////////// +// @tableofcontents @coverpage + +//=----------------------------------------------------------- +//!@section Source code of |camera.cxx|. + +/*!@" + + In this section we show the (commented) code of the program for the + read-out of the output files generated by the simulator of the + reflector, |reflector 0.3|. + + @"*/ + +//=----------------------------------------------------------- +//!@subsection Includes and Global variables definition. + +//!@{ + +// includes for ROOT +// BEWARE: the order matters! + +#include "TROOT.h" + +#include "TFile.h" +#include "TTree.h" +#include "TBranch.h" + +#include "MDiag.h" + +#include "MTrigger.hxx" + +#include "MRawEvt.h" +#include "MMcEvt.h" + +/*!@" + + All the defines are located in the file |camera.h|. + + @"*/ + +#include "camera.h" +//!@} + +/*!@" + + The following set of flags are used in time of compilation. They do + not affect directly the behaviour of the program at run-time + (though, of course, if you disconnected the option for + implementation of the Trigger logic, you will not be able to use any + trigger at all. The 'default' values mean default in the sense of + what you got from the server when you obtained this program. + + @"*/ + +//!@{ + +// flag for debugging (default: OFF ) +#define __DEBUG__ +#undef __DEBUG__ + +// flag for NNT in CT1 camera (default: ON ) +#undef __CT1_NO_NEIGHBOURS__ +#define __CT1_NO_NEIGHBOURS__ + +// flag for calculation of NSB (default: ON ) +#undef __NSB__ +#define __NSB__ + +// flag for calculation of QE for pixels (default: ON ) +#undef __QE__ +#define __QE__ + + +// flag for implementation of DETAIL_TRIGGER (default: ON ) +// +// This is the new implementation of Trigger studies +// It relies on a better simulation of the time stucture +// of the PhotoMultiplier. For more details see the +// documentation of the --> class MTrigger <-- +#define __DETAIL_TRIGGER__ +#undef __DETAIL_TRIGGER__ + +// flag for implementation of TRIGGER (default: ON ) +#undef __TRIGGER__ +#define __TRIGGER__ + +// flag for implementation of Tail-Cut (default: ON ) +#undef __TAILCUT__ +#define __TAILCUT__ + +// flag for calculation of islands stat. (default: ON ) +#define __ISLANDS__ +#undef __ISLANDS__ + +// flag for calculation of image parameters (default: ON ) +#undef __MOMENTS__ +#define __MOMENTS__ + +// flag for ROOT (default: ON ) +#undef __ROOT__ +#define __ROOT__ + +//!@} + +//=----------------------------------------------------------- +//!@subsection Definition of global variables. + +/*!@" + + Now we define some global variables with data about the telescope, + such as "focal distance", number of pixels/mirrors, + "size of the camera", and so on. + + @"*/ + +/*!@" + + Depending on the telescope we are using (CT1 or MAGIC), the + information stored in the definition file is different. + The variable |ct_Type| has the value 0 when we use + CT1, and 1 when we use MAGIC. + + @"*/ + +//!@{ +static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC +//!@} + +/*!@" + + And this is the information about the whole telescope. + + @"*/ + +//!@{ + +// parameters of the CT (from the CT definition file) + +////@: Focal distances [cm] +//static float *ct_Focal; + +//@: Mean Focal distances [cm] +static float ct_Focal_mean; + +//@: STDev. Focal distances [cm] +static float ct_Focal_std; + +//@: Mean Point Spread function [cm] +static float ct_PSpread_mean; + +//@: STDev. Point Spread function [cm] +static float ct_PSpread_std; + +//@: STDev. Adjustmente deviation [cm] +static float ct_Adjustment_std; + +//@: Radius of the Black Spot in mirror [cm] +static float ct_BlackSpot_rad; + +//@: Radius of one mirror [cm] +static float ct_RMirror; + +//@: Camera width [cm] +static float ct_CameraWidth; + +//@: Pixel width [cm] +static float ct_PixelWidth; + +//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60) +static float ct_PixelWidth_corner_2_corner; + +//@: ct_PixelWidth_corner_2_corner / 2 +static float ct_PixelWidth_corner_2_corner_half; + +//@: Number of mirrors +static int ct_NMirrors = 0; + +//@: Number of pixels +static int ct_NPixels; + +//@: Number of pixels +static int ct_NCentralPixels; + +//@: Number of pixels +static int ct_NGapPixels; + +//@: ct_Apot = ct_PixelWidth / 2 +static float ct_Apot; + +//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth +static float ct_2Apot; + +//@: name of the CT definition file to use +static char ct_filename[256]; + +//@: list of showers to be skipped +static int *Skip; + +//@: number of showers to be skipped +static int nSkip=0; + +//@: flag: TRUE: data come from STDIN; FALSE: from file +static int Data_From_STDIN = FALSE; + +//@: flag: TRUE: write all images to output; FALSE: only triggered showers +static int Write_All_Images = FALSE; + +//@: flag: TRUE: write all data to output; FALSE: only triggered showers +static int Write_All_Data = FALSE; + +//@: flag: TRUE: selection on the energy +static int Select_Energy = TRUE; + +//@: Lower edge of the selected energy range (in GeV) +static float Select_Energy_le = 0.0; + +//@: Upper edge of the selected energy range (in GeV) +static float Select_Energy_ue = 100000.0; + +//!@} + +/*!@" + + The following double-pointer is a 2-dimensional table with information + about each pixel. The routine read_pixels will generate + the information for filling it using igen_pixel_coordinates(). + + @"*/ + +//!@{ +// Pointer to a tables/Arrays with information about the pixels +// and data stored on them with information about the pixels + +//@: table for IJ sys. +static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4]; + +//@: coordinates x,y for each pixel +static float **pixary; + +//@: indexes of pixels neighbours of a given one +static int **pixneig; + +//@: number of neighbours a pixel have +static int *npixneig; + +//@: contents of the pixels (ph.e.) +static float *fnpix; + +//@: contents of the pixels (ph.e.) after cleanning +static float *fnpixclean; + +//@: contents of the sum of all ph.e. in one timeslice of 1 ns +static float *fnslicesum ; + + +//!@} + +/*!@" + + The following double-pointer is a 2-dimensional table with the + Quantum Efficiency @$QE@$ of each pixel in the camera, as a function + of the wavelength @$\lambda@$. The routine |read_pixels()| will read + also this information from the file |qe.dat|. + + @"*/ + +//!@{ +// Pointer to a table with QE, number of datapoints, and wavelengths + +//@: table of QE +static float ***QE; + +//@: number of datapoints for the QE curve +static int pointsQE; + +//@: table of QE +static float *QElambda; +//!@} + +/*!@" + + The following double-pointer is a 2-dimensional table with information + about each mirror in the dish. The routine |read_ct_file()| will read + this information from the CT definition file. + + @"*/ + +//!@{ +// Pointer to a table with the following info.: + +static float **ct_data; + +/* + * TYPE=0 (CT1) + * i s rho theta x y z thetan phin xn yn zn + * + * i : number of the mirror + * s : arc length [cm] + * rho : polar rho of the position of the center of the mirror [cm] + * theta : polar angle of the position of the center of the mirror [cm] + * x : x coordinate of the center of the mirror [cm] + * y : y coordinate of the center of the mirror [cm] + * z : z coordinate of the center of the mirror [cm] + * thetan : polar theta angle of the direction where the mirror points to + * phin : polar phi angle of the direction where the mirror points to + * xn : xn coordinate of the normal vector in the center (normalized) + * yn : yn coordinate of the normal vector in the center (normalized) + * zn : zn coordinate of the normal vector in the center (normalized) + * + * TYPE=1 (MAGIC) + * i f sx sy x y z thetan phin + * + * i : number of the mirror + * f : focal distance of that mirror + * sx : curvilinear coordinate of mirror's center in X[cm] + * sy : curvilinear coordinate of mirror's center in X[cm] + * x : x coordinate of the center of the mirror [cm] + * y : y coordinate of the center of the mirror [cm] + * z : z coordinate of the center of the mirror [cm] + * thetan : polar theta angle of the direction where the mirror points to + * phin : polar phi angle of the direction where the mirror points to + * xn : xn coordinate of the normal vector in the center (normalized) + * yn : yn coordinate of the normal vector in the center (normalized) + * zn : zn coordinate of the normal vector in the center (normalized) + */ +//!@} + +/*!@" + + We define a table into where random numbers will be stored. + The routines used for random number generation are provided by + |RANLIB| (taken from NETLIB, |www.netlib.org|), and by + the routine |double drand48(void)| (prototype defined in + |stdlib.h|) through the macro |RandomNumber| defined in + |camera.h|. + + @"*/ + +//!@{ +// table of random numbers + +// (unused) +//static double RandomNumbers[500]; +//!@} + +/*!@" + + The following is a variable to count the number of Cphotons + in the different steps of the simulation. + The definition is as follows: + @[ + \mbox{CountCphotons}[ \mbox{FILTER} ] \equiv + \mbox{\it Number of photons after the filter} \mbox{FILTER} + @] + The filters are defined and can be found in the file |camera.h|. + + @"*/ + +//!@{ +// vector to count photons at any given step of the simulation + +static int CountCphotons[10]; +//!@} + +/*!@" + + The following are the set of parameters calculated for each image. + The routines for their calculations are in |moments.cxx|. + + @"*/ + +//!@{ +// parameters of the images + +static Moments_Info *moments_ptr; +static LenWid_Info *lenwid_ptr; + +static float *maxs; +static int *nmaxs; +static float length, width, dist, xdist, azw, miss, alpha, *conc; +static float phiasym, asymx, asymy; +static float charge, smax, maxtrigthr_phe; + +//!@} + +extern char FileName[]; + + +//=----------------------------------------------------------- +// @subsection Main program. + +//!@{ + +//++++++++++++++++++++++++++++++++++++++++ +// MAIN PROGRAM +//---------------------------------------- + +int main(int argc, char **argv) +{ + + //!@' @#### Definition of variables. + //@' + + char inname[256]; //@< input file name + char outname[256]; //@< output file name + char datname[256]; //@< data (ASCII) output file name + char diagname[256]; //@< diagnistic output file (ROOT format) + char rootname[256] ; //@< ROOT file name + + char parname[256]; //@< parameters file name + + char sign[20]; //@< initialize sign + + char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file + + ifstream inputfile; //@< stream for the input file + ofstream outputfile; //@< stream for the output file + ofstream datafile; //@< stream for the data file + + MCEventHeader mcevth; //@< Event Header class (MC) + MCCphoton cphoton; //@< Cherenkov Photon class (MC) + + float thetaCT, phiCT; //@< parameters of a given shower + float thetashw, phishw; //@< parameters of a given shower + float coreD, coreX, coreY; //@< core position and distance + float impactD; //@< impact parameter + float l1, m1, n1; //@< auxiliary variables + float l2, m2, n2; //@< auxiliary variables + float num, den; //@< auxiliary variables + + int nshow=0; //@< partial number of shower in a given run + int ntshow=0; //@< total number of showers + int ncph=0; //@< partial number of shower in a given run + int ntcph=0; //@< total number of showers + + int i, ii, j, k; //@< simple counters + + float t_ini; //@< time of the first Cphoton read in + float t; //@< time for a single photon + int t_chan ; //@< the bin (channel) in time of a single photon + + int startchan ; //@< the first bin with entries in the time slices + + float cx, cy, ci, cj; //@< coordinates in the XY and IJ systems + int ici, icj, iici, iicj; //@< coordinates in the IJ (integers) + + int nPMT; //@< number of pixel + + float wl, last_wl; //@< wavelength of the photon + float qe; //@< quantum efficiency + float **qeptr; + + int simulateNSB; //@< Will we simulate NSB? + float meanNSB; //@< NSB mean value (per pixel) + float qThreshold; //@< Threshold value + float qTailCut; //@< Tail Cut value + int nIslandsCut; //@< Islands Cut value + int countIslands; //@< Will we count the islands? + int anaPixels; + + float fCorrection; //@< Factor to apply to pixel values (def. 1.) + + float q0; //@< trigger threshold ( intermediate variable ) + float maxcharge; //@< maximum charge in pixels + int noverq0, novq0; //@< number of pixels above threshold + int ngrpq0, mxgrp; //@< number of pixels in a group + + int trigger; //@< trigger flag + int itrigger; //@< index of pixel fired + int ntrigger = 0; //@< number of triggers in the whole file + unsigned char triggerBits; //@< byte for trigger condition check (MAGIC) + int bit; //@< intermediate variable + + float plateScale_cm2deg; //@< plate scale (deg/cm) + float degTriggerZone; //@< trigger area in the camera (radius, in deg.) + + float dtheta, dphi; //@< deviations of CT from shower axis + + int still_in_loop = FALSE; + + char Signature[20]; + + float *image_data; + int nvar, hidt; + + struct camera cam; // structure holding the camera definition + +#ifdef __DETAIL_TRIGGER__ + + MTrigger Trigger ; //@< A instance of the Class MTrigger + +#endif __DETAIL_TRIGGER__ + + //!@' @#### Definition of variables for |getopt()|. + //@' + + int ch, errflg = 0; //@< used by getopt + + /*!@' + + @#### Beginning of the program. + + We start with the main program. First we (could) make some + presentation, and follows the reading of the parameters file (now + from the |stdin|), the reading of the CT parameters file, and the + creation of the output file, where the processed data will be + stored. + + */ + + //++ + // START + //-- + + // make unbuffered output + + cout.setf ( ios::stdio ); + + // parse command line options (see reflector.h) + + parname[0] = '\0'; + + optarg = NULL; + while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) ) + switch (ch) { + case 'f': + strcpy(parname, optarg); + break; + case 'h': + usage(); + break; + default : + errflg++; + } + + // show help if error + + if ( errflg>0 ) + usage(); + + // make some sort of presentation + + present(); + + // read parameters file + + if ( strlen(parname) < 1 ) + readparam(NULL); + else + readparam(parname); + + // read data from file or from STDIN? + + Data_From_STDIN = get_data_from_stdin(); + + // write all images, even those without trigger? + + Write_All_Images = get_write_all_images(); + + // write all data (i.e., ph.e.s in pixels) + + Write_All_Data = get_write_all_data(); + + // get filenames + + strcpy( inname, get_input_filename() ); + strcpy( outname, get_output_filename() ); + strcpy( datname, get_data_filename() ); + strcpy( diagname, get_diag_filename() ); + strcpy( rootname, get_root_filename() ); + strcpy( ct_filename, get_ct_filename() ); + + // get different parameters of the simulation + + qThreshold = get_threshold(); + qTailCut = get_tail_cut(); + simulateNSB = get_nsb( &meanNSB ); + countIslands = get_islands_cut( &nIslandsCut ); + + // get selections on the parameters + + Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue); + + // log filenames information + + log(SIGNATURE, + "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n", + "Filenames", + "In", inname, + "Out", outname, + "Data", datname, + "Diag", diagname, + "ROOT", rootname, + "CT", ct_filename); + + + // log flags information + + log(SIGNATURE, + "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n", + "Flags", + "Data_From_STDIN", ONoff(Data_From_STDIN), + "Write_All_Images", ONoff(Write_All_Images), + "Write_All_Data", ONoff(Write_All_Data)); + + // log parameters information + + log(SIGNATURE, + "%s:\n\t%20s: %f\n\t%20s: %f\n\t%20s: %f %s\n\t%20s: %f %s\n", + "Parameters", + "q0 (Threshold)", qThreshold, + "t0 (Tail-cut)", qTailCut, + "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB), + "i0 (Islands-cut)", nIslandsCut, ONoff(countIslands)); + + // log selections + + log(SIGNATURE, + "%s:\n\t%20s: %s (%f:%f)\n", + "Selections:", + "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue); + + // set all random numbers seeds + + setall( get_seeds(0), get_seeds(1) ); + + // get list of showers to evt. skip + + nSkip = get_nskip_showers(); + + if (nSkip > 0) { + Skip = new int[ nSkip ]; + get_skip_showers( Skip ); + + log(SIGNATURE, "There are some showers to skip:\n"); + for (i=0; iSetAutoSave(100000000); + + Int_t split = 1; + Int_t bsize = 64000; + MDiagEventobject *event = 0; + + // Create one branch. If splitlevel is set, event is a superbranch + // creating a sub branch for each data member of the Eventobject event. + + tree->Branch("event", "MDiagEventobject", &event, bsize, split); + +#ifdef __ROOT__ + + MRawEvt *Evt = new MRawEvt() ; + MMcEvt *McEvt = new MMcEvt (); + + // initalize the ROOT file + // + // erzeuge ein Root file + // + + TFile outfile ( rootname , "RECREATE" ) ; + + // + // create a Tree for the Event data stream + // + + TTree EvtTree("EvtTree","Events of Run"); + + bsize=128000; split=1; + + EvtTree.Branch("MRawEvt","MRawEvt", + &Evt, bsize, split); + + EvtTree.Branch("MMcEvt","MMcEvt", + &McEvt, bsize, split); + + + float flli = 0. ; + unsigned short ulli = 0 ; + +#endif __ROOT__ + + // for safety and for dimensioning image_data: count the elements in the + // diagnostic data branch + + i=0; + i++; // "n" + i++; // "primary" + i++; // "energy" + i++; // "cored" + i++; // "impact" + i++; // "xcore" + i++; // "ycore" + i++; // "theta" + i++; // "phi" + i++; // "deviations" + i++; // "dtheta" + i++; // "dphi" + i++; // "trigger" + i++; // "ncphs" + i++; // "maxpassthr_phe" + i++; // "nphes" + i++; // "nphes2" + i++; // "length" + i++; // "width" + i++; // "dist" + i++; // "xdist" + i++; // "azw" + i++; // "miss" + i++; // "alpha" + i++; // "conc2" + i++; // "conc3" + i++; // "conc4" + i++; // "conc5" + i++; // "conc6" + i++; // "conc7" + i++; // "conc8" + i++; // "conc9" + i++; // "conc10" + i++; // "asymx" + i++; // "asymy" + i++; // "phiasym" + + nvar = i; + image_data = new float[nvar]; + + // set plate scale (deg/cm) and trigger area (deg) + + plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381; + + if ( ! get_trigger_radius( °TriggerZone ) ) + degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0); + + if ( ! get_correction( &fCorrection ) ) + fCorrection = 1.0; + + // number of pixels for parameters + + anaPixels = get_ana_pixels(); + anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels; + + // open input file if we DO read data from a file + + if (! Data_From_STDIN) { + log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname ); + inputfile.open( inname ); + if ( inputfile.bad() ) + error( SIGNATURE, "Cannot open input file: %s\n", inname ); + } + + // get signature, and check it + // NOTE: this part repeats further down in the code; + // if you change something here you probably want to change it + // there as well + + strcpy(Signature, REFL_SIGNATURE); + + strcpy(sign, Signature); + + if ( Data_From_STDIN ) + cin.read( (char *)sign, strlen(Signature)); + else + inputfile.read( (char *)sign, strlen(Signature)); + + if (strcmp(sign, Signature) != 0) { + cerr << "ERROR: Signature of .rfl file is not correct\n"; + cerr << '"' << sign << '"' << '\n'; + cerr << "should be: " << Signature << '\n'; + exit(1); + } + + if ( Data_From_STDIN ) + cin.read( (char *)sign, 1); + else + inputfile.read( (char *)sign, 1); + + // open output file + + log( SIGNATURE, "Openning output \"phe\" file %s\n", outname ); + outputfile.open( outname ); + + if ( outputfile.bad() ) + error( SIGNATURE, "Cannot open output file: %s\n", outname ); + + // open data file + + log( SIGNATURE, "Openning data \"dat\" file %s\n", datname ); + datafile.open( datname ); + + if ( outputfile.bad() ) + error( SIGNATURE, "Cannot open output file: %s\n", outname ); + + // write signature + + outputfile.write( SIGNATURE, sizeof(SIGNATURE) ); + + // initializes flag + + strcpy( flag, " \0" ); + + // allocate space for PMTs numbers of pixels + + fnpix = new float [ ct_NPixels ]; + fnpixclean = new float [ ct_NPixels ]; + +#ifdef __ROOT__ + + fnslicesum = new float [ (2 * SLICES) ] ; + + float slices [ct_NPixels][ (2 * SLICES) ] ; + float slices2 [ct_NPixels][ SLICES ] ; + + float trans [ SLICES ] ; +#endif __ROOT__ + + + moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 ); + + //!@' @#### Main loop. + //@' + + //begin my version + + // get flag + + if ( Data_From_STDIN ) + cin.read( flag, SIZE_OF_FLAGS ); + else + inputfile.read ( flag, SIZE_OF_FLAGS ); + + // loop over the file + + still_in_loop = TRUE; + + while ( + ((! Data_From_STDIN) && (! inputfile.eof())) + || + (Data_From_STDIN && still_in_loop) + ) { + + // reading .rfl files + if(!isA( flag, FLAG_START_OF_RUN )){ + error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag ); + } + else { // found start of run + nshow=0; + // read next flag + + if ( Data_From_STDIN ) + cin.read( flag, SIZE_OF_FLAGS ); + else + inputfile.read ( flag, SIZE_OF_FLAGS ); + + while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event + /*!@' + + For the case |FLAG\_START\_OF\_EVENT|, + we read each Cherenkov photon, and follow these steps: + + @enumerate + + @- Transform XY-coordinates to IJ-coordinates. + + @- With this, we obtain the pixel where the photon hits. + + @- Use the wavelength $\lambda$ and the table of QE, and + calculate the estimated (third order interpolated) quantum + efficiency for that photon. The photon can be rejected. + + @- If accepted, then add to the pixel. + + @endenumerate + + In principle, we should stop here, and use another program to + 'smear' the image, to add the Night Sky Background, and to + simulate the trigger logic, but we will make this program + quick and dirty, and include all here. + + If we are reading PHE files, we jump to the point where the + pixelization process already has finished. + + */ + + ++nshow; + log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow); + + // get MCEventHeader + + if ( Data_From_STDIN ) + cin.read( (char*)&mcevth, mcevth.mysize() ); + else + mcevth.read( inputfile ); + + // calculate core distance and impact parameter + + coreD = mcevth.get_core(&coreX, &coreY); + + // calculate impact parameter (shortest distance betwee the original + // trajectory of the primary (assumed shower-axis) and the + // direction where the telescope points to + // + // we use the following equation, given that the shower core position + // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1), + // and the telescope position and orientation are (x2,y2,z2)=(0,0,0) + // and (l2,m2,n2) + // + // | | + // | x1-x2 y1-y2 z1-z2 | + // | | + // + | l1 m1 n1 | + // - | | + // | l2 m2 n2 | + // | | + // dist = ------------------------------------ ( > 0 ) + // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2 + // [ | | + | | + | | ] + // [ |l2 m2| |m2 n2| |n2 l2| ] + // + // playing a little bit, we get this reduced for in our case: + // + // + // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) / + // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 - + // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2) + + // read the direction of the incoming shower + + thetashw = mcevth.get_theta(); + phishw = mcevth.get_phi(); + + // calculate vector for shower + + l1 = sin(thetashw)*cos(phishw); + m1 = sin(thetashw)*sin(phishw); + n1 = cos(thetashw); + + // read the deviation of the telescope with respect to the shower + + mcevth.get_deviations ( &thetaCT, &phiCT ); + + 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 + + impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. ); + + } 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; + + // fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ", + // thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2, + // num, den); + + } + + // clear camera + + for ( i=0; i (2 * SLICES)){ + log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n", + t_chan, (2 * SLICES), (2 * SLICES)); + t_chan = (2 * SLICES); + } + else if(t_chan < 0){ + log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n", + t_chan, 0, 0); + t_chan = 0; + } + /*!@' + + @#### Pixelization (for the central pixels). + + In order to calculate the coordinates, we use the + change of system described in the documentation + of the source code of |pixel\_coord.cxx|. + Then, we will use simply the matrix of change + from one system to the other. In our case, this is: + + @[ + \begin{bmatrix}X\\Y\\\end{bmatrix} + = + \begin{bmatrix} + 1 & \cos(60^\circ)\\ + 0 & \sin(60^\circ)\\ + \end{bmatrix} + \begin{bmatrix}I\\J\\\end{bmatrix} + @] + + and hence + + @[ + \begin{bmatrix}I\\J\\\end{bmatrix} + = + \begin{bmatrix} + 1 & -\frac{\cos(60^\circ)}{\sin(60^\circ)}\\ + 0 &\frac{1}{\sin(60^\circ)}\\ + \end{bmatrix} + \begin{bmatrix}X\\Y\\\end{bmatrix} + @] + + */ + + //+++ + // Pixelization + //--- + + // calculate ij-coordinates + + // We use a change of coordinate system, using the following + // matrix of change (m^-1) (this is taken from Mathematica output). + /* + * In[1]:= m={{1,cos60},{0,sin60}}; MatrixForm[m] + * + * Out[1]//MatrixForm= 1 cos60 + * + * 0 sin60 + * + * In[2]:= inv=Inverse[m]; MatrixForm[inv] + * + * Out[2]//MatrixForm= cos60 + * -(-----) + * 1 sin60 + * + * 1 + * ----- + * 0 sin60 + * + */ + + // go to IJ-coordinate system + + cx = cphoton.get_x(); + cy = cphoton.get_y(); + + // get wavelength + + last_wl = wl; + wl = cphoton.get_wl(); + + if ( wl < 1.0 ) + break; + + if ( (wl > 600.0) || (wl < 290.0) ) + break; + + // check if photon is inside outermost camera radius + + if(sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)){ + + // read next CPhoton + if ( Data_From_STDIN ) + cin.read( flag, SIZE_OF_FLAGS ); + else + inputfile.read ( flag, SIZE_OF_FLAGS ); + + // go to beginning of loop, the photon is lost + continue; + + } + + // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl; + + ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5); + cj = floor( (cy/SIN60) / ct_2Apot + 0.5); + + ici = (int)(ci); + icj = (int)(cj); + + iici = ici+PIX_ARRAY_HALF_SIDE; + iicj = icj+PIX_ARRAY_HALF_SIDE; + + // is it inside the array? + + if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) && + (iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) { + + // try to put into pixel + + // obtain the pixel number for this photon + + nPMT = (int) + pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM]; + + } + else{ + + nPMT = -1; + + } + + // check if outside the central camera + + if ( (nPMT < 0) || (nPMT >= ct_NCentralPixels) ) { + + // check the outer pixels + nPMT = -1; + + for(i=ct_NCentralPixels; i k) + + qeptr = (float **)QE[nPMT]; + + FindLagrange(qeptr,k,wl); + + // if random > quantum efficiency, reject it + + qe = Lagrange(qeptr,k,wl) / 100.0; + + // fprintf(stdout, "%f\n", qe); + + if ( RandomNumber > qe ) { + + // read next CPhoton + if ( Data_From_STDIN ) + cin.read( flag, SIZE_OF_FLAGS ); + else + inputfile.read ( flag, SIZE_OF_FLAGS ); + + // go to beginning of loop + continue; + + } + +#endif // __QE__ + + //+++ + // Cphoton is accepted + //--- + + // increase the number of Cphs. in the PMT, i.e., + // increase in one unit the counter of the photons + // stored in the pixel nPMT + + fnpix[nPMT] += 1.0; + +#ifdef __ROOT__ + fnslicesum[t_chan] += 1.0 ; + slices[nPMT][t_chan] += 1.0 ; +#endif __ROOT__ + +#ifdef __DETAIL_TRIGGER__ + // + // fill the Trigger class with this phe + // + // + + Trigger.Fill( nPMT, ( t - t_ini ) ) ; +#endif __DETAIL_TRIGGER__ + + // read next CPhoton + + if ( Data_From_STDIN ) + cin.read( flag, SIZE_OF_FLAGS ); + else + inputfile.read ( flag, SIZE_OF_FLAGS ); + + } // end while, i.e. found end of event + + log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n", + ncph, ntcph); + + // show number of photons + + //cout << ncph << " photons read . . . " << endl << flush; + + // skip it ? + + for ( i=0; i Select_Energy_ue )) { + log(SIGNATURE, "select_energy: shower rejected.\n"); + continue; + } + } + +#ifdef __NSB__ + + //!@' @#### NSB (Night Sky Background) simulation. + //@' + + //+++ + // NSB simulation + //--- + + // add NSB "noise" + // TO DO: make meanNSB an array and read the contents from a file! + + if ( simulateNSB ) + for ( i=0; iFillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ; + + // now put out the data of interest + // + // 1. -> look for the first slice with signal + // + + for ( i=0; i<(2 * SLICES); i++ ) + if ( fnslicesum[i] > 0. ) + break ; + + startchan = i ; + + // + // copy the slices out of the big array + // + // put the first slice with signal to slice 4 + // + + for (i=0; i 0 ) { + + for ( ii=0; ii < 15; ii++ ) { + trans [ii] = slices2[i][ii] ; + } + + Evt->FillPixel ( (UShort_t) i , trans ) ; + + } + } + + // + // + // + + McEvt->Fill( (UShort_t) mcevth.get_primary() , + mcevth.get_energy(), + mcevth.get_theta(), + mcevth.get_phi(), + mcevth.get_core(), + mcevth.get_coreX(), + mcevth.get_coreY(), + flli, + ulli, ulli, ulli, ulli, ulli ) ; + + // + // write it out to the file outfile + // + + EvtTree.Fill() ; + + // clear all + + Evt->Clear() ; + McEvt->Clear() ; + +#endif __ROOT__ + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // at this point we have a camera full of + // ph.e.s + // we should first apply the trigger condition, + // and if there's trigger, then clean the image, + // calculate the islands statistics and the + // other parameters of the image (Hillas' parameters + // and so on). + //-------------------------------------------------- + +#ifdef __DEBUG__ + printf("\n"); + + for ( ici=0; ici -1 ) { + + if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) { + + printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow, + (int)pixels[ici][icj][PIXNUM], + pixels[ici][icj][PIXX], + pixels[ici][icj][PIXY], + fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj); + + } + + } + + } + + } + + for (i=0; i degTriggerZone) + continue; + + // 'ngrpq0' is the number of neighbours of pixel i with q > q0 + + ngrpq0 = 0; + + // look at each pixel in the neighborhood, and count + // those above threshold q0 + + for ( j=0 ; j-1 ; ++j ) + if ( fnpix[pixneig[i][j]] > q0 ) + ++ngrpq0; + + // check whether we have trigger + + if ( ct_Type == 0 ) { + + //++ >>>>> CT1 <<<<< + +#ifdef __CT1_NO_NEIGHBOURS__ + + if ( noverq0 > 1 ) + trigger = TRUE; + +#else + + if ( ngrpq0 > 0 ) + trigger = TRUE; + +#endif + + //-- >>>>> CT1 <<<<< + + } else { + + //++ >>>>> MAGIC <<<<< + + // (at least 4 packed with at least q0 phes) + + // there are 3 cases + // 1. zero, one or two neighbours have enough charge: no trigger + // 2. five or six neighbours with enough charge: trigger? sure!! + // 3. three or four neighbours with q > q0 : we must look + // for 'closeness'. + + switch ( ngrpq0 ) { + + case 0: + case 1: + case 2: + + trigger = FALSE; + break; + + case 3: + case 4: + + // if reaches this line, it means three or four neighbours + // around the central pixel + + triggerBits = 1; + + for ( j=0 ; j-1; ++j ) { + + if ( fnpix[pixneig[i][j]] > q0 ) { + + if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) { + + if ( nint(pixary[pixneig[i][j]][1]*10.0) > + nint(pixary[i][1]*10.0) ) + bit = 2; + else if ( nint(pixary[pixneig[i][j]][1]*10.0) < + nint(pixary[i][1]*10.0) ) + bit = 6; + else + bit = 1; + + } else { + + if ( nint(pixary[pixneig[i][j]][1]*10.0) > + nint(pixary[i][1]*10.0) ) + bit = 3; + else if ( nint(pixary[pixneig[i][j]][1]*10.0) < + nint(pixary[i][1]*10.0) ) + bit = 5; + else + bit = 4; + + } + + triggerBits |= (1< 6 !!! Exiting.\n\n"); + break; + + } // switch (ngrpq0) + + } // ct_Type + + } // for each pixel i + + if ( trigger == FALSE ) { + break; + } // end if + + } //end for each threshold + maxtrigthr_phe = (float) k-1; // i.e. maxtrigthr_phe < 0. means that + // the event doesn't even pass threshold 0. + // maxtrigthr_phe >= 0 means, the event passes some threshold + // or (in case the input parameter "threshold" was > 0), the event + // passes the threshold given by the input parameter. + if ( maxtrigthr_phe >= 0. ) { + trigger = TRUE; + } + + + novq0 = noverq0; + + if ( trigger == TRUE ) { + + itrigger = i; + ++ntrigger; + + memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels ); + +#ifdef __TAILCUT__ + + //!@' @#### Tail-cut condition. + //@' + + //++ + // tail-cut + //-- + + // Tail-Cut = 0 : No Tail-Cut + // Tail-Cut > 0 : Make Tail-Cut + // Tail-Cut < 0 : Make Tail-Cut with t_0 = Sqrt[ maximum ] + + if (qTailCut > 0.0) { + + for ( i=0; icharge ; + smax = moments_ptr->smax ; + maxs = moments_ptr->maxs ; + nmaxs = moments_ptr->nmaxs ; + length = moments_ptr->length ; + width = moments_ptr->width ; + dist = moments_ptr->dist ; + xdist = moments_ptr->xdist ; + azw = moments_ptr->azw ; + miss = moments_ptr->miss ; + alpha = moments_ptr->alpha ; + conc = moments_ptr->conc ; + asymx = moments_ptr->asymx ; + asymx = moments_ptr->asymx ; + phiasym= moments_ptr->phi; + + lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary, + plateScale_cm2deg, + ct_PixelWidth_corner_2_corner_half); + + + // fill the diagnostic Tree + + event = new MDiagEventobject(); + + i=0; + image_data[i] = event->n = hidt/10; i++; + image_data[i] = event->primary = mcevth.get_primary(); i++; + image_data[i] = event->energy = mcevth.get_energy(); i++; + image_data[i] = event->cored = coreD; i++; + image_data[i] = event->impact = impactD; i++; + image_data[i] = event->xcore = coreX; i++; + image_data[i] = event->ycore = coreY; i++; + image_data[i] = event->theta = mcevth.get_theta(); i++; + image_data[i] = event->phi = mcevth.get_phi(); i++; + image_data[i] = event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++; + image_data[i] = event->dtheta = dtheta; i++; + image_data[i] = event->dphi = dphi; i++; + image_data[i] = event->trigger = trigger; i++; + image_data[i] = event->ncphs = ncph; i++; + image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++; + image_data[i] = event->nphes = charge; i++; + image_data[i] = event->nphes2 = smax; i++; + image_data[i] = event->length = length; i++; + image_data[i] = event->width = width; i++; + image_data[i] = event->dist = dist; i++; + image_data[i] = event->xdist = xdist; i++; + image_data[i] = event->azw = azw; i++; + image_data[i] = event->miss = miss; i++; + image_data[i] = event->alpha = alpha; i++; + image_data[i] = event->conc2 = conc[0]; i++; + image_data[i] = event->conc3 = conc[1]; i++; + image_data[i] = event->conc4 = conc[2]; i++; + image_data[i] = event->conc5 = conc[3]; i++; + image_data[i] = event->conc6 = conc[4]; i++; + image_data[i] = event->conc7 = conc[5]; i++; + image_data[i] = event->conc8 = conc[6]; i++; + image_data[i] = event->conc9 = conc[7]; i++; + image_data[i] = event->conc10 = conc[8]; i++; + image_data[i] = event->asymx = asymx; i++; + image_data[i] = event->asymy = asymy; i++; + image_data[i] = event->phiasym = phiasym; i++; + + // there should be "nvar" variables + + if ( i != nvar ) + error( SIGNATURE, "Wrong entry number for diagnostic data.\n" ); + + tree->Fill(); + delete event; + + // put information in the data file, + + datafile << ntrigger; + for(i=0;in = hidt/10; i++; + image_data[i] = event->primary = mcevth.get_primary(); i++; + image_data[i] = event->energy = mcevth.get_energy(); i++; + image_data[i] = event->cored = coreD = mcevth.get_core(&coreX, &coreY); i++; + image_data[i] = event->impact = coreD; i++; + image_data[i] = event->xcore = coreX; i++; + image_data[i] = event->ycore = coreY; i++; + image_data[i] = event->theta = mcevth.get_theta(); i++; + image_data[i] = event->phi = mcevth.get_phi(); i++; + image_data[i] = event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++; + image_data[i] = event->dtheta = dtheta; i++; + image_data[i] = event->dphi = dphi; i++; + image_data[i] = event->trigger = trigger; i++; + image_data[i] = event->ncphs = ncph; i++; + image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + image_data[i] = -1.; i++; + + // there should be "nvar" variables + + if ( i != nvar ) + error( SIGNATURE, "Wrong entry length for Ntuple.\n" ); + + tree->Fill(); + delete event; + + // put this information in the data file, + + if ( Write_All_Data ) { + + datafile << ntrigger; + for ( i=0; iWrite(); + + hfile->Close(); + +#ifdef __DETAIL_TRIGGER__ + // Output of Trigger statistics + // + + Trigger.PrintStat() ; +#endif __DETAIL_TRIGGER__ + + // program finished + + log( SIGNATURE, "Done.\n"); + + return( 0 ); +} +//!@} + +// @T \newpage + +//!@subsection Functions definition. + +//!----------------------------------------------------------- +// @name present +// +// @desc Make some presentation +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +present(void) +{ + cout << "##################################################\n" + << SIGNATURE << '\n' << '\n' + << "Processor of the reflector output\n" + << "J C Gonzalez, Jun 1998\n" + << "##################################################\n\n" + << flush ; +} +//!@} + + +//!----------------------------------------------------------- +// @name usage +// +// @desc show help +// +// @date Tue Dec 15 16:23:30 MET 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +usage(void) +{ + present(); + cout << "\nusage ::\n\n" + << "\t camera " + << " [ -@ paramfile ] " + << " [ -h ] " + << "\n\n or \n\n" + << "\t camera < paramfile" + << "\n\n"; + exit(0); +} +//!@} + + +//!----------------------------------------------------------- +// @name log +// +// @desc function to send log information +// +// @var funct Name of the caller function +// @var fmt Format to be used (message) +// @var ... Other information to be shown +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +log(const char *funct, char *fmt, ...) +{ + va_list args; + + // Display the name of the function that called error + printf("[%s]: ", funct); + + // Display the remainder of the message + va_start(args, fmt); + vprintf(fmt, args); + va_end(args); +} +//!@} + + +//!----------------------------------------------------------- +// @name error +// +// @desc function to send an error message, and abort the program +// +// @var funct Name of the caller function +// @var fmt Format to be used (message) +// @var ... Other information to be shown +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +error(const char *funct, char *fmt, ...) +{ + va_list args; + + // Display the name of the function that called error + fprintf(stderr, "ERROR in %s: ", funct); + + // Display the remainder of the message + va_start(args, fmt); + vfprintf(stderr, fmt, args); + va_end(args); + + perror(funct); + + abort(); +} +//!@} + + +//!----------------------------------------------------------- +// @name isA +// +// @desc returns TRUE(FALSE), if the flag is(is not) the given +// +// @var s1 String to be searched +// @var flag Flag to compare with string s1 +// @return TRUE: both strings match; FALSE: oth. +// +// @date Wed Jul 8 15:25:39 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +int +isA( char * s1, const char * flag ) { + return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 ); +} +//!@} + + +//!----------------------------------------------------------- +// @name read_ct_file +// +// @desc read CT definition file +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +read_ct_file(void) +{ + char line[LINE_MAX_LENGTH]; //@< line to get from the ctin + char token[ITEM_MAX_LENGTH]; //@< a single token + int i, j; //@< dummy counters + + log( "read_ct_file", "start.\n" ); + + ifstream ctin ( ct_filename ); + + if ( ctin.bad() ) + error( "read_ct_file", + "Cannot open CT def. file: %s\n", ct_filename ); + + // loop till the "end" directive is reached + + while (!ctin.eof()) { + + // get line from stdin + + ctin.getline(line, LINE_MAX_LENGTH); + + // look for each item at the beginning of the line + + for (i=0; i<=define_mirrors; i++) + if (strstr(line, CT_ITEM_NAMES[i]) == line) + break; + + // if it is not a valid line, just ignore it + + if (i == define_mirrors+1) + continue; + + // case block for each directive + + switch ( i ) { + + case type: // (0:CT1 ¦ 1:MAGIC) + + // get focal distance + + sscanf(line, "%s %d", token, &ct_Type); + + log( "read_ct_file", ": %s\n", + ((ct_Type==0) ? "CT1" : "MAGIC") ); + + break; + + case focal_distance: // [cm] + + // get focal distance + + sscanf(line, "%s %f", token, &ct_Focal_mean); + + log( "read_ct_file", ": %f cm\n", ct_Focal_mean ); + + break; + + case focal_std: // s(focal distance) [cm] + + // get focal distance + + sscanf(line, "%s %f", token, &ct_Focal_std); + + log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std ); + + break; + + case point_spread: // [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_PSpread_mean); + + log( "read_ct_file", ": %f cm\n", ct_PSpread_mean ); + + break; + + case point_std: // s(point spread) [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_PSpread_std); + + log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std ); + + break; + + case adjustment_dev: // s(adjustment_dev) [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_Adjustment_std); + + log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std ); + + break; + + case black_spot: // radius of the black spot in the center [cm] + + // get black spot radius + + sscanf(line, "%s %f", token, &ct_BlackSpot_rad); + + log( "read_ct_file", "Radius of the black spots: %f cm\n", + ct_BlackSpot_rad); + + break; + + case r_mirror: // radius of the mirrors [cm] + + // get radius of mirror + + sscanf(line, "%s %f", token, &ct_RMirror); + + log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror ); + + break; + + case n_mirrors: // number of mirrors + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NMirrors); + + log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors ); + + break; + + case camera_width: // camera width [cm] + + // get the name of the ct_file from the line + + sscanf(line, "%s %f", token, &ct_CameraWidth); + + log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth ); + + break; + + case n_pixels: // number of pixels + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NPixels); + + log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels ); + + break; + + case n_centralpixels: // number of central pixels + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NCentralPixels); + + log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels ); + + break; + + case n_gappixels: // number of gap pixels + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NGapPixels); + + log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels ); + + break; + + case pixel_width: // pixel width [cm] + + // get the name of the ct_file from the line + + sscanf(line, "%s %f", token, &ct_PixelWidth); + + ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0)); + ct_PixelWidth_corner_2_corner_half = + ct_PixelWidth_corner_2_corner * 0.50; + ct_Apot = ct_PixelWidth / 2; + ct_2Apot = ct_Apot * 2.0; + + log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth ); + + break; + + case define_mirrors: // read table with the parameters of the mirrors + + log( "read_ct_file", "Table of mirrors data:\n" ); + + // check whether the number of mirrors was already set + + if ( ct_NMirrors == 0 ) + error( "read_ct_file", "NMirrors was not set.\n" ); + + // allocate memory for paths list + + log( "read_ct_file", "Allocating memory for ct_data\n" ); + + ct_data = new float*[ct_NMirrors]; + + for (i=0; i> ct_data[i][j]; + + break; + + } // switch ( i ) + + } // end while + + // end + + log( "read_ct_file", "done.\n" ); + + return; +} +//!@} + + +//!----------------------------------------------------------- +// @name read_pixels +// +// @desc read pixels data +// +// @date Fri Mar 12 16:33:34 MET 1999 +//------------------------------------------------------------ +// @function + +//!@{ +void +read_pixels(struct camera *pcam) +{ + ifstream qefile; + char line[LINE_MAX_LENGTH]; + int n, i, j, k; + float qe; + + //------------------------------------------------------------ + // first, pixels' coordinates + + pcam->inumpixels = ct_NPixels; + pcam->inumcentralpixels = ct_NCentralPixels; + pcam->inumgappixels = ct_NGapPixels; + pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels; + pcam->dpixdiameter_cm = ct_PixelWidth; + + // initialize pixel numbers + + for ( i=0; idi[k]; + j = (int) pcam->dj[k]; + + pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k; + pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k]; + pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k]; + + pixary[k][0] = pcam->dxc[k]; + pixary[k][1] = pcam->dyc[k]; + } + + // calculate tables of neighbours + +#ifdef __DEBUG__ + for ( n=0 ; n -1) && + ((j-1) < pointsQE) && ((j-1) > -1) ) { + QE[i-1][0][j-1] = QElambda[j-1]; + QE[i-1][1][j-1] = qe; + } + + } + + // close file + + qefile.close(); + + // end + + log("read_pixels", "Done.\n"); + +} +//!@} + + +//!----------------------------------------------------------- +// @name pixels_are_neig +// +// @desc check whether two pixels are neighbours +// +// @var pix1 Number of the first pixel +// @var pix2 Number of the second pixel +// @return TRUE: both pixels are neighbours; FALSE: oth. +// +// @date Wed Sep 9 17:58:37 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +int +pixels_are_neig(int pix1, int pix2) +{ + if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) + + SQR( pixary[pix1][1] - pixary[pix2][1] ) ) + > ct_PixelWidth_corner_2_corner ) + return ( FALSE ); + else + return ( TRUE ); +} +//!@} + +//!----------------------------------------------------------- +// @name igen_pixel_coordinates +// +// @desc generate the pixel center coordinates +// +// @var *pcam structure camera containing all the +// camera information +// @return total number of pixels +// +// DP +// +// @date Thu Oct 14 10:41:03 CEST 1999 +//------------------------------------------------------------ +// @function + +//!@{ +/******** igen_pixel_coordinates() *********************************/ + +int igen_pixel_coordinates(struct camera *pcam) { + /* generate pixel coordinates, return value is number of pixels */ + + int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment; + float fsegment_fract; + double dtsize; + double dhsize; + double dpsize; + double dxfirst_pix; + double dyfirst_pix; + double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6; + double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6; + + + double dstartx, dstarty; /* for the gap pixels and outer pixels */ + int j, nrow; + + dpsize = pcam->dpixdiameter_cm; + dtsize = dpsize * sqrt(3.) / 2.; + dhsize = dpsize / 2.; + + /* Loop over central pixels to generate co-ordinates */ + + for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){ + + /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */ + + pcam->dpixsizefactor[ipixno] = 1.; + + in = 0; + + i = 0; + itot_inside_ring = 0; + iring_no = 0; + + /* Calculate the number of pixels out to and including the ring containing pixel number */ + /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */ + + while (itot_inside_ring == 0){ + + iN = 3*(i*(i+1)) + 1; + + if (ipixno <= iN){ + iring_no = i; + itot_inside_ring = iN; + } + + i++; + } + + + /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */ + + ipix_in_ring = 0; + for (i = 0; i < iring_no; ++i){ + + ipix_in_ring = ipix_in_ring + 6; + } + + /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */ + /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */ + /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */ + + isegment = 0; + fsegment_fract = 0.; + if (iring_no > 0) { + + isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */ + + fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ; + + } + + /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */ + /* pixel width (flat to flat) dpsize. */ + + dxfirst_pix = dpsize*iring_no; + dyfirst_pix = 0.; + + /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */ + + ddxseg1 = - dhsize*iring_no; + ddyseg1 = dtsize*iring_no; + ddxseg2 = -dpsize*iring_no; + ddyseg2 = 0.; + ddxseg3 = ddxseg1; + ddyseg3 = -ddyseg1; + ddxseg4 = -ddxseg1; + ddyseg4 = -ddyseg1; + ddxseg5 = -ddxseg2; + ddyseg5 = 0.; + ddxseg6 = -ddxseg1; + ddyseg6 = ddyseg1; + + /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */ + /* anti-clockwise around the ring by adding the segment to segment vectors. */ + + switch (isegment) { + + case 0: + + pcam->dxc[ipixno-1] = 0.; + pcam->dyc[ipixno-1] = 0.; + + case 1: + pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract; + + break; + + case 2: + + pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.; + + break; + + case 3: + + pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract; + + break; + + case 4: + + pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract; + + break; + + case 5: + + pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.; + + break; + + case 6: + + pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract; + pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract; + + break; + + default: + + fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno); + return(0); + + } /* end switch */ + + } /* end for */ + + dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize; + dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize; + + if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */ + + j = pcam->inumcentralpixels; + + for(i=0; iinumgappixels; i=i+6){ + pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize; + pcam->dyc[j + i ] = dstarty; + pcam->dpixsizefactor[j + i] = 1.; + pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.; + pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1]; + pcam->dpixsizefactor[j + i + 1] = 1.; + pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1]; + pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1]; + pcam->dpixsizefactor[j + i+ 2] = 1.; + pcam->dxc[j + i + 3] = - pcam->dxc[j + i]; + pcam->dyc[j + i + 3] = dstarty; + pcam->dpixsizefactor[j + i+ 3] = 1.; + pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2]; + pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2]; + pcam->dpixsizefactor[j + i+ 4] = 1.; + pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1]; + pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1]; + pcam->dpixsizefactor[j + i + 5] = 1.; + } /* end for */ + } /* end if */ + + /* generate positions of the outer pixels */ + + if( pcam->inumbigpixels > 0 ){ + + j = pcam->inumcentralpixels + pcam->inumgappixels; + + for(i=0; iinumbigpixels; i++){ + pcam->dpixsizefactor[j + i] = 2.; + } + + in = 0; + + nrow = (int) ceil(dstartx / 2. / dpsize); + + while(in < pcam->inumbigpixels){ + + pcam->dxc[j + in] = dstartx + dpsize; + pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.); + pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.; + pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.); + pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in]; + pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in]; + pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in]; + pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in]; + pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow]; + pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow]; + pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in]; + pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in]; + for(i=1; idxc[j + in + i] = pcam->dxc[j + in] - i * dpsize; + pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.); + pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize; + pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow]; + pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i]; + pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i]; + pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i]; + pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i]; + pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow]; + pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow]; + pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i]; + pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i]; + } + in = in + 6 * nrow; + dstartx = dstartx + 2. * dpsize; + nrow = nrow + 1; + } /* end while */ + + } /* end if */ + + /* generate the ij coordinates */ + + for(i=0; iinumpixels; i++){ + pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize; + pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60; + + // fprintf(stdout, "%d %f %f %f %f %f\n", + // i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i], + // pcam->dpixsizefactor[i]); + + } + + return(pcam->inumpixels); + +} +//!@} + +//!----------------------------------------------------------- +// @name bpoint_is_in_pix +// +// @desc check if a point (x,y) in camera coordinates is inside a given pixel +// +// @var *pcam structure camera containing all the +// camera information +// @var dx, dy point coordinates in centimeters +// @var ipixnum pixel number (starting at 0) +// @return TRUE if the point is inside the pixel, FALSE otherwise +// +// DP +// +// @date Thu Oct 14 16:59:04 CEST 1999 +//------------------------------------------------------------ +// @function + +//!@{ + +/******** 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 */ + + b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum]; + a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum]; + c = 1. - 1./sqrt(3.); + if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){ + fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum); + fprintf(stderr, "Exiting.\n"); + exit(203); + } + xx = dx - pcam->dxc[ipixnum]; + yy = dy - pcam->dyc[ipixnum]; + + 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 */ + } + else{ + return(FALSE); /* point is outside */ + } +} + +//!@} + +//------------------------------------------------------------ +// @name dist_r_P +// +// @desc distance straight line r - point P +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +// @function @code +//------------------------------------------------------------ +// dist_r_P +// +// distance straight line r - point P +//------------------------------------------------------------ + +float +dist_r_P(float a, float b, float c, + float l, float m, float n, + float x, float y, float z) +{ + return ( + sqrt((SQR((a-x)*m-(b-y)*l) + + SQR((b-y)*n-(c-z)*m) + + SQR((c-z)*l-(a-x)*n))/ + (SQR(l)+SQR(m)+SQR(n)) + ) + ); +} +// @endcode + + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ +// +// $Log: not supported by cvs2svn $ +// Revision 1.3 1999/10/22 15:01:28 petry +// version sent to H.K. and N.M. on Fri Oct 22 1999 +// +// Revision 1.2 1999/10/22 09:44:23 petry +// first synthesized version which compiles and runs without crashing; +// +// Revision 1.1.1.1 1999/10/21 16:35:10 petry +// first synthesised version +// +// Revision 1.13 1999/03/15 14:59:05 gonzalez +// camera-1_1 +// +// Revision 1.12 1999/03/02 09:56:10 gonzalez +// *** empty log message *** +// +// +//!@} + +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.h (revision 308) @@ -0,0 +1,302 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= camera +//= +//= @file camera.h +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: camera.h,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-11-05 11:59:31 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |camera.h|. + +/*!@" + + This section shows the include file |camera.h| + + @"*/ + +//=----------------------------------------------------------- +//!@subsection Include files. + +//!@{ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "camera-v.h" + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "creadparam.h" +#include "atm.h" +#include "moments.h" + +#include "lagrange.h" + +#include "MCEventHeader.hxx" +#include "MCCphoton.hxx" + +// command line options available + +#define COMMAND_LINE_OPTIONS "f:h" + +/*@' + + This is C++, but RANLIB routines are written in pure ANSI C. + In order to read easily these routines, we must include + the following directive + +*/ + +extern "C" { +#include "ranlib.h" +} + +// version of the reflector program that can read + +#define REFL_PROGRAM reflector +#define REFL_VERSION 0.4 + +const char REFL_SIGNATURE[] = GLUE_postp( REFL_PROGRAM, REFL_VERSION ); + +//!@} + +//=------------------------------------------------------------ +//!@subsection Macro-definitions, and constants. + +//!@{ + +#define SLICES 100 +#define WIDTH_TIMESLICE 3.3 + +#define SIN60 0.866025403784439 +#define COS60 0.500000000000000 + +#define RandomNumber drand48() + +#define PIX_ARRAY_SIDE 40 +#define PIX_ARRAY_HALF_SIDE 20 +#define PIXNUM 0 +#define PIXX 1 +#define PIXY 2 + +#define iMAXNUMPIX 595 // total maximum possible number of pixels in the camera + +//@ the trigger threshold up to which the maximum passable threshold is tested +#define iMAX_THRESHOLD_PHE 50 + +//@ number of the 1st. pixel of a sector s in a ring r (central pixel: ring=0) +#define FIRST_PIXEL(r,s) ( ((r)>0) ? (3*(r)*((r)-1) + (r)*(s) + 1) : 0 ) + +//@ number of the pixels include in a camera of r pixels +#define NUMBER_PIXELS(r) ( ((r)>0) ? FIRST_PIXEL((r)+1,0) : 1 ) + +//@ now we define the list CT_ITEM_LIST of possible items in the CT def. file +#define CT_ITEM_LIST /* LIST OF ITEMS IN THE CT DEFINITION FILE */ \ +T(type), /* type of definition file */ \ +T(focal_distance), /* std(focal distance) */ \ +T(focal_std), /* focal distance */ \ +T(point_spread), /* std(point spread) */ \ +T(point_std), /* point spread */ \ +T(adjustment_dev), /* std of adjustment deviation */ \ +T(black_spot), /* radius of the black spot in center of mirror */ \ +T(n_mirrors), /* number of mirrors */ \ +T(r_mirror), /* radius of one mirror */ \ +T(camera_width), /* camera width */ \ +T(n_pixels), /* total number of pixels in the camera */ \ +T(n_centralpixels), /* number of central pixels in the camera */ \ +T(n_gappixels), /* number of gap pixels in the camera */ \ +T(pixel_width), /* pixel width */ \ +T(define_mirrors) /* this entry is followed by the def. of pixels */ + +#define T(x) x //@< define T() as the name as it is + +enum CT_ITEM_TYPE { + CT_ITEM_LIST +}; + +#undef T + +#define T(x) #x //@< define T() as the string of x + +const char *const CT_ITEM_NAMES[] = { + CT_ITEM_LIST +}; + +#undef T + + +// TYPE=0 (CT1) +// i s rho theta x y z thetan phin xn yn zn +// +// i : number of the mirror +// s : arc length [cm] +// rho : polar rho of the position of the center of the mirror [cm] +// theta : polar angle of the position of the center of the mirror [cm] +// x : x coordinate of the center of the mirror [cm] +// y : y coordinate of the center of the mirror [cm] +// z : z coordinate of the center of the mirror [cm] +// thetan : polar theta angle of the direction where the mirror points to +// phin : polar phi angle of the direction where the mirror points to +// xn : xn coordinate of the normal vector in the center (normalized) +// yn : yn coordinate of the normal vector in the center (normalized) +// zn : zn coordinate of the normal vector in the center (normalized) +// +// TYPE=1 (MAGIC) +// i f sx sy x y z thetan phin +// +// i : number of the mirror +// f : focal distance of that mirror +// sx : curvilinear coordinate of mirror's center in X[cm] +// sy : curvilinear coordinate of mirror's center in X[cm] +// x : x coordinate of the center of the mirror [cm] +// y : y coordinate of the center of the mirror [cm] +// z : z coordinate of the center of the mirror [cm] +// thetan : polar theta angle of the direction where the mirror points to +// phin : polar phi angle of the direction where the mirror points to +// xn : xn coordinate of the normal vector in the center (normalized) +// yn : yn coordinate of the normal vector in the center (normalized) +// zn : zn coordinate of the normal vector in the center (normalized) + +#define CT_I 0 + +#define CT_S 1 +#define CT_RHO 2 +#define CT_THETA 3 + +#define CT_FOCAL 1 +#define CT_SX 2 +#define CT_SY 3 + +#define CT_X 4 +#define CT_Y 5 +#define CT_Z 6 +#define CT_THETAN 7 +#define CT_PHIN 8 +#define CT_XN 9 +#define CT_YN 10 +#define CT_ZN 11 + +#define CT_NDATA 12 + +//!@} + +//=------------------------------------------------------------ +//!@subsection data types + +struct camera { /* camera parameters for imaging */ + int inumpixels; + int inumcentralpixels; + int inumgappixels; + int inumbigpixels; + double dpixdiameter_cm; /* diameter of the central and gap pixels in centimeters */ + double dpixsizefactor[iMAXNUMPIX]; /* size of the pixel relative to dpixdiameter_deg */ + double dxc[iMAXNUMPIX]; /* Pixel coordinates in camera coordinate system (x points from pixel 1 to 2). */ + double dyc[iMAXNUMPIX]; /* The numbering of the pixels in these arrays starts at 0! */ + double dxpointcorr_deg; /* correction of the pixel coordinates; to be added to dxc[] to get correct value */ + double dypointcorr_deg; /* correction of the pixel coordinates; to be added to dxc[] to get correct value */ + double di[iMAXNUMPIX]; /* i coordinate in JCs bi-axis hexagonal coordinate system */ + double dj[iMAXNUMPIX]; /* j coordinate in JCs bi-axis hexagonal coordinate system */ + +}; + + +//=------------------------------------------------------------ +//!@subsection Pre-defined file names. + +//!@{ + +#define QE_FILE "qe.dat" + +//!@} + +//=------------------------------------------------------------ +//!@subsection Prototypes of functions. + +//!@{ + +//++ +// prototypes +//-- + +#define ONoff(x) ((x==TRUE) ? "[ ON ]" : "[ off]") + +// Under Linux, the nint function does not exist, so we have to define it. +#define nint(x) ((int)floor((x)+0.5)) + +void present(void); +void usage(void); +void clean(void); +void log(const char *funct, char *fmt, ...); +void error(const char *funct, char *fmt, ...); +int isA( char * s1, const char * flag ); +void read_ct_file(void); +int igen_pixel_coordinates(struct camera *cam); +void read_pixels(struct camera *cam); +int pixels_are_neig(int pix1, int pix2); +int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam); +float dist_r_P(float a, float b, float c, + float l, float m, float n, + float x, float y, float z); + +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + *$Log: not supported by cvs2svn $ + *Revision 1.3 1999/10/22 15:32:56 petry + *tidied-up version, really sent to H.K. and N.M., 22-10-99 + * + *Revision 1.2 1999/10/22 15:01:28 petry + *version sent to H.K. and N.M. on Fri Oct 22 1999 + * + *Revision 1.1.1.1 1999/10/21 16:35:10 petry + *first synthesised version + * + * Revision 1.8 1999/03/15 14:59:06 gonzalez + * camera-1_1 + * + * Revision 1.7 1999/03/02 09:56:11 gonzalez + * *** empty log message *** + * + * Revision 1.6 1999/01/14 17:32:40 gonzalez + * Added to camera the STDIN input option (data_from_input) + * + */ + +//!@} +//=EOF + Index: trunk/MagicSoft/Simulation/Detector/Camera/config.mk.linux-gnu =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/config.mk.linux-gnu (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/config.mk.linux-gnu (revision 308) @@ -0,0 +1,69 @@ +################################################################## +# +# config.mk +# +# @file config.mk +# @title small configuration file for Makefile +# @author J C Gonz\'alez +# @email gonzalez@mppmu.mpg.de +# @date Fri Mar 12 11:51:11 MET 1999 +# +#_______________________________________________________________ +# +# Created: Fri Mar 12 11:51:11 MET 1999 +# Author: Jose Carlos Gonzalez +# Purpose: Makefile for the compilation of the camera program +# Notes: +# +#--------------------------------------------------------------- +# $RCSfile: config.mk.linux-gnu,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 1999-11-05 11:59:32 $ +################################################################## +# @maintitle + +# @code + +# program + +PROGRAM = camera + +# compilers + +CC = g++ +CXX = g++ +F77 = f77 + +DOCUM = ${HOME}/detector/sus/sus + +# includes + +INCLUDE = ../include-GENERAL +INCLUDE_COR = ../include-CORSIKA +INCLUDE_MC = ../include-MC +INCLUDE_TRIGGER = ../include-MTrigger +INCLUDE_EVITA = ../../../include-Classes +INCLUDE_REFL = ../Reflector +INCLUDE_ROOT = ${ROOTSYS}/include +INCLUDE_g++ = /usr/include/g++ + +#OPTIM = -O2 -ieee +OPTIM = -O2 -Wall -fno-rtti -fno-exceptions -fPIC +DEBUG = -g + +# libraries + +RANLIBDIR = ../lib +CERNDIR = /cern + +# system + +SYSTEM = linux + +# uncomment this for quiet compilation + +.SILENT: + +# @endcode +##EOF Index: trunk/MagicSoft/Simulation/Detector/Camera/creadparam.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/creadparam.cxx (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/creadparam.cxx (revision 308) @@ -0,0 +1,856 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= creadparam +//= +//= @file creadparam.cxx +//= @desc Reading of parameters file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: creadparam.cxx,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-11-05 11:59:34 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |creadparam.cxx|. + +/*!@" + + This section describes briefly the source code for the file + |creadparam.cxx|. This file is very closely related to the file + |readparams.cxx| from the |reflector| program. Actually, this later + file was the ancestror of the file you are looking at. + + All the defines it uses are located in the file |creadparam.h|. In + the first one we can see the definitions of the commands available + for the parameters file. We describe these commands in a later + section. + + @"*/ + +//!@subsection Includes and Global variables definition. + +/*!@" + + All the defines are located in the file {\tt creadparam.h}. + + @"*/ + +//!@{ + +#include "creadparam.h" + +//!@} + +//!@subsection Definition of global variables. + +/*!@" + + Here we define the global variables where the values from the + parameters file are stored. + + @"*/ + +//!@{ + +static char Input_filename[PATH_MAX_LENGTH]; //@< input filename +static char Output_filename[PATH_MAX_LENGTH]; //@< output filename +static char Data_filename[PATH_MAX_LENGTH]; //@< data filename +static char DIAG_filename[PATH_MAX_LENGTH]; //@< data filename +static char ROOT_filename[PATH_MAX_LENGTH]; //@< data filename +static char CT_filename[PATH_MAX_LENGTH]; //@< name of the CT def. file +static int simulateNSB = TRUE; //@< Will we simulate NSB? +static int anaPixels = -1; //@< number of pixels for the analysis +static float meanNSB; //@< NSB mean value (per pixel) +static float qThreshold; //@< Threshold value +static float qTailCut; //@< Tail Cut value +static int nIslandsCut; //@< Islands Cut value +static int countIslands = TRUE; //@< Will we count the islands? +static long int Seeds[2]; +static int *Skip; +static int nSkip=0; +static int Data_From_STDIN = FALSE; +static int Read_Phe = FALSE; +static int Read_Phe_All = FALSE; +static int Write_All_Images = FALSE; +static int Write_All_Data = FALSE; +static int Select_Energy = TRUE; +static float Select_Energy_le = 0.0; //@< GeV +static float Select_Energy_ue = 100000.0; //@< GeV +static float Trigger_Radius; +static int Set_Trigger_Radius=FALSE; +static float fCorrection; +static int Apply_Correction=FALSE; +//!@} + +//!@subsection The function |readparam()|. + +//!----------------------------------------------------------- +// @name creadparam +// +// @desc read parameters from the stdin / parameters file +// +// @var *filename Name of the parameters file (NULL->STDIN) +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +readparam(char * filename) +{ + char sign[] = GLUE_postp( PROGRAM, VERSION ); //@< initialize sign + char line[LINE_MAX_LENGTH]; //@< line to get from the stdin + char token[ITEM_MAX_LENGTH]; //@< a single token + int i, j; //@< dummy counters + ifstream ifile; + + // use cin or ifile (reading from STDIN or from parameters file? + if ( filename != NULL ) + ifile.open( filename ); + + // get signature + if ( filename != NULL ) + ifile.getline(line, LINE_MAX_LENGTH); + else + cin.getline(line, LINE_MAX_LENGTH); + line[strlen(SIGNATURE)] = '\0'; + strcpy(line, sign); + if (strcmp(sign, SIGNATURE) != 0) { + cerr << "ERROR: Signature of parameters file is not correct\n"; + cerr << '"' << sign << '"' << '\n'; + cerr << "should be: " << SIGNATURE << '\n'; + exit(1); + } + + // loop till the "end" directive is reached + int is_end = FALSE; + while (! is_end) { + + // get line from file or stdin + if ( filename != NULL ) + ifile.getline(line, LINE_MAX_LENGTH); + else + cin.getline(line, LINE_MAX_LENGTH); + + // skip comments (start with '#') + if (line[0] == '#') + continue; + + // show user comments (start with '>') + if (line[0] == '>') { + cout << line << endl << flush; + continue; + } + + // look for each item at the beginning of the line + for (i=0; i<=end_file; i++) + if (strstr(line, ITEM_NAMES[i]) == line) + break; + + // if it is not a valid line, just ignore it + if (i == end_file+1) { + cerr << "Skipping unknown token in [" << line << "]\n"; + continue; + } + + // case block for each directive + switch ( i ) { + + case input_file: //@< name of the input file + + // get the name of the input_file from the line + sscanf(line, "%s %s", token, Input_filename); + + break; + + case output_file: //@< name of the output file + + // get the name of the output_file from the line + sscanf(line, "%s %s", token, Output_filename); + + break; + + case data_file: //@< name of the data file + + // get the name of the data_file from the line + sscanf(line, "%s %s", token, Data_filename); + + break; + + case diag_file: //@< name of the DIAG file + + // get the name of the data_file from the line + sscanf(line, "%s %s", token, DIAG_filename); + cout << '[' << DIAG_filename << ']' << endl << flush; + + break; + + case root_file: //@< name of the ROOT file + + // get the name of the data_file from the line + sscanf(line, "%s %s", token, ROOT_filename); + cout << '[' << ROOT_filename << ']' << endl << flush; + + break; + + case ct_file: //@< name of the telescope file + + // get the name of the ct_file from the line + sscanf(line, "%s %s", token, CT_filename); + + break; + + case nsb_on: //@< simulate NSB? + + // we will simulate NSB + simulateNSB = TRUE; + + break; + + case nsb_off: //@< simulate NSB? + + // we will NOT simulate NSB + simulateNSB = FALSE; + + break; + + case nsb_mean: //@< value of per pixel + + // get value of (in photons) + sscanf(line, "%s %f", token, &meanNSB); + simulateNSB = TRUE; + + break; + + case ana_pixels: //@< number of pixels for analysis + + // number of pixels for analysis + sscanf(line, "%s %d", token, &anaPixels); + + break; + + case threshold: //@< value of threshold for trigger (q0) + + // get value of threshold (in ph.e.) + sscanf(line, "%s %f", token, &qThreshold); + + break; + + case tail_cut: //@< value of tail_cut (t0) + + // get value of tail_cut (in ph.e.) + sscanf(line, "%s %f", token, &qTailCut); + + break; + + case islands_on: //@< DO count islands + + // DO count islands + countIslands = TRUE; + + break; + + case islands_off: //@< do NOT count islands + + // do NOT count islands + countIslands = FALSE; + + break; + + case islands_cut: //@< value of islands_cut (i0) + + // get value of islands_cut (in ph.e.) + sscanf(line, "%s %d", token, &nIslandsCut); + countIslands = TRUE; + + break; + + case seeds: //@< values of seeds for random numbers + + // get seeds + sscanf(line, "%s %ld %ld", token, &Seeds[0], &Seeds[1]); + + break; + + case skip: //@< skip pathological showers + + // get showers to skip + cin >> nSkip; + Skip = new int[nSkip]; + for (j=0; j> Skip[j]; + cout << Skip[j] << endl << flush; + } + + break; + + case data_from_stdin: //@< to read data from stdin + + // change boolean value + Data_From_STDIN = TRUE; + + break; + + case read_phe: //@< to read PHE files + + // change boolean value + Read_Phe = TRUE; + + break; + + case read_phe_all: //@< to read PHE files from write_all_images + + // change boolean value + Read_Phe_All = TRUE; + + break; + + case write_all_images: //@< to write ALL the images + + // change boolean value + Write_All_Images = TRUE; + + break; + + case write_all_data: //@< to write single pixel data + + // change boolean value + Write_All_Data = TRUE; + + break; + + case select_energy: //@< value of islands_cut (i0) + + // get energy range + sscanf(line, "%s %f %f", token, &Select_Energy_le, &Select_Energy_ue); + Select_Energy = TRUE; + + break; + + case trigger_radius: //@< set radius of trigger area + + // get trigger radius + sscanf(line, "%s %f", token, &Trigger_Radius); + Set_Trigger_Radius = TRUE; + + break; + + case correction: //@< apply a kind of correction to pixel values + + // get trigger radius + sscanf(line, "%s %f", token, &fCorrection); + Apply_Correction = TRUE; + + break; + + case end_file: //@< end of the parameters file + + // show a short message + is_end = TRUE; + + break; + + } // switch ( i ) + + } // while (! is_end) + + // after the loop is finished, return to the main function + return; +} +//!@} + + +//!----------------------------------------------------------- +// @name get_input_filename +// +// @desc get name of the input file +// +// @return Name of the Input file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_input_filename(void) +{ + return (Input_filename); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_output_filename +// +// @desc get name of the output file +// +// @return Name of the Output file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_output_filename(void) +{ + return (Output_filename); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_data_filename +// +// @desc get name of the data file +// +// @return Name of the Data file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_data_filename(void) +{ + return (Data_filename); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_diag_filename +// +// @desc get name of the diagnostic output file +// +// @return Name of the DIAG file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_diag_filename(void) +{ + return (DIAG_filename); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_root_filename +// +// @desc get name of the ROOT file +// +// @return Name of the ROOT file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_root_filename(void) +{ + return (ROOT_filename); +} +//!@} + + + +//!----------------------------------------------------------- +// @name get_ct_filename +// +// @desc get name of CT definition file +// +// @return Name of the CT definition file +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +char * +get_ct_filename(void) +{ + return (CT_filename); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_nsb +// +// @desc are we going to simulate NSB ? +// +// @var *n Mean value for the NSB (ph.e./pixel) +// @return TRUE: we'll simulate NSB; FALSE: we won't +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +int +get_nsb(float *n) +{ + *n = meanNSB; + return ( simulateNSB ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_threshold +// +// @desc get threshold value +// +// @return Value of the threshold q$_0$ +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +float +get_threshold(void) +{ + return( qThreshold ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_tail_cut +// +// @desc get tail cut value +// +// @return Value for the Tail-cut +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +float +get_tail_cut(void) +{ + return( qTailCut ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_islands_cut +// +// @desc are we going to count the islands ? +// +// @var *n Cut on islands number +// @return TRUE: we'll count the islands; FALSE: we won't +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +int +get_islands_cut(int *n) +{ + *n = nIslandsCut; + return ( countIslands ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_ana_pixels +// +// @desc number of pixels for the analysis +// +// @return Number of pixels to use in the image parameters +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +int +get_ana_pixels(void) +{ + return ( anaPixels ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_seeds +// +// @desc are we going to count the islands ? +// +// @var *n Number of the seed +// @return N-th random-number Seed +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +long int +get_seeds(int n) +{ + return ( Seeds[n] ); +} +//!@} + + +//!----------------------------------------------------------- +// @name get_skip_showers +// +// @desc get list of showers to skip +// +// @var *s1 Pointer to a vector of number of showers +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +void +get_skip_showers( int *s ) +{ + int i; + for (i=0; i +#include +#include +#include +#include +#include +#include + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "camera-v.h" + +//!@} + +//!@subsection Macro-definitions, and constants. + +//!@{ + +// now we define the list ITEM_LIST of possible items in +// the parameters file. note that they are functions of +// T(x). we will change T(x) to be the real item or the +// string of this item when needed + +#define ITEM_LIST /* LIST OF ITEMS IN THE PARAMETERS FILE */ \ +T(input_file), /* input file */ \ +T(output_file), /* output file */ \ +T(data_file), /* data file */ \ +T(diag_file), /* diagnostic output file (ROOT format) */ \ +T(root_file), /* ROOT file */ \ +T(ct_file), /* file with the characteristics of the CT */ \ +T(ana_pixels), /* size of the camera for parameters calculation */ \ +T(nsb_on), /* activates NSB simulation */ \ +T(nsb_off), /* de-activates NSB simulation */ \ +T(nsb_mean), /* mean value of NSB contribution per pixel */ \ +T(threshold), /* value of q0 for trigger */ \ +T(tail_cut), /* value of tail cut (t0) */ \ +T(islands_on), /* DO count islands */ \ +T(islands_off), /* do NOT count islands */ \ +T(islands_cut), /* value of islands cut (i0) */ \ +T(seeds), /* seeds for random number generation */ \ +T(data_from_stdin), /* to read data from STDIN */ \ +T(skip), /* skip pathological showers */ \ +T(read_phe_all), /* id., but was processed with write_all_images */ \ +T(read_phe), /* read an already camera processed file */ \ +T(write_all_images),/* write to file .phe ALL images (even w.o. trigger)*/ \ +T(write_all_data), /* write to file .dat ALL image data */ \ +T(select_energy), /* energy range to read: only for .phe files */ \ +T(trigger_radius), /* trigger radius for the camera */ \ +T(correction), /* factor for correction in the pixel values */ \ +T(end_file) /* end of the parameters file */ + +#define T(x) x // define T() as the name as it is + +enum ITEM_TYPE { + ITEM_LIST +}; + +#undef T + +#define T(x) #x // define T() as the string of x + +const char *const ITEM_NAMES[] = { + ITEM_LIST +}; + +#undef T + +#define LINE_MAX_LENGTH 400 +#define ITEM_MAX_LENGTH 40 +#define PATH_MAX_LENGTH 120 + +// mean values of NSB contribution per pixel + +static const float Mean_NSB_MAGIC = 5.0; //@< for MAGIC +static const float Mean_NSB_CT1 = 5.0; //@< for CT1 + +//!@} + +//!@subsection Prototypes of functions. + +//!@{ + +//++ +// prototypes +//-- + +void readparam(char * filename); +char *get_input_filename(void); +char *get_output_filename(void); +char *get_data_filename(void); +char *get_diag_filename(void); +char *get_root_filename(void); +char *get_ct_filename(void); +int get_nsb(float *n); +float get_threshold(void); +float get_tail_cut(void); +int get_islands_cut(int *n); +long int get_seeds(int n); +int get_ana_pixels(void); +void get_skip_showers( int *s ); +int get_nskip_showers( void ); +int get_data_from_stdin(void); +int get_read_phe(void); +int get_read_phe_all(void); +int get_write_all_images(void); +int get_write_all_data(void); +int get_select_energy(float *le, float *ue); +int get_trigger_radius(float *radius); +int get_correction(float *corr); +//!@} + +//!@{ + +#endif // ! _creadparam_ + +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + * $Log: not supported by cvs2svn $ + * Revision 1.2 1999/10/22 15:01:29 petry + * version sent to H.K. and N.M. on Fri Oct 22 1999 + * + * Revision 1.1.1.1 1999/10/21 16:35:11 petry + * first synthesised version + * + * Revision 1.7 1999/03/15 14:59:09 gonzalez + * camera-1_1 + * + * Revision 1.6 1999/03/02 09:56:13 gonzalez + * *** empty log message *** + * + * Revision 1.5 1999/01/14 17:32:43 gonzalez + * Added to camera the STDIN input option (data_from_input) + * + */ + +//!@} +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Camera/moments.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/moments.cxx (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/moments.cxx (revision 308) @@ -0,0 +1,728 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= moments +//= +//= @file moments.cxx +//= @desc Calculation of image parameters +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: moments.cxx,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-11-05 11:59:33 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |moments.cxx|. + +/*!@{ + + This section describes briefly the source code for the file + |moments.cxx|. All the defines it uses are located in the file + |moments.h|. + + @"*/ + +//!@subsection Includes and Global variables definition. + +/*!@" + + All the defines are located in the file |moments.h|. + + @"*/ + +//!@{ + +#include "moments.h" + +//!@} + +//!@subsection Definition of global variables. + +//!@{ + +static int npix; //@< number of pixels +static float *q; //@< charges in the pixels +static float xm, ym; //@< centroid (used in moments and lenwid) + +//@: structure with information about the image +static Moments_Info m; + +//@: structure with information about islands +static Islands_Info is; + +//@: structure with information about lenwid +static LenWid_Info lw; + +//!@} + +//!@subsection The function |moments()|. + +//!----------------------------------------------------------- +// @name moments +// +// @desc calculate moments on the image +// +// @var n Number of pixels +// @var *image Vector of ph.e.s in pixels +// @var **pix Array with information about the pixels +// @var plateScale Plate scale for the CT in use +// @var flag 1: initialize; other: normal +// +// @return Pointer to structure Moments_Info +// +// @date Mon Sep 14 15:22:44 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +Moments_Info * +moments( int n, float *image, float **pix, + float plateScale, int flag ) +{ + register int i, k; + + float x, y; + float x2m, xym, y2m; + float x3m, x2ym, xy2m, y3m; + float zz, zd, zu, zv; + float ax, ay, unitx, unity, sigmaax; + float sx2, sxy, sy2; + float sx3, sx2y, sxy2, sy3; + + if ( flag == 1 ) { + q = new float[n]; + is.fi = new float[n]; + is.vislands = new float[n]; + is.islands = new int[n]; + is.isl = new int[n]; + for (i=1; i 0.0 ) { + x = pix[i][0]; + y = pix[i][1]; + xm += x * q[i]; + ym += y * q[i]; + x2m += x * x * q[i]; + xym += x * y * q[i]; + y2m += y * y * q[i]; + x3m += x * x * x * q[i]; + x2ym += x * x * y * q[i]; + xy2m += x * y * y * q[i]; + y3m += y * y * y * q[i]; + m.charge += q[i]; + } + + xm *= plateScale; + ym *= plateScale; + x2m *= plateScale * plateScale; + xym *= plateScale * plateScale; + y2m *= plateScale * plateScale; + x3m *= plateScale * plateScale * plateScale; + x2ym *= plateScale * plateScale * plateScale; + xy2m *= plateScale * plateScale * plateScale; + y3m *= plateScale * plateScale * plateScale; + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // extremes and charges + //-------------------------------------------------- + + for (i=0; i<10; ++i) + m.maxs[i] = 0.0; + + for (i=0; i m.maxs[0] ) { + for (k=9; k>0; --k) + m.maxs[k] = m.maxs[k-1]; + for (k=9; k>0; --k) + m.nmaxs[k] = m.nmaxs[k-1]; + m.maxs[0] = q[i]; + m.nmaxs[0] = i; + } + } + + // calculates weighted position of the maximum (6 pixels) + + m.xmax = m.ymax = m.smax = 0.0; + + for (i=0; i<6; ++i) { + m.xmax += pix[m.nmaxs[i]][0] * q[m.nmaxs[i]]; + m.ymax += pix[m.nmaxs[i]][1] * q[m.nmaxs[i]]; + m.smax += q[m.nmaxs[i]]; + } + + if (m.smax==0.) m.smax=1.; + if (m.charge==0.) m.charge=1.; + + m.xmax = m.xmax * plateScale / m.smax; + m.ymax = m.ymax * plateScale / m.smax; + + // calculate concentrations with 2,3,4...10 pixels + + m.conc[0] = q[ m.nmaxs[0] ] + q[ m.nmaxs[1] ]; + + for (i=2; i<10; ++i) + m.conc[i-1] = m.conc[i-2] + q[ m.nmaxs[i] ]; + + for (i=0; i<9; ++i) + m.conc[i] /= m.charge; + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // 1st moments + //-------------------------------------------------- + + xm /= m.charge; + ym /= m.charge; + + m.m1x = xm; + m.m1y = ym; + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // 2nd moments + //-------------------------------------------------- + + x2m /= m.charge; + xym /= m.charge; + y2m /= m.charge; + + // around the origin + m.m2xx = x2m; + m.m2xy = xym; + m.m2yy = y2m; + + // around the mean + sx2 = x2m - SQR(xm); + sxy = xym - xm * ym; + sy2 = y2m - SQR(ym); + + m.m2cxx = sx2; + m.m2cxy = sxy; + m.m2cyy = sy2; + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // 3rd moments + //-------------------------------------------------- + + x3m /= m.charge; + x2ym /= m.charge; + xy2m /= m.charge; + y3m /= m.charge; + + // around the origin + m.m3xxx = x3m; + m.m3xxy = x2ym; + m.m3xyy = xy2m; + m.m3yyy = y3m; + + // around the mean + sx3 = x3m - 3 * x2m * xm + 2 * xm * xm * xm; + sx2y = x2ym - 2 * xym * xm + 2 * xm * xm * ym - x2m * ym; + sxy2 = xy2m - 2 * xym * ym + 2 * xm * ym * ym - y2m * xm; + sy3 = y3m - 3 * y2m * ym + 2 * ym * ym * ym; + + m.m3cxxx = x3m; + m.m3cxxy = x2ym; + m.m3cxyy = xy2m; + m.m3cyyy = y3m; + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // hillas' parameters + //-------------------------------------------------- + + zd = sy2 - sx2; + zz = sqrt( SQR(zd) + 4.*SQR( sxy ));; + if ( (zz < 1.e-6) || (sxy == 0.) ) { + m.dist = -1.; + return &m; + } + zu = 1.0 + zd / zz; + zv = 2.0 - zu; + + /* + a = (zd + zz) / (2 * sxy); + b = ym - a * xm; + + m.length = sqrt( fabs( sx2 + 2 * a * sxy + a * a * sy2 ) / (1+a*a) ); + m.width = sqrt( fabs( sx2 - 2 * a * sxy + a * a * sy2 ) / (1+a*a) ); + m.dist = sqrt( SQR(xm) + SQR(ym) ); + m.xdist = sqrt( SQR( m.xmax ) + SQR( m.ymax ) ); + m.azw = sqrt( fabs( SQR(xm)* y2m - 2.* xm * ym * xym + x2m*SQR(ym) ) ); + m.miss = fabs( b / sqrt(1+a*a) ); + m.alpha = DEG( asin( m.miss / m.dist ) ); + */ + + m.length = sqrt( fabs(sx2 + sy2 + zz) / 2. ); + m.width = sqrt( fabs(sx2 + sy2 - zz) / 2. ); + m.dist = sqrt( SQR(xm) + SQR(ym) ); + m.xdist = sqrt( SQR(m.xmax) + SQR(m.ymax) ); + m.azw = sqrt( fabs( SQR(xm)*y2m - 2.*xm*ym*xym + x2m*SQR(ym) ) ) / m.dist; + m.miss = sqrt( fabs( (SQR(xm)*zu + SQR(ym)*zv)/2. - (2.*sxy*xm*ym/zz) ) ); + m.alpha = DEG( asin( m.miss/m.dist ) ); + + + /* + length = sqrt( fabs(sx2 + sy2 + zz) /2. ); + width = sqrt( fabs(sx2 + sy2 - zz) / 2. ); + dist = sqrt( SQR(xm) + SQR(ym) ); + xdist = sqrt( SQR(m.xmax) + SQR(m.ymax) ); + azw = sqrt( fabs( SQR(xm)*y2m - 2.*xm*ym*xym + x2m*SQR(ym) ) ); + miss = sqrt( fabs( (SQR(xm)*zu + SQR(ym)*zv)/2. - (2.*sxy*xm*ym/zz) ) ); + alpha = DEG( asin(miss/dist) ); + */ + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // asymetry + //-------------------------------------------------- + + unitx = sqrt(0.5*zv); + unity = SGN( sxy )*sqrt(0.5*zu); + + if ( m.xdist > 0.0 ) { + + m.phi = acos((unitx*m.xmax + unity*m.ymax )/m.xdist); + + sigmaax = sx3*CUB(cos(m.phi)) + 3.0*sx2y*SQR(cos(m.phi))*sin(m.phi) + + 3.0*sxy2*cos(m.phi)*SQR(sin(m.phi)) + sy3*CUB(sin(m.phi)); + sigmaax = pow(fabs(sigmaax),0.3333333)*SGN(sigmaax); + + ax = sigmaax*unitx; + ay = sigmaax*unity; + m.asymx = (ax*m.xmax + ay*m.ymax)/(m.xdist*m.length*cos(m.phi)); + m.asymy = 0.0; + + } else { + + m.phi=-1000.0; + m.asymx = -1000.0; + m.asymy = -1000.0; + + } + + /* + cout + << "length "<< length << endl + << "width "<< width << endl + << "dist "<< dist << endl + << "xdist "<< xdist << endl + << "azw "<< azw << endl + << "miss "<< miss << endl + << "alpha "<< alpha << endl << flush; + */ + + return &m; +} +//!@} + +// @T \newpage + +//!@subsection The function |islands()|. + +//!----------------------------------------------------------- +// @name islands +// +// @desc implementation of the "islands" algorithm +// +// @var n Number of pixels +// @var f Vector with the image (ph.e.s in pixels) +// @var **pixneig Array with indices of neighbour pixels +// @var *npixneig Vector with number of neighbours per pixel +// @var cleanning TRUE: remove spurious islands +// @var ipixcut Islands number cut +// +// @return Pointer to structure Islands_Info +// +// @date Mon Sep 14 15:22:44 MET DST 1998 +//------------------------------------------------------------ +// @function + +//!@{ +Islands_Info * +islands( int n, float *f, int **pixneig, int *npixneig, + int cleanning, int ipixcut) +{ + register int i; + int j, k; + int haschanged; + + is.numisl = 0; + + memcpy( is.fi, f, sizeof(float) * n ); + + // be aware: here we use the side effect of ++ + // there are two possibilities of using the operator ++: + // 1) a = ++i => is.first increments i, then evaluates expresion + // 2) a = i++ => is.first evaluates expresion, then increments i + // we INTENTIONALLY use the second form + + // algorithm to isolate/detect is.islands + j=1; + for (i=0; i0.0 ) { + is.isl[i] = j; + ++j; + } else { + is.isl[i] = 0; + } + + haschanged = TRUE; + while ( haschanged ) { + haschanged = FALSE; + for (i=0; i 0 ) + for (j=0; (j-1); ++j) + if ( (k=is.isl[pixneig[i][j]]) > 0 ) + if ( is.isl[i] > is.isl[k] ) { + is.isl[i] = is.isl[k]; + haschanged=TRUE; + } + } + + // count is.islands + + for (i=0;i0) { + is.islands[is.isl[i]]++; + is.vislands[is.isl[i]] += is.fi[i]; + } + + for (i=0,j=0,is.numisl=0; i0) { + j++; + //cout << '#' << j << ':' << is.islands[i] << " q=" << is.vislands[i] + // << endl; + + if (is.islands[i] > ipixcut) + is.numisl++; + } + + } + + cout << j << '[' << is.numisl << "] is.islands\n" << flush; + + if ( cleanning ) { + + // cleanning image: pixcut = 3 (any is.island with <= 3 pixels is removed + for (i=0;i forget that pixel + + if ( (fabs(dist_to_axis) > max_distance) && + (SGN(dist_to_axis) != sign_of_semiplane) ) + continue; + + // (B) + + // if distance to the axis if larger than pixel diameter, + // AND + // the semiplane is the GOOD one -> add this pixel + + if ( (fabs(dist_to_axis) > max_distance) && + (SGN(dist_to_axis) == sign_of_semiplane) ) { + + // here the sum + + weight = image[i]; + j = k+sign_of_semiplane; + wsum[j] += weight * dist_to_axis * dist_to_axis; + sum[j] += weight; + + continue; + } + + // (C) + // if we reach this point, that means that the center + // of our pixel is too close to the axis, and we have + // to feed it into the routine to check if the pixel + // crosses the axis + + // ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE + // simplified algorithm + // assume the pixels are circular, and takes + // the fraction of the surface lying on the semiplane + // ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE ** NOTE + + // alpha + alpha = 2*asin(sqrt(2*(radius-dist_to_axis)*radius - + radius2) / radius); + + // here the sum + // the fraction is the fraction of the area inside the semiplane + weight = image[i] * ( (alpha * radius2 / 2.0) / (M_PI * radius2)); + j = k+sign_of_semiplane; + wsum[j] += weight * dist_to_axis * dist_to_axis; + sum[j] += weight; + + } // foreach pixel pixels + + } // foreach semiplane + + } // foreach axis + + lw.length1 = (sum[0] > 0.) ? sqrt(wsum[0] / sum[0]) : -1; + lw.width1 = (sum[1] > 0.) ? sqrt(wsum[1] / sum[1]) : -1; + lw.length2 = (sum[2] > 0.) ? sqrt(wsum[2] / sum[2]) : -1; + lw.width2 = (sum[3] > 0.) ? sqrt(wsum[3] / sum[3]) : -1; + + lw.length1 *= plateScale; + lw.width1 *= plateScale; + lw.length2 *= plateScale; + lw.width2 *= plateScale; + + return &lw; +} +//!@} + + +//!@subsection Auxiliary functions. + +//!----------------------------------------------------------- +// @name crosspt +// +// @desc calculate cross point of segments AB and CD +// +// @var ax Coor. X of point A +// @var ay Coor. Y of point A +// @var bx Coor. X of point A +// @var by Coor. Y of point A +// @var cx Coor. X of point A +// @var cy Coor. Y of point A +// @var dx Coor. X of point A +// @var dy Coor. Y of point A +// @var *pcrossx Coor. X of cross point +// @var *pcrossy Coor. Y of cross point +// +// @date Mon Mar 8 13:35:54 MET 1999 +//------------------------------------------------------------ +// @function + +//!@{ +void +crosspt( float ax, float ay, + float bx, float by, + float cx, float cy, + float dx, float dy, + float * pcrossx, float * pcrossy) +{ + float w, r; + + // the points A and B, and C and D define two segments (AB and CD) + // the coordinates of these points are + // A(ax,ay), B(bx,by), C(cx,cy), D(dx,dy) + + w=(bx-ax)*(dy-cy)-(by-ay)*(dx-cx); + r=(ay-cy)*(dx-cx)-(ax-cx)*(dy-cy); + + *pcrossx = ax + r*(bx-ax)/w; + *pcrossy = ay + r*(by-ay)/w; + +} +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ +// +// $Log: not supported by cvs2svn $ +// Revision 1.2 1999/10/22 15:01:29 petry +// version sent to H.K. and N.M. on Fri Oct 22 1999 +// +// Revision 1.1.1.1 1999/10/21 16:35:10 petry +// first synthesised version +// +// Revision 1.1 1999/03/08 10:04:06 gonzalez +// *** empty log message *** +// +// Revision 1.4 1999/03/02 09:56:14 gonzalez +// *** empty log message *** +// +// Revision 1.5 1999/03/15 14:59:10 gonzalez +// camera-1_1 +// +//!@} + +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Camera/moments.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/moments.h (revision 308) +++ trunk/MagicSoft/Simulation/Detector/Camera/moments.h (revision 308) @@ -0,0 +1,158 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= moments +//= +//= @file moments.h +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: moments.h,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-11-05 11:59:33 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |moments.h|. + +/*!@" + + In this section you can find the source code for the file + |moments.h|. This file is mainly needed by |moments.cxx|. + + @"*/ + +//!@{ + +#ifndef _moments_ +#define _moments_ + +#ifndef _this_ +#define _this_ moments +#endif + +//!@} + +//!@subsection Include files. + +//!@{ + +#include +#include +#include +#include +#include +#include +#include + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "camera-v.h" + +//!@} + +//!@subsection Macro-definitions, and constants. + +//!@{ + +typedef struct { + + // moments + float m1x, m1y; // first moments (mean) + float m2xx, m2xy, m2yy; // second moments (around origin) + float m2cxx, m2cxy, m2cyy; // second moments (around mean) + float m3xxx, m3xxy, m3xyy, m3yyy; // third moments (around origin) + float m3cxxx, m3cxxy, m3cxyy, m3cyyy; // third moments (around mean) + + // charges + float charge; // total charge in the image + float xmax, ymax; // position of the maximum + float smax; // charge of the block of maximum + float maxs[10]; // charges of the first 10 max. + int nmaxs[10]; // number of pixels of 10 max. + + // parameters of the image + float length, width, dist, xdist, azw, miss, alpha, conc[9]; + float phi, asymx, asymy; + +} Moments_Info; + + +typedef struct { + float *fi; + int *isl, *islands; + float *vislands; + int numisl; +} Islands_Info; + + +typedef struct { + float length1, length2; + float width1, width2; +} LenWid_Info; + +//!@} + +//!@subsection Prototypes of functions. + +//!@{ + +//++ +// prototypes +//-- + +Moments_Info * moments( int n, float *image, float **pix, + float plateScale, int flag); +Islands_Info * islands( int n, float *f, int **pixneig, int *npixneig, + int cleanning, int ipixcut ); +LenWid_Info * lenwid( int n, float *image, float **pix, + float plateScale, float max_distance); + +void crosspt( float ax, float ay, + float bx, float by, + float cx, float cy, + float dx, float dy, + float * pcrossx, float * pcrossy); + +//!@} + +//!@{ + +#endif // ! _moments_ + +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + *$Log: not supported by cvs2svn $ + *Revision 1.1.1.1 1999/10/21 16:35:10 petry + *first synthesised version + * + * Revision 1.4 1999/03/15 14:59:10 gonzalez + * camera-1_1 + * + * Revision 1.3 1999/03/02 09:56:15 gonzalez + * *** empty log message *** + * + */ + +//!@} +//=EOF