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