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 00417 //Set the Zaxis scale; 00418 // histMap->GetZaxis()->SetRangeUser(-1,1); 00419 00420 return histMap; 00421 } 00422 00423 TH2D *AraEventCorrelator::getInterferometricMap(UsefulAtriStationEvent *evPtr, AraAntPol::AraAntPol_t polType,AraCorrelatorType::AraCorrelatorType_t corType) 00424 { 00425 static int counter=0; 00426 fillDeltaTArrays(corType); 00427 //Now need to get correlations for the fNumPairs antenna pairs and then look up the correlation values of each one 00428 Double_t scale=1./fNumPairs; 00429 TH2D *histMap = new TH2D("histMap","histMap",NUM_PHI_BINS,-180,180,NUM_THETA_BINS,-90,90); 00430 TGraph *grRaw[fNumAnts]; 00431 TGraph *grInt[fNumAnts]; 00432 TGraph *grNorm[fNumAnts]; 00433 TGraph *grCor[fNumPairs]; 00434 for(int i=0;i<fNumAnts;i++) { 00435 grRaw[i]=0; 00436 grInt[i]=0; 00437 grNorm[i]=0; 00438 } 00439 for(int i=0;i<fNumPairs;i++) 00440 grCor[i]=0; 00441 00442 if(polType==AraAntPol::kVertical) { 00443 for(int ind=0;ind<fNumAnts;ind++) { 00444 std::cerr << ind << "\t" << fRfChanVPol[ind] << "\n"; 00445 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanVPol[ind]); 00446 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00447 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00448 } 00449 std::cerr << "Got graphs and made int maps\n"; 00450 00451 for(int pair=0;pair<fNumPairs;pair++) { 00452 int ind1=0; 00453 int ind2=0; 00454 getPairIndices(pair,ind1,ind2); 00455 std::cerr << pair << "\t" << ind1 << "\t" << ind2 << "\n"; 00456 00457 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00458 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00459 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00460 //I think this is the correct equation to work out the bin number 00461 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00462 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00463 Double_t dt=fVPolDeltaT[pair][phiBin][thetaBin]; 00464 //Double_t corVal=grCor[pair]->Eval(dt); 00465 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00466 corVal*=scale; 00467 Double_t binVal=histMap->GetBinContent(globalBin); 00468 histMap->SetBinContent(globalBin,binVal+corVal); 00469 } 00470 } 00471 } 00472 } 00473 else { 00474 for(int ind=0;ind<fNumAnts;ind++) { 00475 grRaw[ind]=evPtr->getGraphFromRFChan(fRfChanHPol[ind]); 00476 grInt[ind]=FFTtools::getInterpolatedGraph(grRaw[ind],0.5); 00477 grNorm[ind]=getNormalisedGraph(grInt[ind]); 00478 } 00479 00480 for(int pair=0;pair<fNumPairs;pair++) { 00481 int ind1=0; 00482 int ind2=0; 00483 getPairIndices(pair,ind1,ind2); 00484 grCor[pair]=FFTtools::getCorrelationGraph(grNorm[ind1],grNorm[ind2]); 00485 for(int phiBin=0;phiBin<NUM_PHI_BINS;phiBin++) { 00486 for(int thetaBin=0;thetaBin<NUM_THETA_BINS;thetaBin++) { 00487 //I think this is the correct equation to work out the bin number 00488 //Could just use TH2::GetBin(binx,biny) but the below should be faster 00489 Int_t globalBin=(phiBin+1)+(thetaBin+1)*(NUM_PHI_BINS+2); 00490 Double_t dt=fHPolDeltaT[pair][phiBin][thetaBin]; 00491 // Double_t corVal=grCor[pair]->Eval(dt); 00492 Double_t corVal=fastEvalForEvenSampling(grCor[pair],dt); 00493 corVal*=scale; 00494 Double_t binVal=histMap->GetBinContent(globalBin); 00495 histMap->SetBinContent(globalBin,binVal+corVal); 00496 } 00497 } 00498 } 00499 } 00500 if(fDebugMode) { 00501 char histName[180]; 00502 for(int i=0;i<fNumAnts;i++) { 00503 sprintf(histName,"grRaw%d_%d",i,counter); 00504 grRaw[i]->SetName(histName); 00505 grRaw[i]->SetTitle(histName); 00506 grRaw[i]->Write(); 00507 sprintf(histName,"grInt%d_%d",i,counter); 00508 grInt[i]->SetName(histName); 00509 grInt[i]->SetTitle(histName); 00510 grInt[i]->Write(); 00511 sprintf(histName,"grNorm%d_%d",i,counter); 00512 grNorm[i]->SetName(histName); 00513 grNorm[i]->SetTitle(histName); 00514 grNorm[i]->Write(); 00515 } 00516 for(int i=0;i<fNumPairs;i++) { 00517 sprintf(histName,"grCor%d_%d",i,counter); 00518 grCor[i]->SetName(histName); 00519 grCor[i]->SetTitle(histName); 00520 grCor[i]->Write(); 00521 } 00522 counter++; 00523 } 00524 for(int i=0;i<fNumAnts;i++) { 00525 if(grRaw[i]) delete grRaw[i]; 00526 if(grInt[i]) delete grInt[i]; 00527 if(grNorm[i]) delete grNorm[i]; 00528 } 00529 for(int i=0;i<fNumPairs;i++) { 00530 if(grCor[i]) delete grCor[i]; 00531 } 00532 //Moved the scaling to earlier for optimisation reasons 00533 // histMap->Scale(scale); 00534 return histMap; 00535 } 00536 00537 TGraph* AraEventCorrelator::getNormalisedGraph(TGraph *grIn) 00538 { 00539 Double_t rms=grIn->GetRMS(2); 00540 Double_t mean=grIn->GetMean(2); 00541 Double_t *xVals = grIn->GetX(); 00542 Double_t *yVals = grIn->GetY(); 00543 Int_t numPoints = grIn->GetN(); 00544 Double_t *newY = new Double_t [numPoints]; 00545 for(int i=0;i<numPoints;i++) { 00546 newY[i]=(yVals[i]-mean)/rms; 00547 } 00548 TGraph *grOut = new TGraph(numPoints,xVals,newY); 00549 delete [] newY; 00550 return grOut; 00551 } 00552 00553 Double_t AraEventCorrelator::fastEvalForEvenSampling(TGraph *grIn, Double_t xvalue) 00554 { 00555 Int_t numPoints=grIn->GetN(); 00556 if(numPoints<2) return 0; 00557 Double_t *xVals=grIn->GetX(); 00558 Double_t *yVals=grIn->GetY(); 00559 Double_t dx=xVals[1]-xVals[0]; 00560 if(dx<=0) return 0; 00561 00562 Int_t p0=Int_t((xvalue-xVals[0])/dx); 00563 if(p0<0) p0=0; 00564 if(p0>=numPoints) p0=numPoints-2; 00565 return FFTtools::simpleInterploate(xVals[p0],yVals[p0],xVals[p0+1],yVals[p0+1],xvalue); 00566 }
Generated on Mon Dec 9 13:20:21 2013 for ARA ROOT v3.13 Software by
1.4.7