Changeset 616
- Timestamp:
- Nov 28, 2011, 3:59:31 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/trunk/mergeAnaFiles.cc
r615 r616 63 63 //-------------------------------------------------------------- 64 64 //Utility functions 65 // template <class T> 66 // void MeanSigma(TArray<T> const & a, T & mean, T & sig){ 67 // if (a.NbDimensions() < 1) 68 // throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !"); 69 // const T * pe; 70 // sa_size_t j,k; 71 // mean=0.; 72 // sig = 0.; 73 // T valok; 74 // if (a.AvgStep() > 0) { // regularly spaced elements 75 // sa_size_t step = a.AvgStep(); 76 // sa_size_t maxx = a.Size()*step; 77 // pe = a.Data(); 78 // for(k=0; k<maxx; k+=step ) { 79 // valok = pe[k]; 80 // mean += valok; sig += valok*valok; 81 // } 82 // } 83 // else { // Non regular data spacing ... 84 // int_4 ka = a.MaxSizeKA(); 85 // sa_size_t step = a.Step(ka); 86 // sa_size_t gpas = a.Size(ka)*step; 87 // sa_size_t naxa = a.Size()/a.Size(ka); 88 // for(j=0; j<naxa; j++) { 89 // pe = a.DataBlock().Begin()+a.Offset(ka,j); 90 // for(k=0; k<gpas; k+=step) { 91 // valok = pe[k]; 92 // mean += valok; sig += valok*valok; 93 // } 94 // } 95 // } 96 // T dsz = (T)(a.Size()); 97 // mean /= dsz; 98 // if (dsz > 1.5) { 99 // sig = sig/dsz - mean*mean; 100 // sig *= (dsz/(dsz-1)); 101 // if (sig >= 0.) sig = sqrt(sig); 102 // } 103 // else sig = 0.; 104 // } 105 65 106 66 107 sa_size_t freqToChan(r_4 f){ … … 79 120 for (sa_size_t ir=0; ir<nr; ir++){ 80 121 TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false); 81 doublemean, sigma;122 r_8 mean, sigma; 82 123 MeanSigma(tmp,mean,sigma); 83 124 vec(ir) = mean; … … 273 314 }//eo loop on slices 274 315 } 316 275 317 //AST 18/11/11 276 318 //------------------------------------------------------- … … 449 491 } 450 492 ifs.close(); 451 452 493 // 453 494 calibFileName = directoryName + "/" … … 621 662 } 622 663 ifs.close(); 623 624 664 //link <cycle> - <calibration coefficient> 625 665 //We cannot rely on identical cycle list of the OFF and ON calibration … … 770 810 771 811 } 812 772 813 //------------------------------------------------------- 773 814 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra … … 812 853 sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5)); 813 854 //Lower and Higher freq. just arround 1420.4MHz Freq. bin to perform mean follow up 814 sa_size_t ch1420Low; 815 sa_size_t ch1420High; 816 if (para.sourceName_ == "Abell85") { 817 ch1420Low = freqToChan(1420.4-0.2); // Abell85 818 ch1420High = freqToChan(1420.4+0.2); 819 } else if (para.sourceName_ == "Abell1205") { 820 ch1420Low = freqToChan(1420.4-0.3); // Abell1205 821 ch1420High = freqToChan(1420.4+0.2); 822 } else if (para.sourceName_ == "Abell2440") { 823 ch1420Low = freqToChan(1420.4); // Abell2440 824 ch1420High = freqToChan(1420.4+0.3); 825 } else { 826 ch1420Low = freqToChan(1420.4-0.2); // Abell85 827 ch1420High = freqToChan(1420.4+0.2); 828 } 855 sa_size_t ch1420Low = freqToChan(1420.4-0.2); 856 sa_size_t ch1420High = freqToChan(1420.4+0.2); 857 829 858 //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up 830 859 sa_size_t ch1420aLow = freqToChan(1418); … … 865 894 string srcLower = tokens2[0]; 866 895 867 868 869 896 PInPersist fin(*iFile); 870 897 vector<string> vec = fin.GetNameTags(); … … 1082 1109 nr = calibBAOfactors_Off_Run_Ch0.NRows(); 1083 1110 for (sa_size_t ir=0; ir<nr; ++ir){ 1084 //cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "1085 //<< calibBAOfactors_Off_Run_Ch0(ir,1) << endl;1111 cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val " 1112 << calibBAOfactors_Off_Run_Ch0(ir,1) << endl; 1086 1113 1087 1114 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)] … … 1134 1161 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl; 1135 1162 } else { 1136 cout << "Cycle NOT calibrated "<< endl;1163 cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl; 1137 1164 } 1138 1165 }//debug … … 1178 1205 //JEC 29/10/11 add ON-OFF/OFF 1179 1206 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data 1180 TMatrix<r_4> aSpecOffFi tltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);1207 TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 1181 1208 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window 1182 medianFiltering(aSpecOff,halfWidth,aSpecOffFi tltered);1209 medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered); 1183 1210 1184 diffOnOffOvOff_noCalib.Div(aSpecOffFi tltered); //division in place1211 diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place 1185 1212 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib; 1186 1213 … … 1191 1218 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib; 1192 1219 1193 //1194 1220 totalNumberCycles++; 1221 1195 1222 //Fill NTuple 1196 1223 xnt[0] = totalNumberCycles; … … 1266 1293 //JEC 18/11/11 follow up the 1400-1420MHz OFF only 1267 1294 TMatrix<r_4> aSpecOffovOff(aSpecOff,false); 1268 aSpecOffovOff.Div(aSpecOffFi tltered);1295 aSpecOffovOff.Div(aSpecOffFiltered); 1269 1296 1270 1297 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS); 1298 // meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib); 1271 1299 meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib); 1272 1300 … … 1275 1303 1276 1304 TMatrix<r_4> aSpecOnovOff(aSpecOn,false); 1277 aSpecOnovOff.Div(aSpecOffFi tltered);1305 aSpecOnovOff.Div(aSpecOffFiltered); 1278 1306 1279 1307 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS); 1308 //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib); 1280 1309 meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib); 1281 1310 … … 1352 1381 1353 1382 tag = "onoffevol"; 1354 fos << PPFNameTag(tag) << onoffevolution; 1383 fos << PPFNameTag(tag) << onoffevolution; 1384 1355 1385 }//end of save 1386 1387 1356 1388 } 1357 1389 //JEC 14/11/11 New meanRawDiffOnOffCycles START … … 1376 1408 TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // ON/Filtered_OFF 1377 1409 TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF 1378 1379 1410 //Tuple only for RAW things to follow 1380 1411 static const int NINFO=11; … … 1417 1448 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { 1418 1449 if (para.debuglev_>90){ 1419 } 1450 cout << "load file <" << *iFile << ">" << endl; 1451 } 1452 1453 vector<string> tokens; 1454 split(*iFile,"_",tokens); 1455 string dateOfRun = tokens[1]; 1456 if (para.debuglev_>90){ 1457 cout << "date <" << dateOfRun << ">" << endl; 1458 } 1459 vector<string> tokens2; 1460 split(tokens[2],".",tokens2); 1461 string srcLower = tokens2[0]; 1462 1420 1463 PInPersist fin(*iFile); 1421 1464 vector<string> vec = fin.GetNameTags(); … … 1447 1490 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e); 1448 1491 if (para.debuglev_>90){ 1449 cout << " sp actra <" << *iSpec << ">" << endl;1492 cout << " spectra <" << *iSpec << ">" << endl; 1450 1493 cout << " cycle " << iSpec->substr(b) << endl; 1451 1494 } … … 1475 1518 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){ 1476 1519 1520 // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440 1521 if ( *ic == 1 && srcLower == "abell1205" ) { 1522 if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) { 1523 cout << "Skipping non-calibrated cycle " << *ic << endl; 1524 continue; 1525 } 1526 } else if ( *ic == 1 && srcLower == "abell2440" ) { 1527 if ( dateOfRun == "20110516" ) { 1528 cout << "Skipping non-calibrated cycle " << *ic << endl; 1529 continue; 1530 } 1531 } else if ( *ic == 3 && srcLower == "abell1205" ) { 1532 if ( dateOfRun == "20110810" ) { 1533 cout << "Skipping non-calibrated cycle " << *ic << endl; 1534 continue; 1535 } 1536 } 1537 1477 1538 string ppftag; 1478 1539 //load ON phase … … 1490 1551 //Perform the difference ON-OFF 1491 1552 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff; 1553 1492 1554 meanDiffONOFF_noCalib += diffOnOff_noCalib; 1493 1555 1494 1556 //JEC 29/10/11 add ON-OFF/OFF 1495 1557 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data 1496 TMatrix<r_4> aSpecOffFi tltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);1558 TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 1497 1559 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window 1498 medianFiltering(aSpecOff,halfWidth,aSpecOffFi tltered);1560 medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered); 1499 1561 1500 diffOnOffOvOff_noCalib.Div(aSpecOffFi tltered); //division in place1562 diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place 1501 1563 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib; 1502 1503 1564 1504 1565 //JEC 15/11/11 add ON/OFF and OFF/OFF 1505 1566 TMatrix<r_4> onOvOff(aSpecOn,false); 1506 onOvOff.Div(aSpecOffFi tltered);1567 onOvOff.Div(aSpecOffFiltered); 1507 1568 meanONovOFF_noCalib += onOvOff; 1508 1569 1509 1570 TMatrix<r_4> offOvOff(aSpecOff,false); 1510 offOvOff.Div(aSpecOffFi tltered);1571 offOvOff.Div(aSpecOffFiltered); 1511 1572 meanOFFovOFF_noCalib += offOvOff; 1512 1573 1513 1574 totalNumberCycles++; 1514 1515 1575 1516 1576 //Fill NTuple 1517 1577 xnt[0] = totalNumberCycles; 1518 1578 1519 1520 1579 //Follow up arround the 1420.4MHz Freq. 1521 1580 TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS); … … 1535 1594 1536 1595 1537 //JEC 25/10/11 follow 14 00-1420MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie1538 TVector<r_4> meanInRange_14 00a1420Freq_noCalib(NUMBER_OF_CHANNELS);1539 meanInRange(diffOnOffOvOff_noCalib,ch14 00,ch1420,meanInRange_1400a1420Freq_noCalib);1540 xnt[5] = meanInRange_14 00a1420Freq_noCalib(0);1541 xnt[6] = meanInRange_14 00a1420Freq_noCalib(1);1542 1543 //JEC 18/11/11 follow up the 14 00-1420MHz OFF only1596 //JEC 25/10/11 follow 1410-1415 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie 1597 TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS); 1598 meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib); 1599 xnt[5] = meanInRange_1410a1415Freq_noCalib(0); 1600 xnt[6] = meanInRange_1410a1415Freq_noCalib(1); 1601 1602 //JEC 18/11/11 follow up the 1410-1415MHz OFF only 1544 1603 TMatrix<r_4> aSpecOffovOff(aSpecOff,false); 1545 aSpecOffovOff.Div(aSpecOffFi tltered);1604 aSpecOffovOff.Div(aSpecOffFiltered); 1546 1605 1547 1606 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS); … … 1552 1611 1553 1612 TMatrix<r_4> aSpecOnovOff(aSpecOn,false); 1554 aSpecOnovOff.Div(aSpecOffFi tltered);1613 aSpecOnovOff.Div(aSpecOffFiltered); 1555 1614 1556 1615 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS); … … 1600 1659 1601 1660 tag = "onoffevol"; 1602 fos << PPFNameTag(tag) << onoffevolution; 1661 fos << PPFNameTag(tag) << onoffevolution; 1603 1662 1604 1663 }//end save 1664 1665 1605 1666 } 1606 1667 //JEC 14/11/11 New meanRawDiffOnOffCycles END … … 1962 2023 } else if (action == "meanCalibBAODiffOnOff") { 1963 2024 meanCalibBAODiffOnOffCycles(); 1964 } else if (action == "calibCoeffNtp") {1965 calibCoeffNtp();1966 2025 } else { 1967 2026 msg = "Unknown action " + action;
Note: See TracChangeset
for help on using the changeset viewer.