AraCorrelator/AraEventCorrelator.cxx
00001 #include <iostream> 00002 #include "AraEventCorrelator.h" 00003 #include "AraGeomTool.h" 00004 #include "AraAntennaInfo.h" 00005 #include "AraStationInfo.h" 00006 #include "UsefulIcrrStationEvent.h" 00007 #include "UsefulAtriStationEvent.h" 00008 #include "FFTtools.h" 00009 #include "TH2D.h" 00010 #include "TMath.h" 00011 00012 ClassImp(AraEventCorrelator); 00013 00014 AraEventCorrelator * AraEventCorrelator::fgInstance=0; 00015 00016 00017 AraEventCorrelator::AraEventCorrelator(Int_t numAnts, Int_t stationId) 00018 { 00019 //Default constructor 00020 fNumAnts=numAnts; 00021 fNumPairs=0; 00022 fStationId=stationId; 00023 for(int first=0;first<(fNumAnts-1);first++) { 00024 for(int second=first+1;second<fNumAnts;second++) { 00025 fFirstAnt[fNumPairs]=first; 00026 fSecondAnt[fNumPairs]=second; 00027 fNumPairs++; 00028 } 00029 } 00030 00031 // fillDeltaTArrays(); 00032 fillAntennaPositions(stationId); 00033 setupDeltaTInfinity(); 00034 setupDeltaT40m(); 00035 } 00036 00037 AraEventCorrelator::~AraEventCorrelator() 00038 { 00039 //Default destructor 00040 } 00041 00042 00043 //______________________________________________________________________________ 00044 AraEventCorrelator* AraEventCorrelator::Instance(Int_t numAnts, Int_t stationId) 00045 { 00046 //static function 00047 if(fgInstance) 00048 return fgInstance; 00049 00050 fgInstance = new AraEventCorrelator(numAnts, stationId); 00051 return fgInstance; 00052 } 00053 00054 void AraEventCorrelator::fillAntennaPositions(Int_t stationId){ 00055 if(AraGeomTool::isAtriStation(stationId)) fillAntennaPositionsAtri(); 00056 else if(AraGeomTool::isIcrrStation(stationId)) fillAntennaPositionsIcrr(); 00057 else fprintf(stderr, "%s station %i is neither ATRI nor ICRR\n", __FUNCTION__, (int)stationId); 00058 00059 } 00060 00061 void AraEventCorrelator::fillAntennaPositionsAtri() 00062 { 00063 AraGeomTool *araGeom=AraGeomTool::Instance(); 00064 for(int ant=0;ant<ANTS_PER_ATRI;ant++){ 00065 int antPolNum=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->antPolNum; 00066 if(araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->polType==AraAntPol::kVertical) { 00067 if(antPolNum<8) { 00068 fRfChanVPol[antPolNum]=ant; 00069 fVPolPos[antPolNum][0]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[0]; 00070 fVPolPos[antPolNum][1]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[1]; 00071 fVPolPos[antPolNum][2]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[2]; 00072 fVPolRho[antPolNum]=TMath::Sqrt(fVPolPos[antPolNum][0]*fVPolPos[antPolNum][0]+ 00073 fVPolPos[antPolNum][1]*fVPolPos[antPolNum][1]); 00074 fVPolPhi[antPolNum]=TMath::ATan2(fVPolPos[antPolNum][1],fVPolPos[antPolNum][0]); 00075 } 00076 } 00077 if(araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->polType==AraAntPol::kHorizontal) { 00078 if(antPolNum<8) { 00079 fRfChanHPol[antPolNum]=ant; 00080 fHPolPos[antPolNum][0]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[0]; 00081 fHPolPos[antPolNum][1]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[1]; 00082 fHPolPos[antPolNum][2]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[2]; 00083 fHPolRho[antPolNum]=TMath::Sqrt(fHPolPos[antPolNum][0]*fHPolPos[antPolNum][0]+ 00084 fHPolPos[antPolNum][1]*fHPolPos[antPolNum][1]); 00085 fHPolPhi[antPolNum]=TMath::ATan2(fHPolPos[antPolNum][1],fHPolPos[antPolNum][0]); 00086 } 00087 } 00088 } 00089 std::cout << "\n"; 00090 for(int i=0;i<8;i++) { 00091 std::cout << "V\t" << i << "\t" << fVPolPos[i][0] << "\t" << fVPolPos[i][1] << "\t" << fVPolPos[i][2] << "\n"; 00092 } 00093 for(int i=0;i<8;i++) { 00094 std::cout << "H\t" << i << "\t" << fHPolPos[i][0] << "\t" << fHPolPos[i][1] << "\t" << fHPolPos[i][2] << "\n"; 00095 } 00096 std::cout << "\n"; 00097 //Now fill the arrays with angles 00098 Double_t deltaPhi=360./NUM_PHI_BINS; 00099 Double_t deltaTheta=180./NUM_THETA_BINS; 00100 00101 for(int i=0;i<NUM_PHI_BINS;i++) { 00102 fPhiWaveDeg[i]=-180+0.5*deltaPhi+deltaPhi*i; 00103 fPhiWave[i]=fPhiWaveDeg[i]*TMath::DegToRad(); 00104 } 00105 for(int i=0;i<NUM_THETA_BINS;i++) { 00106 fThetaWaveDeg[i]=-90+0.5*deltaTheta+deltaTheta*i; 00107 fThetaWave[i]=fThetaWaveDeg[i]*TMath::DegToRad(); 00108 } 00109 00110 } 00111 00112 void AraEventCorrelator::fillAntennaPositionsIcrr() 00113 { 00114 AraGeomTool *araGeom=AraGeomTool::Instance(); 00115 00116 for(int ant=0;ant<ANTS_PER_ICRR;ant++) { 00117 int antPolNum=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->antPolNum; 00118 std::cerr << ant << "\t" << antPolNum << "\t" << araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->polType << "\n"; 00119 if(araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->polType==AraAntPol::kVertical) { 00120 if(antPolNum<7) { 00121 fRfChanVPol[antPolNum]=ant; 00122 fVPolPos[antPolNum][0]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[0]; 00123 fVPolPos[antPolNum][1]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[1]; 00124 fVPolPos[antPolNum][2]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[2]; 00125 fVPolRho[antPolNum]=TMath::Sqrt(fVPolPos[antPolNum][0]*fVPolPos[antPolNum][0]+ 00126 fVPolPos[antPolNum][1]*fVPolPos[antPolNum][1]); 00127 fVPolPhi[antPolNum]=TMath::ATan2(fVPolPos[antPolNum][1],fVPolPos[antPolNum][0]); 00128 } 00129 } 00130 if(araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->polType==AraAntPol::kHorizontal) { 00131 if(antPolNum<7) { 00132 fRfChanHPol[antPolNum]=ant; 00133 fHPolPos[antPolNum][0]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[0]; 00134 fHPolPos[antPolNum][1]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[1]; 00135 fHPolPos[antPolNum][2]=araGeom->getStationInfo(fStationId)->getAntennaInfo(ant)->getLocationXYZ()[2]; 00136 fHPolRho[antPolNum]=TMath::Sqrt(fHPolPos[antPolNum][0]*fHPolPos[antPolNum][0]+ 00137 fHPolPos[antPolNum][1]*fHPolPos[antPolNum][1]); 00138 fHPolPhi[antPolNum]=TMath::ATan2(fHPolPos[antPolNum][1],fHPolPos[antPolNum][0]); 00139 } 00140 } 00141 } 00142 std::cout << "\n"; 00143 for(int i=0;i<7;i++) { 00144 std::cout << "V\t" << i << "\t" << fVPolPos[i][0] << "\t" << fVPolPos[i][1] << "\t" << fVPolPos[i][2] << "\n"; 00145 } 00146 for(int i=0;i<7;i++) { 00147 std::cout << "H\t" << i << "\t" << fHPolPos[i][0] << "\t" << fHPolPos[i][1] << "\t" << fHPolPos[i][2] << "\n"; 00148 } 00149 std::cout << "\n"; 00150 00151 //Now fill the arrays with angles 00152 Double_t deltaPhi=360./NUM_PHI_BINS; 00153 Double_t deltaTheta=180./NUM_THETA_BINS; 00154 00155 for(int i=0;i<NUM_PHI_BINS;i++) { 00156 fPhiWaveDeg[i]=-180+0.5*deltaPhi+deltaPhi*i; 00157 fPhiWave[i]=fPhiWaveDeg[i]*TMath::DegToRad(); 00158 } 00159 for(int i=0;i<NUM_THETA_BINS;i++) { 00160 fThetaWaveDeg[i]=-90+0.5*deltaTheta+deltaTheta*i; 00161 fThetaWave[i]=fThetaWaveDeg[i]*TMath::DegToRad(); 00162 } 00163 } 00164 00165 00166 void AraEventCorrelator::setupDeltaTInfinity() 00167 { 00168 00169 00170 for(int pair=0;pair<fNumPairs;pair++) { 00171 int ind1=0; 00172 int ind2=0; 00173 getPairIndices(pair,ind1,ind2); 00174 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00175 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00176 fVPolDeltaTInfinity[pair][phiBin][thetaBin]= 00177 calcDeltaTInfinity(fVPolPos[ind1],fVPolRho[ind1],fVPolPhi[ind1], 00178 fVPolPos[ind2],fVPolRho[ind2],fVPolPhi[ind2], 00179 fPhiWave[phiBin],fThetaWave[thetaBin]); 00180 fHPolDeltaTInfinity[pair][phiBin][thetaBin]= 00181 calcDeltaTInfinity(fHPolPos[ind1],fHPolRho[ind1],fHPolPhi[ind1], 00182 fHPolPos[ind2],fHPolRho[ind2],fHPolPhi[ind2], 00183 fPhiWave[phiBin],fThetaWave[thetaBin]); 00184 } 00185 } 00186 } 00187 00188 } 00189 00190 00191 void AraEventCorrelator::setupDeltaT40m() 00192 { 00193 Double_t R=41.8; 00194 00195 00196 for(int pair=0;pair<fNumPairs;pair++) { 00197 int ind1=0; 00198 int ind2=0; 00199 getPairIndices(pair,ind1,ind2); 00200 // std::cout << "Starting pair " << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00201 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00202 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00203 fVPolDeltaT40m[pair][phiBin][thetaBin]= 00204 calcDeltaTR(fVPolPos[ind1],fVPolRho[ind1],fVPolPhi[ind1], 00205 fVPolPos[ind2],fVPolRho[ind2],fVPolPhi[ind2], 00206 fPhiWave[phiBin],fThetaWave[thetaBin],R); 00207 fHPolDeltaT40m[pair][phiBin][thetaBin]= 00208 calcDeltaTR(fHPolPos[ind1],fHPolRho[ind1],fHPolPhi[ind1], 00209 fHPolPos[ind2],fHPolRho[ind2],fHPolPhi[ind2], 00210 fPhiWave[phiBin],fThetaWave[thetaBin],R); 00211 } 00212 } 00213 } 00214 00215 } 00216 00217 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) 00218 { 00219 00220 Double_t d1=TMath::Cos(thetaWave)*(ant1[2]*TMath::Tan(thetaWave)+rho1*TMath::Cos(phi1-phiWave)); 00221 Double_t d2=TMath::Cos(thetaWave)*(ant2[2]*TMath::Tan(thetaWave)+rho2*TMath::Cos(phi2-phiWave)); 00222 Double_t t1t2=(d2-d1)*AraGeomTool::nTopOfIce/TMath::C(); 00223 t1t2*=1e9; 00224 return t1t2; 00225 // Double_t t2t1=(d1-d2)*AraGeomTool::nTopOfIce/TMath::C(); 00226 // t2t1*=1e9; 00227 // return t2t1; 00228 } 00229 00230 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) 00231 { 00232 00233 Double_t xs=R*TMath::Cos(thetaWave)*TMath::Cos(phiWave); 00234 Double_t ys=R*TMath::Cos(thetaWave)*TMath::Sin(phiWave); 00235 Double_t zs=R*TMath::Sin(thetaWave); 00236 00237 // if(phiWave*TMath::RadToDeg()>30 && phiWave*TMath::RadToDeg()<31) { 00238 // if(thetaWave*TMath::RadToDeg()>-45 && thetaWave*TMath::RadToDeg()<-44) { 00239 // std::cout << phiWave*TMath::RadToDeg() << "\t" << thetaWave*TMath::RadToDeg() << "\t" << R << "\t" << xs << "\t" << ys << "\t" << zs << "\n"; 00240 // } 00241 // } 00242 00243 Double_t d1=TMath::Sqrt((xs-ant1[0])*(xs-ant1[0])+(ys-ant1[1])*(ys-ant1[1])+(zs-ant1[2])*(zs-ant1[2])); 00244 Double_t d2=TMath::Sqrt((xs-ant2[0])*(xs-ant2[0])+(ys-ant2[1])*(ys-ant2[1])+(zs-ant2[2])*(zs-ant2[2])); 00245 00246 00247 // if(phiWave*TMath::RadToDeg()>30 && phiWave*TMath::RadToDeg()<31) { 00248 // if(thetaWave*TMath::RadToDeg()>-45 && thetaWave*TMath::RadToDeg()<-44) { 00249 // std::cout << d1 << "\t" << d2 << "\t" << d1-d2 << "\n"; 00250 // } 00251 // } 00252 00253 00254 Double_t t1t2=(d1-d2)*AraGeomTool::nTopOfIce/TMath::C(); 00255 t1t2*=1e9; 00256 return t1t2; 00257 } 00258 00259 void AraEventCorrelator::fillDeltaTArrays(AraCorrelatorType::AraCorrelatorType_t corType) { 00260 static AraCorrelatorType::AraCorrelatorType_t lastCorType=AraCorrelatorType::kNotACorType; 00261 if(lastCorType==corType) return; 00262 lastCorType=corType; 00263 00264 switch (corType) { 00265 case AraCorrelatorType::kSphericalDist40: 00266 for(int pair=0;pair<fNumPairs;pair++) { 00267 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00268 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00269 fVPolDeltaT[pair][phiBin][thetaBin]=fVPolDeltaT40m[pair][phiBin][thetaBin]; 00270 fHPolDeltaT[pair][phiBin][thetaBin]=fHPolDeltaT40m[pair][phiBin][thetaBin]; 00271 } 00272 } 00273 } 00274 break; 00275 case AraCorrelatorType::kPlaneWave: 00276 default: 00277 for(int pair=0;pair<fNumPairs;pair++) { 00278 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00279 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00280 fVPolDeltaT[pair][phiBin][thetaBin]=fVPolDeltaTInfinity[pair][phiBin][thetaBin]; 00281 fHPolDeltaT[pair][phiBin][thetaBin]=fHPolDeltaTInfinity[pair][phiBin][thetaBin]; 00282 } 00283 } 00284 } 00285 break; 00286 } 00287 } 00288 00289 00290 00291 void AraEventCorrelator::getPairIndices(int pair, int &ant1, int &ant2) 00292 { 00293 00294 if(pair>=0 && pair<fNumPairs) { 00295 ant1=fFirstAnt[pair]; 00296 ant2=fSecondAnt[pair]; 00297 } 00298 } 00299 00300 TH2D *AraEventCorrelator::getInterferometricMap(UsefulIcrrStationEvent *evPtr, AraAntPol::AraAntPol_t polType,AraCorrelatorType::AraCorrelatorType_t corType) 00301 { 00302 static int counter=0; 00303 fillDeltaTArrays(corType); 00304 //Now need to get correlations for the fNumPairs antenna pairs and then look up the correlation values of each one 00305 Double_t scale=1./fNumPairs; 00306 TH2D *histMap = new TH2D("histMap","histMap",NUM_PHI_BINS,-180,180,NUM_THETA_BINS,-90,90); 00307 TGraph *grRaw[fNumAnts]; 00308 TGraph *grInt[fNumAnts]; 00309 TGraph *grNorm[fNumAnts]; 00310 TGraph *grCor[fNumPairs]; 00311 for(int i=0;i<fNumAnts;i++) { 00312 grRaw[i]=0; 00313 grInt[i]=0; 00314 grNorm[i]=0; 00315 } 00316 for(int i=0;i<fNumPairs;i++) 00317 grCor[i]=0; 00318 00319 if(polType==AraAntPol::kVertical) { 00320 for(int ind=0;ind<fNumAnts;ind++) { 00321 // std::cerr << ind << "\t" << fRfChanVPol[ind] << "\n"; 00322 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanVPol[ind]); 00323 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00324 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00325 } 00326 // std::cerr << "Got graphs and made int maps\n"; 00327 00328 for(int pair=0;pair<fNumPairs;pair++) { 00329 int ind1=0; 00330 int ind2=0; 00331 getPairIndices(pair,ind1,ind2); 00332 // std::cerr << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00333 00334 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00335 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00336 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00337 //I think this is the correct equation to work out the bin number 00338 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00339 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00340 Double_t dt=fVPolDeltaT[pair][phiBin][thetaBin]; 00341 //Double_t corVal=grCor[pair]->Eval(dt); 00342 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00343 corVal*=scale; 00344 Double_t binVal=histMap->GetBinContent(globalBin); 00345 histMap->SetBinContent(globalBin,binVal+corVal); 00346 } 00347 } 00348 } 00349 } 00350 else { 00351 for(int ind=0;ind<fNumAnts;ind++) { 00352 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanHPol[ind]); 00353 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00354 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00355 } 00356 00357 for(int pair=0;pair<fNumPairs;pair++) { 00358 int ind1=0; 00359 int ind2=0; 00360 getPairIndices(pair,ind1,ind2); 00361 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00362 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00363 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00364 //I think this is the correct equation to work out the bin number 00365 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00366 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00367 Double_t dt=fHPolDeltaT[pair][phiBin][thetaBin]; 00368 // Double_t corVal=grCor[pair]->Eval(dt); 00369 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00370 corVal*=scale; 00371 Double_t binVal=histMap->GetBinContent(globalBin); 00372 histMap->SetBinContent(globalBin,binVal+corVal); 00373 } 00374 } 00375 } 00376 } 00377 if(fDebugMode) { 00378 char histName[180]; 00379 for(int i=0;i<fNumAnts;i++) { 00380 sprintf(histName,"grRaw%d_%d",i,counter); 00381 grRaw[i]->SetName(histName); 00382 grRaw[i]->SetTitle(histName); 00383 grRaw[i]->Write(); 00384 sprintf(histName,"grInt%d_%d",i,counter); 00385 grInt[i]->SetName(histName); 00386 grInt[i]->SetTitle(histName); 00387 grInt[i]->Write(); 00388 sprintf(histName,"grNorm%d_%d",i,counter); 00389 grNorm[i]->SetName(histName); 00390 grNorm[i]->SetTitle(histName); 00391 grNorm[i]->Write(); 00392 } 00393 for(int i=0;i<fNumPairs;i++) { 00394 sprintf(histName,"grCor%d_%d",i,counter); 00395 grCor[i]->SetName(histName); 00396 grCor[i]->SetTitle(histName); 00397 grCor[i]->Write(); 00398 } 00399 counter++; 00400 } 00401 for(int i=0;i<fNumAnts;i++) { 00402 if(grRaw[i]) delete grRaw[i]; 00403 if(grInt[i]) delete grInt[i]; 00404 if(grNorm[i]) delete grNorm[i]; 00405 } 00406 for(int i=0;i<fNumPairs;i++) { 00407 if(grCor[i]) delete grCor[i]; 00408 } 00409 //Moved the scaling to earlier for optimisation reasons 00410 // histMap->Scale(scale); 00411 return histMap; 00412 } 00413 00414 TH2D *AraEventCorrelator::getInterferometricMap(UsefulAtriStationEvent *evPtr, AraAntPol::AraAntPol_t polType,AraCorrelatorType::AraCorrelatorType_t corType) 00415 { 00416 static int counter=0; 00417 fillDeltaTArrays(corType); 00418 //Now need to get correlations for the fNumPairs antenna pairs and then look up the correlation values of each one 00419 Double_t scale=1./fNumPairs; 00420 TH2D *histMap = new TH2D("histMap","histMap",NUM_PHI_BINS,-180,180,NUM_THETA_BINS,-90,90); 00421 TGraph *grRaw[fNumAnts]; 00422 TGraph *grInt[fNumAnts]; 00423 TGraph *grNorm[fNumAnts]; 00424 TGraph *grCor[fNumPairs]; 00425 for(int i=0;i<fNumAnts;i++) { 00426 grRaw[i]=0; 00427 grInt[i]=0; 00428 grNorm[i]=0; 00429 } 00430 for(int i=0;i<fNumPairs;i++) 00431 grCor[i]=0; 00432 00433 if(polType==AraAntPol::kVertical) { 00434 for(int ind=0;ind<fNumAnts;ind++) { 00435 std::cerr << ind << "\t" << fRfChanVPol[ind] << "\n"; 00436 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanVPol[ind]); 00437 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00438 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00439 } 00440 std::cerr << "Got graphs and made int maps\n"; 00441 00442 for(int pair=0;pair<fNumPairs;pair++) { 00443 int ind1=0; 00444 int ind2=0; 00445 getPairIndices(pair,ind1,ind2); 00446 std::cerr << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00447 00448 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00449 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00450 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00451 //I think this is the correct equation to work out the bin number 00452 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00453 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00454 Double_t dt=fVPolDeltaT[pair][phiBin][thetaBin]; 00455 //Double_t corVal=grCor[pair]->Eval(dt); 00456 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00457 corVal*=scale; 00458 Double_t binVal=histMap->GetBinContent(globalBin); 00459 histMap->SetBinContent(globalBin,binVal+corVal); 00460 } 00461 } 00462 } 00463 } 00464 else { 00465 for(int ind=0;ind<fNumAnts;ind++) { 00466 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanHPol[ind]); 00467 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00468 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00469 } 00470 00471 for(int pair=0;pair<fNumPairs;pair++) { 00472 int ind1=0; 00473 int ind2=0; 00474 getPairIndices(pair,ind1,ind2); 00475 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00476 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00477 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00478 //I think this is the correct equation to work out the bin number 00479 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00480 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00481 Double_t dt=fHPolDeltaT[pair][phiBin][thetaBin]; 00482 // Double_t corVal=grCor[pair]->Eval(dt); 00483 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00484 corVal*=scale; 00485 Double_t binVal=histMap->GetBinContent(globalBin); 00486 histMap->SetBinContent(globalBin,binVal+corVal); 00487 } 00488 } 00489 } 00490 } 00491 if(fDebugMode) { 00492 char histName[180]; 00493 for(int i=0;i<fNumAnts;i++) { 00494 sprintf(histName,"grRaw%d_%d",i,counter); 00495 grRaw[i]->SetName(histName); 00496 grRaw[i]->SetTitle(histName); 00497 grRaw[i]->Write(); 00498 sprintf(histName,"grInt%d_%d",i,counter); 00499 grInt[i]->SetName(histName); 00500 grInt[i]->SetTitle(histName); 00501 grInt[i]->Write(); 00502 sprintf(histName,"grNorm%d_%d",i,counter); 00503 grNorm[i]->SetName(histName); 00504 grNorm[i]->SetTitle(histName); 00505 grNorm[i]->Write(); 00506 } 00507 for(int i=0;i<fNumPairs;i++) { 00508 sprintf(histName,"grCor%d_%d",i,counter); 00509 grCor[i]->SetName(histName); 00510 grCor[i]->SetTitle(histName); 00511 grCor[i]->Write(); 00512 } 00513 counter++; 00514 } 00515 for(int i=0;i<fNumAnts;i++) { 00516 if(grRaw[i]) delete grRaw[i]; 00517 if(grInt[i]) delete grInt[i]; 00518 if(grNorm[i]) delete grNorm[i]; 00519 } 00520 for(int i=0;i<fNumPairs;i++) { 00521 if(grCor[i]) delete grCor[i]; 00522 } 00523 //Moved the scaling to earlier for optimisation reasons 00524 // histMap->Scale(scale); 00525 return histMap; 00526 } 00527 00528 TGraph* AraEventCorrelator::getNormalisedGraph(TGraph *grIn) 00529 { 00530 Double_t rms=grIn->GetRMS(2); 00531 Double_t mean=grIn->GetMean(2); 00532 Double_t *xVals = grIn->GetX(); 00533 Double_t *yVals = grIn->GetY(); 00534 Int_t numPoints = grIn->GetN(); 00535 Double_t *newY = new Double_t [numPoints]; 00536 for(int i=0;i<numPoints;i++) { 00537 newY[i]=(yVals[i]-mean)/rms; 00538 } 00539 TGraph *grOut = new TGraph(numPoints,xVals,newY); 00540 delete [] newY; 00541 return grOut; 00542 } 00543 00544 Double_t AraEventCorrelator::fastEvalForEvenSampling(TGraph *grIn, Double_t xvalue) 00545 { 00546 Int_t numPoints=grIn->GetN(); 00547 if(numPoints<2) return 0; 00548 Double_t *xVals=grIn->GetX(); 00549 Double_t *yVals=grIn->GetY(); 00550 Double_t dx=xVals[1]-xVals[0]; 00551 if(dx<=0) return 0; 00552 00553 Int_t p0=Int_t((xvalue-xVals[0])/dx); 00554 if(p0<0) p0=0; 00555 if(p0>=numPoints) p0=numPoints-2; 00556 return FFTtools::simpleInterploate(xVals[p0],yVals[p0],xVals[p0+1],yVals[p0+1],xvalue); 00557 }
Generated on Mon Jun 3 16:10:04 2013 for ARA ROOT v3.9 Software by
1.4.7