AraCorrelator/AraEventCorrelator.cxx
00001 #include <iostream> 00002 #include "AraEventCorrelator.h" 00003 #include "AraGeomTool.h" 00004 #include "AraAntennaInfo.h" 00005 #include "UsefulAraEvent.h" 00006 #include "FFTtools.h" 00007 #include "TH2D.h" 00008 #include "TMath.h" 00009 00010 ClassImp(AraEventCorrelator); 00011 00012 AraEventCorrelator * AraEventCorrelator::fgInstance=0; 00013 00014 00015 AraEventCorrelator::AraEventCorrelator(Int_t numAnts) 00016 { 00017 //Default constructor 00018 fNumAnts=numAnts; 00019 fNumPairs=0; 00020 for(int first=0;first<(fNumAnts-1);first++) { 00021 for(int second=first+1;second<fNumAnts;second++) { 00022 fFirstAnt[fNumPairs]=first; 00023 fSecondAnt[fNumPairs]=second; 00024 fNumPairs++; 00025 } 00026 } 00027 00028 00029 // fillDeltaTArrays(); 00030 fillAntennaPositions(); 00031 setupDeltaTInfinity(); 00032 setupDeltaT40m(); 00033 } 00034 00035 AraEventCorrelator::~AraEventCorrelator() 00036 { 00037 //Default destructor 00038 } 00039 00040 00041 //______________________________________________________________________________ 00042 AraEventCorrelator* AraEventCorrelator::Instance(Int_t numAnts) 00043 { 00044 //static function 00045 if(fgInstance) 00046 return fgInstance; 00047 00048 fgInstance = new AraEventCorrelator(numAnts); 00049 return fgInstance; 00050 } 00051 00052 void AraEventCorrelator::fillAntennaPositions() 00053 {AraGeomTool *araGeom=AraGeomTool::Instance(); 00054 for(int ant=0;ant<ANTS_PER_STATION;ant++) { 00055 int antPolNum=araGeom->fAntInfo[ant].antPolNum; 00056 std::cerr << ant << "\t" << antPolNum << "\t" << araGeom->fAntInfo[ant].polType << "\n"; 00057 if(araGeom->fAntInfo[ant].polType==AraAntPol::kVertical) { 00058 if(antPolNum<7) { 00059 fRfChanVPol[antPolNum]=ant; 00060 fVPolPos[antPolNum][0]=araGeom->fAntInfo[ant].antLocation[0]; 00061 fVPolPos[antPolNum][1]=araGeom->fAntInfo[ant].antLocation[1]; 00062 fVPolPos[antPolNum][2]=araGeom->fAntInfo[ant].antLocation[2]; 00063 fVPolRho[antPolNum]=TMath::Sqrt(fVPolPos[antPolNum][0]*fVPolPos[antPolNum][0]+ 00064 fVPolPos[antPolNum][1]*fVPolPos[antPolNum][1]); 00065 fVPolPhi[antPolNum]=TMath::ATan2(fVPolPos[antPolNum][1],fVPolPos[antPolNum][0]); 00066 } 00067 } 00068 if(araGeom->fAntInfo[ant].polType==AraAntPol::kHorizontal) { 00069 if(antPolNum<7) { 00070 fRfChanHPol[antPolNum]=ant; 00071 fHPolPos[antPolNum][0]=araGeom->fAntInfo[ant].antLocation[0]; 00072 fHPolPos[antPolNum][1]=araGeom->fAntInfo[ant].antLocation[1]; 00073 fHPolPos[antPolNum][2]=araGeom->fAntInfo[ant].antLocation[2]; 00074 fHPolRho[antPolNum]=TMath::Sqrt(fHPolPos[antPolNum][0]*fHPolPos[antPolNum][0]+ 00075 fHPolPos[antPolNum][1]*fHPolPos[antPolNum][1]); 00076 fHPolPhi[antPolNum]=TMath::ATan2(fHPolPos[antPolNum][1],fHPolPos[antPolNum][0]); 00077 } 00078 } 00079 } 00080 for(int i=0;i<7;i++) { 00081 std::cout << "V\t" << i << "\t" << fVPolPos[i][0] << "\t" << fVPolPos[i][1] << "\t" << fVPolPos[i][2] << "\n"; 00082 } 00083 for(int i=0;i<7;i++) { 00084 std::cout << "H\t" << i << "\t" << fHPolPos[i][0] << "\t" << fHPolPos[i][1] << "\t" << fHPolPos[i][2] << "\n"; 00085 } 00086 00087 //Now fill the arrays with angles 00088 Double_t deltaPhi=360./NUM_PHI_BINS; 00089 Double_t deltaTheta=180./NUM_THETA_BINS; 00090 00091 for(int i=0;i<NUM_PHI_BINS;i++) { 00092 fPhiWaveDeg[i]=-180+0.5*deltaPhi+deltaPhi*i; 00093 fPhiWave[i]=fPhiWaveDeg[i]*TMath::DegToRad(); 00094 } 00095 for(int i=0;i<NUM_THETA_BINS;i++) { 00096 fThetaWaveDeg[i]=-90+0.5*deltaTheta+deltaTheta*i; 00097 fThetaWave[i]=fThetaWaveDeg[i]*TMath::DegToRad(); 00098 } 00099 } 00100 00101 00102 void AraEventCorrelator::setupDeltaTInfinity() 00103 { 00104 00105 00106 for(int pair=0;pair<fNumPairs;pair++) { 00107 int ind1=0; 00108 int ind2=0; 00109 getPairIndices(pair,ind1,ind2); 00110 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00111 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00112 fVPolDeltaTInfinity[pair][phiBin][thetaBin]= 00113 calcDeltaTInfinity(fVPolPos[ind1],fVPolRho[ind1],fVPolPhi[ind1], 00114 fVPolPos[ind2],fVPolRho[ind2],fVPolPhi[ind2], 00115 fPhiWave[phiBin],fThetaWave[thetaBin]); 00116 fHPolDeltaTInfinity[pair][phiBin][thetaBin]= 00117 calcDeltaTInfinity(fHPolPos[ind1],fHPolRho[ind1],fHPolPhi[ind1], 00118 fHPolPos[ind2],fHPolRho[ind2],fHPolPhi[ind2], 00119 fPhiWave[phiBin],fThetaWave[thetaBin]); 00120 } 00121 } 00122 } 00123 00124 } 00125 00126 00127 void AraEventCorrelator::setupDeltaT40m() 00128 { 00129 Double_t R=41.8; 00130 00131 for(int pair=0;pair<fNumPairs;pair++) { 00132 int ind1=0; 00133 int ind2=0; 00134 getPairIndices(pair,ind1,ind2); 00135 // std::cout << "Starting pair " << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00136 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00137 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00138 fVPolDeltaT40m[pair][phiBin][thetaBin]= 00139 calcDeltaTR(fVPolPos[ind1],fVPolRho[ind1],fVPolPhi[ind1], 00140 fVPolPos[ind2],fVPolRho[ind2],fVPolPhi[ind2], 00141 fPhiWave[phiBin],fThetaWave[thetaBin],R); 00142 fHPolDeltaT40m[pair][phiBin][thetaBin]= 00143 calcDeltaTR(fHPolPos[ind1],fHPolRho[ind1],fHPolPhi[ind1], 00144 fHPolPos[ind2],fHPolRho[ind2],fHPolPhi[ind2], 00145 fPhiWave[phiBin],fThetaWave[thetaBin],R); 00146 } 00147 } 00148 } 00149 00150 } 00151 00152 Double_t AraEventCorrelator::calcDeltaTInfinity(Double_t ant1[3],Double_t rho1, Double_t phi1, Double_t ant2[3],Double_t rho2, Double_t phi2, Double_t phiWave, Double_t thetaWave) 00153 { 00154 00155 Double_t d1=TMath::Cos(thetaWave)*(ant1[2]*TMath::Tan(thetaWave)+rho1*TMath::Cos(phi1-phiWave)); 00156 Double_t d2=TMath::Cos(thetaWave)*(ant2[2]*TMath::Tan(thetaWave)+rho2*TMath::Cos(phi2-phiWave)); 00157 Double_t t1t2=(d2-d1)*AraGeomTool::nTopOfIce/TMath::C(); 00158 t1t2*=1e9; 00159 return t1t2; 00160 // Double_t t2t1=(d1-d2)*AraGeomTool::nTopOfIce/TMath::C(); 00161 // t2t1*=1e9; 00162 // return t2t1; 00163 } 00164 00165 Double_t AraEventCorrelator::calcDeltaTR(Double_t ant1[3],Double_t rho1, Double_t phi1, Double_t ant2[3],Double_t rho2, Double_t phi2, Double_t phiWave, Double_t thetaWave, Double_t R) 00166 { 00167 00168 Double_t xs=R*TMath::Cos(thetaWave)*TMath::Cos(phiWave); 00169 Double_t ys=R*TMath::Cos(thetaWave)*TMath::Sin(phiWave); 00170 Double_t zs=R*TMath::Sin(thetaWave); 00171 00172 // if(phiWave*TMath::RadToDeg()>30 && phiWave*TMath::RadToDeg()<31) { 00173 // if(thetaWave*TMath::RadToDeg()>-45 && thetaWave*TMath::RadToDeg()<-44) { 00174 // std::cout << phiWave*TMath::RadToDeg() << "\t" << thetaWave*TMath::RadToDeg() << "\t" << R << "\t" << xs << "\t" << ys << "\t" << zs << "\n"; 00175 // } 00176 // } 00177 00178 Double_t d1=TMath::Sqrt((xs-ant1[0])*(xs-ant1[0])+(ys-ant1[1])*(ys-ant1[1])+(zs-ant1[2])*(zs-ant1[2])); 00179 Double_t d2=TMath::Sqrt((xs-ant2[0])*(xs-ant2[0])+(ys-ant2[1])*(ys-ant2[1])+(zs-ant2[2])*(zs-ant2[2])); 00180 00181 00182 // if(phiWave*TMath::RadToDeg()>30 && phiWave*TMath::RadToDeg()<31) { 00183 // if(thetaWave*TMath::RadToDeg()>-45 && thetaWave*TMath::RadToDeg()<-44) { 00184 // std::cout << d1 << "\t" << d2 << "\t" << d1-d2 << "\n"; 00185 // } 00186 // } 00187 00188 00189 Double_t t1t2=(d1-d2)*AraGeomTool::nTopOfIce/TMath::C(); 00190 t1t2*=1e9; 00191 return t1t2; 00192 } 00193 00194 void AraEventCorrelator::fillDeltaTArrays(AraCorrelatorType::AraCorrelatorType_t corType) { 00195 static AraCorrelatorType::AraCorrelatorType_t lastCorType=AraCorrelatorType::kNotACorType; 00196 if(lastCorType==corType) return; 00197 lastCorType=corType; 00198 00199 switch (corType) { 00200 case AraCorrelatorType::kSphericalDist40: 00201 for(int pair=0;pair<fNumPairs;pair++) { 00202 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00203 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00204 fVPolDeltaT[pair][phiBin][thetaBin]=fVPolDeltaT40m[pair][phiBin][thetaBin]; 00205 fHPolDeltaT[pair][phiBin][thetaBin]=fHPolDeltaT40m[pair][phiBin][thetaBin]; 00206 } 00207 } 00208 } 00209 break; 00210 case AraCorrelatorType::kPlaneWave: 00211 default: 00212 for(int pair=0;pair<fNumPairs;pair++) { 00213 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00214 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00215 fVPolDeltaT[pair][phiBin][thetaBin]=fVPolDeltaTInfinity[pair][phiBin][thetaBin]; 00216 fHPolDeltaT[pair][phiBin][thetaBin]=fHPolDeltaTInfinity[pair][phiBin][thetaBin]; 00217 } 00218 } 00219 } 00220 break; 00221 } 00222 } 00223 00224 00225 00226 void AraEventCorrelator::getPairIndices(int pair, int &ant1, int &ant2) 00227 { 00228 00229 if(pair>=0 && pair<fNumPairs) { 00230 ant1=fFirstAnt[pair]; 00231 ant2=fSecondAnt[pair]; 00232 } 00233 } 00234 00235 TH2D *AraEventCorrelator::getInterferometricMap(UsefulAraEvent *evPtr, AraAntPol::AraAntPol_t polType,AraCorrelatorType::AraCorrelatorType_t corType) 00236 { 00237 static int counter=0; 00238 fillDeltaTArrays(corType); 00239 //Now need to get correlations for the fNumPairs antenna pairs and then look up the correlation values of each one 00240 Double_t scale=1./fNumPairs; 00241 TH2D *histMap = new TH2D("histMap","histMap",NUM_PHI_BINS,-180,180,NUM_THETA_BINS,-90,90); 00242 TGraph *grRaw[fNumAnts]; 00243 TGraph *grInt[fNumAnts]; 00244 TGraph *grNorm[fNumAnts]; 00245 TGraph *grCor[fNumPairs]; 00246 for(int i=0;i<fNumAnts;i++) { 00247 grRaw[i]=0; 00248 grInt[i]=0; 00249 grNorm[i]=0; 00250 } 00251 for(int i=0;i<fNumPairs;i++) 00252 grCor[i]=0; 00253 00254 if(polType==AraAntPol::kVertical) { 00255 for(int ind=0;ind<fNumAnts;ind++) { 00256 // std::cerr << ind << "\t" << fRfChanVPol[ind] << "\n"; 00257 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanVPol[ind]); 00258 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00259 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00260 } 00261 // std::cerr << "Got graphs and made int maps\n"; 00262 00263 for(int pair=0;pair<fNumPairs;pair++) { 00264 int ind1=0; 00265 int ind2=0; 00266 getPairIndices(pair,ind1,ind2); 00267 // std::cerr << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00268 00269 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00270 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00271 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00272 //I think this is the correct equation to work out the bin number 00273 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00274 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00275 Double_t dt=fVPolDeltaT[pair][phiBin][thetaBin]; 00276 //Double_t corVal=grCor[pair]->Eval(dt); 00277 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00278 corVal*=scale; 00279 Double_t binVal=histMap->GetBinContent(globalBin); 00280 histMap->SetBinContent(globalBin,binVal+corVal); 00281 } 00282 } 00283 } 00284 } 00285 else { 00286 for(int ind=0;ind<fNumAnts;ind++) { 00287 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanHPol[ind]); 00288 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00289 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00290 } 00291 00292 for(int pair=0;pair<fNumPairs;pair++) { 00293 int ind1=0; 00294 int ind2=0; 00295 getPairIndices(pair,ind1,ind2); 00296 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00297 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00298 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00299 //I think this is the correct equation to work out the bin number 00300 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00301 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00302 Double_t dt=fHPolDeltaT[pair][phiBin][thetaBin]; 00303 // Double_t corVal=grCor[pair]->Eval(dt); 00304 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00305 corVal*=scale; 00306 Double_t binVal=histMap->GetBinContent(globalBin); 00307 histMap->SetBinContent(globalBin,binVal+corVal); 00308 } 00309 } 00310 } 00311 } 00312 if(fDebugMode) { 00313 char histName[180]; 00314 for(int i=0;i<fNumAnts;i++) { 00315 sprintf(histName,"grRaw%d_%d",i,counter); 00316 grRaw[i]->SetName(histName); 00317 grRaw[i]->SetTitle(histName); 00318 grRaw[i]->Write(); 00319 sprintf(histName,"grInt%d_%d",i,counter); 00320 grInt[i]->SetName(histName); 00321 grInt[i]->SetTitle(histName); 00322 grInt[i]->Write(); 00323 sprintf(histName,"grNorm%d_%d",i,counter); 00324 grNorm[i]->SetName(histName); 00325 grNorm[i]->SetTitle(histName); 00326 grNorm[i]->Write(); 00327 } 00328 for(int i=0;i<fNumPairs;i++) { 00329 sprintf(histName,"grCor%d_%d",i,counter); 00330 grCor[i]->SetName(histName); 00331 grCor[i]->SetTitle(histName); 00332 grCor[i]->Write(); 00333 } 00334 counter++; 00335 } 00336 for(int i=0;i<fNumAnts;i++) { 00337 if(grRaw[i]) delete grRaw[i]; 00338 if(grInt[i]) delete grInt[i]; 00339 if(grNorm[i]) delete grNorm[i]; 00340 } 00341 for(int i=0;i<fNumPairs;i++) { 00342 if(grCor[i]) delete grCor[i]; 00343 } 00344 //Moved the scaling to earlier for optimisation reasons 00345 // histMap->Scale(scale); 00346 return histMap; 00347 } 00348 00349 TGraph* AraEventCorrelator::getNormalisedGraph(TGraph *grIn) 00350 { 00351 Double_t rms=grIn->GetRMS(2); 00352 Double_t mean=grIn->GetMean(2); 00353 Double_t *xVals = grIn->GetX(); 00354 Double_t *yVals = grIn->GetY(); 00355 Int_t numPoints = grIn->GetN(); 00356 Double_t *newY = new Double_t [numPoints]; 00357 for(int i=0;i<numPoints;i++) { 00358 newY[i]=(yVals[i]-mean)/rms; 00359 } 00360 TGraph *grOut = new TGraph(numPoints,xVals,newY); 00361 delete [] newY; 00362 return grOut; 00363 } 00364 00365 Double_t AraEventCorrelator::fastEvalForEvenSampling(TGraph *grIn, Double_t xvalue) 00366 { 00367 Int_t numPoints=grIn->GetN(); 00368 if(numPoints<2) return 0; 00369 Double_t *xVals=grIn->GetX(); 00370 Double_t *yVals=grIn->GetY(); 00371 Double_t dx=xVals[1]-xVals[0]; 00372 if(dx<=0) return 0; 00373 00374 Int_t p0=Int_t((xvalue-xVals[0])/dx); 00375 if(p0<0) p0=0; 00376 if(p0>=numPoints) p0=numPoints-2; 00377 return FFTtools::simpleInterploate(xVals[p0],yVals[p0],xVals[p0+1],yVals[p0+1],xvalue); 00378 }
Generated on Wed Aug 8 16:18:55 2012 for ARA ROOT Test Bed Software by
1.4.7