| 1 |
|
|---|
| 2 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 3 | // //
|
|---|
| 4 | // MDCA.cc //
|
|---|
| 5 | // Container to store Hillas parameters and DCA stuff //
|
|---|
| 6 | // //
|
|---|
| 7 | // Author(s): S.C. Commichau, L.S. Stark, 7/2003 //
|
|---|
| 8 | // //
|
|---|
| 9 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 10 |
|
|---|
| 11 | #include "MDCA.h"
|
|---|
| 12 |
|
|---|
| 13 | using namespace std;
|
|---|
| 14 |
|
|---|
| 15 | ClassImp(MDCA);
|
|---|
| 16 |
|
|---|
| 17 | // Default constructor
|
|---|
| 18 | MDCA::MDCA(const char *name, const char *title)
|
|---|
| 19 | {
|
|---|
| 20 | fName = name ? name : "MDCA";
|
|---|
| 21 | fTitle = title ? title : "Storage container for Hillas parameters and DCA of one event";
|
|---|
| 22 |
|
|---|
| 23 | Reset();
|
|---|
| 24 |
|
|---|
| 25 | fEllipse = new TEllipse;
|
|---|
| 26 | fRefCircle = new TEllipse;
|
|---|
| 27 |
|
|---|
| 28 | fLineL = new TLine;
|
|---|
| 29 | fLineW = new TLine;
|
|---|
| 30 | fLineX = new TLine;
|
|---|
| 31 | fLineY = new TLine;
|
|---|
| 32 | fLineDCA = new TLine;
|
|---|
| 33 | fLineMean = new TLine;
|
|---|
| 34 |
|
|---|
| 35 | }
|
|---|
| 36 |
|
|---|
| 37 | // Destructor: Deletes ellipse and lines if they do exist
|
|---|
| 38 |
|
|---|
| 39 | MDCA::~MDCA()
|
|---|
| 40 | {
|
|---|
| 41 | Clear();
|
|---|
| 42 | }
|
|---|
| 43 |
|
|---|
| 44 | // Initialize parameters with default values
|
|---|
| 45 |
|
|---|
| 46 | void MDCA::Reset()
|
|---|
| 47 | {
|
|---|
| 48 |
|
|---|
| 49 | fLength = -1;
|
|---|
| 50 | fWidth = -1;
|
|---|
| 51 | fDelta0 = 0;
|
|---|
| 52 | fMeanX = 0;
|
|---|
| 53 | fMeanY = 0;
|
|---|
| 54 | fDelta1 = 0;
|
|---|
| 55 | fDCA = -1000;
|
|---|
| 56 | fX1W = 0;
|
|---|
| 57 | fY1W = 0;
|
|---|
| 58 | fX2W = 0;
|
|---|
| 59 | fY2W = 0;
|
|---|
| 60 | fX1L = 0;
|
|---|
| 61 | fY1L = 0;
|
|---|
| 62 | fX2L = 0;
|
|---|
| 63 | fY2L = 0;
|
|---|
| 64 | fXDCA = 0;
|
|---|
| 65 | fYDCA = 0;
|
|---|
| 66 |
|
|---|
| 67 | }
|
|---|
| 68 |
|
|---|
| 69 |
|
|---|
| 70 | // Print parameters to *fLog
|
|---|
| 71 |
|
|---|
| 72 | void MDCA::Print(Option_t *) const
|
|---|
| 73 | {
|
|---|
| 74 | Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;
|
|---|
| 75 |
|
|---|
| 76 | if (atg<0)
|
|---|
| 77 | atg += 180;
|
|---|
| 78 |
|
|---|
| 79 | *fLog << all;
|
|---|
| 80 | *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
|
|---|
| 81 | *fLog << " - Length [mm] = " << fLength << endl;
|
|---|
| 82 | *fLog << " - Width [mm] = " << fWidth << endl;
|
|---|
| 83 | *fLog << " - Delta0 [deg] = " << fDelta0*kRad2Deg << endl;
|
|---|
| 84 | *fLog << " - Meanx [mm] = " << fMeanX << endl;
|
|---|
| 85 | *fLog << " - Meany [mm] = " << fMeanY << endl;
|
|---|
| 86 | *fLog << " - atg(y/x) [deg] = " << atg << endl;
|
|---|
| 87 | *fLog << " - DCA [mm] = " << fDCA << endl;
|
|---|
| 88 | *fLog << " - Delta1 [deg] = " << fDelta1*kRad2Deg << endl;
|
|---|
| 89 |
|
|---|
| 90 |
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | void MDCA::Paint(Option_t *opt)
|
|---|
| 94 | {
|
|---|
| 95 | Clear();
|
|---|
| 96 |
|
|---|
| 97 | if (fLength<=0 || fWidth<=0) //fLength<0 || fWidth<0 doesn't look nice...
|
|---|
| 98 | return; //We get a circle!
|
|---|
| 99 |
|
|---|
| 100 | // Length line
|
|---|
| 101 | fLineL = new TLine(fX1L, fY1L, fX2L, fY2L);
|
|---|
| 102 | fLineL->SetLineWidth(2);
|
|---|
| 103 | fLineL->SetLineColor(2);
|
|---|
| 104 | fLineL->Draw();
|
|---|
| 105 |
|
|---|
| 106 | // Width line
|
|---|
| 107 | fLineW = new TLine(fX1W, fY1W, fX2W, fY2W);
|
|---|
| 108 | fLineW->SetLineWidth(2);
|
|---|
| 109 | fLineW->SetLineColor(2);
|
|---|
| 110 | fLineW->Draw();
|
|---|
| 111 |
|
|---|
| 112 | // Coordinate system
|
|---|
| 113 | fLineX = new TLine(-600,fYRef,600,fYRef);
|
|---|
| 114 | fLineY = new TLine(fXRef,-600,fXRef,600);
|
|---|
| 115 | fLineX->SetLineWidth(1);
|
|---|
| 116 | fLineX->SetLineColor(108);
|
|---|
| 117 | fLineY->SetLineWidth(1);
|
|---|
| 118 | fLineY->SetLineColor(108);
|
|---|
| 119 | fLineX->Draw();
|
|---|
| 120 | fLineY->Draw();
|
|---|
| 121 |
|
|---|
| 122 | // DCA line
|
|---|
| 123 | fLineDCA = new TLine(fXRef,fYRef,fXDCA+fXRef,fYDCA+fYRef);
|
|---|
| 124 | fLineDCA->SetLineWidth(2);
|
|---|
| 125 | fLineDCA->SetLineColor(2);
|
|---|
| 126 | fLineDCA->Draw();
|
|---|
| 127 |
|
|---|
| 128 | // COG line
|
|---|
| 129 | fLineMean = new TLine(fXRef,fYRef,fMeanX,fMeanY);
|
|---|
| 130 | fLineMean->SetLineWidth(2);
|
|---|
| 131 | fLineMean->SetLineColor(2);
|
|---|
| 132 | fLineMean->Draw();
|
|---|
| 133 |
|
|---|
| 134 | // Reference point marker
|
|---|
| 135 | fRefCircle = new TEllipse(fXRef, fYRef, 10, 10, 0, 360, 0);
|
|---|
| 136 | fRefCircle->SetLineColor(108);
|
|---|
| 137 | fRefCircle->SetFillColor(108);
|
|---|
| 138 | fRefCircle->Draw();
|
|---|
| 139 |
|
|---|
| 140 | // Hillas ellipse
|
|---|
| 141 | fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta0*kRad2Deg+180);
|
|---|
| 142 | fEllipse->SetLineWidth(2);
|
|---|
| 143 | fEllipse->SetLineColor(2);
|
|---|
| 144 | fEllipse->Draw();
|
|---|
| 145 |
|
|---|
| 146 | }
|
|---|
| 147 |
|
|---|
| 148 |
|
|---|
| 149 | // If an ellipse and lines exist they will be deleted
|
|---|
| 150 |
|
|---|
| 151 | void MDCA::Clear(Option_t *)
|
|---|
| 152 | {
|
|---|
| 153 | if (!fEllipse && !fRefCircle && !fLineL && !fLineW && !fLineX && !fLineY && !fLineDCA && !fLineMean)
|
|---|
| 154 | return;
|
|---|
| 155 |
|
|---|
| 156 | delete fEllipse;
|
|---|
| 157 | delete fRefCircle;
|
|---|
| 158 | delete fLineL;
|
|---|
| 159 | delete fLineW;
|
|---|
| 160 | delete fLineX;
|
|---|
| 161 | delete fLineY;
|
|---|
| 162 | delete fLineDCA;
|
|---|
| 163 | delete fLineMean;
|
|---|
| 164 |
|
|---|
| 165 | fLineL = NULL;
|
|---|
| 166 | fLineX = NULL;
|
|---|
| 167 | fLineY = NULL;
|
|---|
| 168 | fLineW = NULL;
|
|---|
| 169 | fLineDCA = NULL;
|
|---|
| 170 | fLineMean = NULL;
|
|---|
| 171 | fEllipse = NULL;
|
|---|
| 172 | fRefCircle = NULL;
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 |
|
|---|
| 176 | Int_t MDCA::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil)
|
|---|
| 177 | {
|
|---|
| 178 | // Get basic Hillas parameters from MHillas
|
|---|
| 179 | fDelta0 = hil.GetDelta();
|
|---|
| 180 | fMeanX = hil.GetMeanX();
|
|---|
| 181 | fMeanY = hil.GetMeanY();
|
|---|
| 182 | fLength = hil.GetLength();
|
|---|
| 183 | fWidth = hil.GetWidth();
|
|---|
| 184 |
|
|---|
| 185 | // The Length Line - rotation and shift
|
|---|
| 186 | fX1L = - (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm]
|
|---|
| 187 | fY1L = - (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm]
|
|---|
| 188 | fX2L = (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm]
|
|---|
| 189 | fY2L = (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm]
|
|---|
| 190 |
|
|---|
| 191 | // The Width Line - rotation and shift
|
|---|
| 192 | fX1W = (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm]
|
|---|
| 193 | fY1W = - (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm]
|
|---|
| 194 | fX2W = - (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm]
|
|---|
| 195 | fY2W = (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm]
|
|---|
| 196 |
|
|---|
| 197 | // Vector of orientation of the shower axis
|
|---|
| 198 | fr1 = fX2L-fX1L;
|
|---|
| 199 | fr2 = fY2L-fY1L;
|
|---|
| 200 |
|
|---|
| 201 | // Determine parameters to calculate coordinates of the DCA vector
|
|---|
| 202 | flambda = (fr1*(fXRef-fMeanX) + fr2*(fYRef-fMeanY))/(fr1*fr1 + fr2*fr2);
|
|---|
| 203 | fmu = (fMeanY-fYRef)/fr1 + flambda*fr2/fr1;
|
|---|
| 204 |
|
|---|
| 205 | // Components of the DCA vector
|
|---|
| 206 | fXDCA = -fmu*fr2;
|
|---|
| 207 | fYDCA = fmu*fr1;
|
|---|
| 208 |
|
|---|
| 209 | // Components of vector going from intersection point of the DCA vector
|
|---|
| 210 | // with the shower axis to the COG
|
|---|
| 211 | fd1 = fMeanX + fmu*fr2 - fXRef;
|
|---|
| 212 | fd2 = fMeanY - fmu*fr1 - fYRef;
|
|---|
| 213 |
|
|---|
| 214 | // Calculate DCA value
|
|---|
| 215 | fDCA = sqrt(fXDCA*fXDCA + fYDCA*fYDCA);
|
|---|
| 216 |
|
|---|
| 217 | // Calculate angle of the shower axis with respect to the x-axis
|
|---|
| 218 | fDelta1 = acos(fd1/sqrt(fd1*fd1 + fd2*fd2)); //botch
|
|---|
| 219 |
|
|---|
| 220 | // Calculate angle of the shower axis with respect to the y-axis
|
|---|
| 221 | //fDelta1 = acos(fd2/sqrt(fd1*fd1 + fd2*fd2));
|
|---|
| 222 |
|
|---|
| 223 | // Determine the correct sign of the DCA (cross product of DCA vector and
|
|---|
| 224 | // vector going from the intersection point of the DCA vector with the shower axis
|
|---|
| 225 | // to the COG)
|
|---|
| 226 | if((fmu*(-fr2*(fMeanY-fYRef)-fr1*(fMeanX-fXRef)))<0)
|
|---|
| 227 | fDCA = -fDCA;
|
|---|
| 228 |
|
|---|
| 229 | gRandom->Rannor(gx,gy);
|
|---|
| 230 | gx = fabs(gx);
|
|---|
| 231 |
|
|---|
| 232 | // This is nice but does not remove the systematics in the profile plot...
|
|---|
| 233 | //if(((1-0.6*gx)*(180-kRad2Deg*fDelta1)>120) || ((1-0.6*gx)*kRad2Deg*fDelta1>120))
|
|---|
| 234 | // fDCA = -1;
|
|---|
| 235 |
|
|---|
| 236 | // Enlarge the interval of Delta correctly...
|
|---|
| 237 | if((fMeanY-fYRef-fmu*fr1)<0)
|
|---|
| 238 | fDelta1 = TwoPi-fDelta1;
|
|---|
| 239 |
|
|---|
| 240 | // Enlarge the interval of Delta correctly... (Delta with respect to the y-axis)
|
|---|
| 241 | // if(-(fMeanX-fXRef+fmu*fr2)<0)
|
|---|
| 242 | // fDelta1 = TwoPi-fDelta1;
|
|---|
| 243 |
|
|---|
| 244 | // This has to be improved...
|
|---|
| 245 | if(fr1 == 0 || fr2 == 0 || (fr1*fr1+fr2*fr2) == 0 || sqrt(fd1*fd1+fd2*fd2) == 0)
|
|---|
| 246 | fDCA = -fDCA;
|
|---|
| 247 |
|
|---|
| 248 | SetReadyToSave();
|
|---|
| 249 |
|
|---|
| 250 | return 0;
|
|---|
| 251 | }
|
|---|
| 252 |
|
|---|
| 253 | void MDCA::SetRefPoint(const Float_t fXRef0, const Float_t fYRef0)
|
|---|
| 254 | {
|
|---|
| 255 | fXRef = fXRef0;
|
|---|
| 256 | fYRef = fYRef0;
|
|---|
| 257 | }
|
|---|
| 258 |
|
|---|
| 259 |
|
|---|
| 260 |
|
|---|
| 261 |
|
|---|
