Changeset 615 for BAORadio/AmasNancay/trunk
- Timestamp:
- Nov 28, 2011, 3:47:28 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/trunk/mergeAnaFiles.cc
r612 r615 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 elements75 // 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 106 65 107 66 sa_size_t freqToChan(r_4 f){ … … 120 79 for (sa_size_t ir=0; ir<nr; ir++){ 121 80 TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false); 122 r_8mean, sigma;81 double mean, sigma; 123 82 MeanSigma(tmp,mean,sigma); 124 83 vec(ir) = mean; … … 314 273 }//eo loop on slices 315 274 } 275 //AST 18/11/11 276 //------------------------------------------------------- 277 //Make n-tuple with calibration coefficients 278 // 279 void calibCoeffNtp() throw(string) { 280 list<string> listOfFiles; 281 string directoryName; 282 directoryName = para.inPath_ + "/" + para.sourceName_; 283 284 //Make the listing of the directory 285 listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_); 286 287 list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd; 288 iFileEnd = listOfFiles.end(); 289 290 static const int NINFO=19; 291 char* calibTupleName[NINFO]={"date","run","cycle" 292 ,"meancoeffoff0","meancoeffoff1" 293 ,"meancoeffon0","meancoeffon1" 294 ,"coeffoff0","coeffoff1" 295 ,"coeffon0","coeffon1" 296 ,"invmeancoeffoff0","invmeancoeffoff1" 297 ,"invmeancoeffon0","invmeancoeffon1" 298 ,"invcoeffoff0","invcoeffoff1" 299 ,"invcoeffon0","invcoeffon1" 300 }; 301 NTuple calibcoeffs(NINFO,calibTupleName); 302 r_4 xnt[NINFO]; 303 304 int totalNumberRuns=0; //total number of runs 305 int totalNumberCycles=0; //total number of cycles for normalisation 306 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { 307 if (para.debuglev_>90){ 308 cout << "load file <" << *iFile << ">" << endl; 309 } 310 311 vector<string> tokens; 312 split(*iFile,"_",tokens); 313 string dateOfRun = tokens[1]; 314 r_4 rdateOfRun = atoi(dateOfRun.c_str()) - 2011*1.e4; 315 316 if (para.debuglev_>90){ 317 cout << "date <" << dateOfRun << ">" << endl; 318 cout << "rdate <" << rdateOfRun << ">" << endl; 319 } 320 vector<string> tokens2; 321 split(tokens[2],".",tokens2); 322 string srcLower = tokens2[0]; 323 324 PInPersist fin(*iFile); 325 vector<string> vec = fin.GetNameTags(); 326 327 vector<string> modeList; 328 modeList.push_back("On"); 329 modeList.push_back("Off"); 330 vector<string>::const_iterator iMode; 331 332 map<string, list<int> > cycleModeCollect; 333 334 for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) { 335 list<string> listOfSpectra; 336 //Keep only required PPF objects 337 string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 338 std::remove_copy_if( 339 vec.begin(), vec.end(), back_inserter(listOfSpectra), 340 not1(StringMatch(matchstr)) 341 ); 342 343 listOfSpectra.sort(stringCompare); 344 iSpecEnd = listOfSpectra.end(); 345 346 matchstr = "[0-9]+"; 347 //Loop of spectra matrix 348 list<int> listOfCycles; 349 for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){ 350 int b,e; 351 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e); 352 if (para.debuglev_>90){ 353 cout << " spectra <" << *iSpec << ">" << endl; 354 cout << " cycle " << iSpec->substr(b) << endl; 355 } 356 listOfCycles.push_back(atoi((iSpec->substr(b)).c_str())); 357 }//end loop spectra 358 cycleModeCollect[*iMode] = listOfCycles; 359 }//end of mode 360 361 //Take the Intersection of the list Of cycles in mode Off and On 362 list<int> commonCycles; 363 set_intersection(cycleModeCollect["On"].begin(), 364 cycleModeCollect["On"].end(), 365 cycleModeCollect["Off"].begin(), 366 cycleModeCollect["Off"].end(), 367 back_inserter(commonCycles) 368 ); 369 370 if (para.debuglev_>90){ 371 cout << "Liste of cycles common to On & Off: <"; 372 for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){ 373 cout << *i << " "; 374 } 375 cout << ">" << endl; 376 } 377 378 // 379 //Load BAO Calibration factors "per Cycle and Channels" 380 //Compute the mean per Cycle to 381 // fill factors "per Run and Channels" with the same cycle list 382 // 383 // 384 //TODO improve the code.... 385 386 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0; 387 TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1; 388 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0; 389 TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1; 390 391 string calibFileName; 392 ifstream ifs; 393 sa_size_t nr,nc; //values read 394 395 //OFF Cycle per Channel 396 calibFileName = directoryName + "/" 397 + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 398 + para.calibFreq_ +"MHz-Ch0Cycles.txt"; 399 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl; 400 ifs.open(calibFileName.c_str()); 401 if ( ! ifs.is_open() ) { 402 throw calibFileName + " cannot be opened..."; 403 } 404 calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc); 405 if(para.debuglev_>9){ 406 cout << "(nr,nc): "<< nr << "," << nc << endl; 407 calibBAOfactors_Off_Cycle_Ch0.Print(cout); 408 cout << endl; 409 } 410 411 TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc); 412 calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0); 413 {//Compute the mean 414 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false); 415 double mean,sigma; 416 MeanSigma(coef,mean,sigma); 417 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean); 418 } 419 if(para.debuglev_>9){ 420 cout << "Fill calib. with mean value " << endl; 421 calibBAOfactors_Off_Run_Ch0.Print(cout); 422 cout << endl; 423 } 424 425 TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch0(nr,nc); 426 invcalibBAOfactors_Off_Cycle_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0); 427 {//Compute the inverse value 428 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false); 429 invcalibBAOfactors_Off_Cycle_Ch0.Column(1).SetCst(1); 430 invcalibBAOfactors_Off_Cycle_Ch0.Column(1).Div(coef); 431 if(para.debuglev_>9){ 432 cout << "Fill calib. with inverse value " << endl; 433 invcalibBAOfactors_Off_Cycle_Ch0.Print(cout); 434 cout << endl; 435 } 436 } 437 TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch0(nr,nc); 438 invcalibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0); 439 {//Compute the inverse mean 440 TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false); 441 double mean,sigma; 442 MeanSigma(coef,mean,sigma); 443 invcalibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean); 444 } 445 if(para.debuglev_>9){ 446 cout << "Fill calib. with inverse mean value " << endl; 447 invcalibBAOfactors_Off_Run_Ch0.Print(cout); 448 cout << endl; 449 } 450 ifs.close(); 451 452 // 453 calibFileName = directoryName + "/" 454 + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 455 + para.calibFreq_ +"MHz-Ch1Cycles.txt"; 456 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl; 457 ifs.open(calibFileName.c_str()); 458 if ( ! ifs.is_open() ) { 459 460 throw calibFileName + " cannot be opened..."; 461 } 462 calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc); 463 if(para.debuglev_>9){ 464 cout << "(nr,nc): "<< nr << "," << nc << endl; 465 calibBAOfactors_Off_Cycle_Ch1.Print(cout); 466 cout << endl; 467 } 468 TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc); 469 calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0); 470 {//Compute the mean 471 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false); 472 double mean,sigma; 473 MeanSigma(coef,mean,sigma); 474 // cout << "Mean: " << mean << " sigma " << sigma << endl; 475 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean); 476 } 477 if(para.debuglev_>9){ 478 cout << "Fill calib. with mean value " << endl; 479 calibBAOfactors_Off_Run_Ch1.Print(cout); 480 cout << endl; 481 } 482 TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch1(nr,nc); 483 invcalibBAOfactors_Off_Cycle_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0); 484 {//Compute the inverse value 485 TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false); 486 invcalibBAOfactors_Off_Cycle_Ch1.Column(1).SetCst(1); 487 invcalibBAOfactors_Off_Cycle_Ch1.Column(1).Div(coef); 488 if(para.debuglev_>9){ 489 cout << "Fill calib. with inverse value " << endl; 490 invcalibBAOfactors_Off_Cycle_Ch1.Print(cout); 491 cout << endl; 492 } 493 } 494 TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch1(nr,nc); 495 invcalibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0); 496 {//Compute the inverse mean 497 TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false); 498 double mean,sigma; 499 MeanSigma(coef,mean,sigma); 500 invcalibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean); 501 } 502 if(para.debuglev_>9){ 503 cout << "Fill calib. with inverse mean value " << endl; 504 invcalibBAOfactors_Off_Run_Ch1.Print(cout); 505 cout << endl; 506 } 507 ifs.close(); 508 509 //ON Cycle per Channel 510 calibFileName = directoryName + "/" 511 + "calib_" + dateOfRun + "_" + srcLower + "_On_" 512 + para.calibFreq_ +"MHz-Ch0Cycles.txt"; 513 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl; 514 ifs.open(calibFileName.c_str()); 515 if ( ! ifs.is_open() ) { 516 517 throw calibFileName + " cannot be opened..."; 518 } 519 calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc); 520 if(para.debuglev_>9){ 521 cout << "(nr,nc): "<< nr << "," << nc << endl; 522 calibBAOfactors_On_Cycle_Ch0.Print(cout); 523 cout << endl; 524 } 525 526 TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc); 527 calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0); 528 {//Compute the mean 529 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false); 530 double mean,sigma; 531 MeanSigma(coef,mean,sigma); 532 // cout << "Mean: " << mean << " sigma " << sigma << endl; 533 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean); 534 } 535 if(para.debuglev_>9){ 536 cout << "Fill calib. with mean value " << endl; 537 calibBAOfactors_On_Run_Ch0.Print(cout); 538 cout << endl; 539 } 540 541 TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch0(nr,nc); 542 invcalibBAOfactors_On_Cycle_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0); 543 {//Compute the inverse value 544 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false); 545 invcalibBAOfactors_On_Cycle_Ch0.Column(1).SetCst(1); 546 invcalibBAOfactors_On_Cycle_Ch0.Column(1).Div(coef); 547 if(para.debuglev_>9){ 548 cout << "Fill calib. with inverse value " << endl; 549 invcalibBAOfactors_On_Cycle_Ch0.Print(cout); 550 cout << endl; 551 } 552 } 553 TMatrix<r_4> invcalibBAOfactors_On_Run_Ch0(nr,nc); 554 invcalibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0); 555 {//Compute the inverse mean 556 TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false); 557 double mean,sigma; 558 MeanSigma(coef,mean,sigma); 559 invcalibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean); 560 } 561 if(para.debuglev_>9){ 562 cout << "Fill calib. with inverse mean value " << endl; 563 invcalibBAOfactors_On_Run_Ch0.Print(cout); 564 cout << endl; 565 } 566 ifs.close(); 567 568 calibFileName = directoryName + "/" 569 + "calib_" + dateOfRun + "_" + srcLower + "_On_" 570 + para.calibFreq_ +"MHz-Ch1Cycles.txt"; 571 if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl; 572 ifs.open(calibFileName.c_str()); 573 if ( ! ifs.is_open() ) { 574 throw calibFileName + " cannot be opened..."; 575 } 576 calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc); 577 if(para.debuglev_>9){ 578 cout << "(nr,nc): "<< nr << "," << nc << endl; 579 calibBAOfactors_On_Cycle_Ch1.Print(cout); 580 cout << endl; 581 } 582 TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc); 583 calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0); 584 {//Compute the mean 585 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false); 586 double mean,sigma; 587 MeanSigma(coef,mean,sigma); 588 // cout << "Mean: " << mean << " sigma " << sigma << endl; 589 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean); 590 } 591 if(para.debuglev_>9){ 592 cout << "Fill calib. with mean value " << endl; 593 calibBAOfactors_On_Run_Ch1.Print(cout); 594 cout << endl; 595 } 596 597 TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch1(nr,nc); 598 invcalibBAOfactors_On_Cycle_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0); 599 {//Compute the inverse value 600 TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false); 601 invcalibBAOfactors_On_Cycle_Ch1.Column(1).SetCst(1); 602 invcalibBAOfactors_On_Cycle_Ch1.Column(1).Div(coef); 603 if(para.debuglev_>9){ 604 cout << "Fill calib. with inverse value " << endl; 605 invcalibBAOfactors_On_Cycle_Ch1.Print(cout); 606 cout << endl; 607 } 608 } 609 TMatrix<r_4> invcalibBAOfactors_On_Run_Ch1(nr,nc); 610 invcalibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0); 611 {//Compute the inverse mean 612 TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false); 613 double mean,sigma; 614 MeanSigma(coef,mean,sigma); 615 invcalibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean); 616 } 617 if(para.debuglev_>9){ 618 cout << "Fill calib. with inverse mean value " << endl; 619 invcalibBAOfactors_On_Run_Ch1.Print(cout); 620 cout << endl; 621 } 622 ifs.close(); 623 624 //link <cycle> - <calibration coefficient> 625 //We cannot rely on identical cycle list of the OFF and ON calibration 626 map<int,r_4> calibBAO_Off_Run_Ch0; 627 map<int,r_4> calibBAO_Off_Run_Ch1; 628 map<int,r_4> calibBAO_On_Run_Ch0; 629 map<int,r_4> calibBAO_On_Run_Ch1; 630 631 map<int,r_4> calibBAO_Off_Cycle_Ch0; 632 map<int,r_4> calibBAO_Off_Cycle_Ch1; 633 map<int,r_4> calibBAO_On_Cycle_Ch0; 634 map<int,r_4> calibBAO_On_Cycle_Ch1; 635 636 map<int,r_4> invcalibBAO_Off_Run_Ch0; 637 map<int,r_4> invcalibBAO_Off_Run_Ch1; 638 map<int,r_4> invcalibBAO_On_Run_Ch0; 639 map<int,r_4> invcalibBAO_On_Run_Ch1; 640 641 map<int,r_4> invcalibBAO_Off_Cycle_Ch0; 642 map<int,r_4> invcalibBAO_Off_Cycle_Ch1; 643 map<int,r_4> invcalibBAO_On_Cycle_Ch0; 644 map<int,r_4> invcalibBAO_On_Cycle_Ch1; 645 646 //per Run based BAO coefficients 647 nr = calibBAOfactors_Off_Run_Ch0.NRows(); 648 for (sa_size_t ir=0; ir<nr; ++ir){ 649 // cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val " 650 // << calibBAOfactors_Off_Run_Ch0(ir,1) << endl; 651 652 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)] 653 = calibBAOfactors_Off_Run_Ch0(ir,1); 654 calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)] 655 = calibBAOfactors_Off_Cycle_Ch0(ir,1); 656 calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)] 657 = calibBAOfactors_Off_Run_Ch1(ir,1); 658 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)] 659 = calibBAOfactors_Off_Cycle_Ch1(ir,1); 660 661 invcalibBAO_Off_Run_Ch0[(int)invcalibBAOfactors_Off_Run_Ch0(ir,0)] 662 = invcalibBAOfactors_Off_Run_Ch0(ir,1); 663 invcalibBAO_Off_Cycle_Ch0[(int)invcalibBAOfactors_Off_Cycle_Ch0(ir,0)] 664 = invcalibBAOfactors_Off_Cycle_Ch0(ir,1); 665 invcalibBAO_Off_Run_Ch1[(int)invcalibBAOfactors_Off_Run_Ch1(ir,0)] 666 = invcalibBAOfactors_Off_Run_Ch1(ir,1); 667 invcalibBAO_Off_Cycle_Ch1[(int)invcalibBAOfactors_Off_Cycle_Ch1(ir,0)] 668 = invcalibBAOfactors_Off_Cycle_Ch1(ir,1); 669 }//eo loop on coef Off 670 671 nr = calibBAOfactors_On_Run_Ch0.NRows(); 672 for (sa_size_t ir=0; ir<nr; ++ir){ 673 calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)] 674 = calibBAOfactors_On_Run_Ch0(ir,1); 675 calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)] 676 = calibBAOfactors_On_Cycle_Ch0(ir,1); 677 calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)] 678 = calibBAOfactors_On_Run_Ch1(ir,1); 679 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)] 680 = calibBAOfactors_On_Cycle_Ch1(ir,1); 681 682 invcalibBAO_On_Run_Ch0[(int)invcalibBAOfactors_On_Run_Ch0(ir,0)] 683 = invcalibBAOfactors_On_Run_Ch0(ir,1); 684 invcalibBAO_On_Cycle_Ch0[(int)invcalibBAOfactors_On_Cycle_Ch0(ir,0)] 685 = invcalibBAOfactors_On_Cycle_Ch0(ir,1); 686 invcalibBAO_On_Run_Ch1[(int)invcalibBAOfactors_On_Run_Ch1(ir,0)] 687 = invcalibBAOfactors_On_Run_Ch1(ir,1); 688 invcalibBAO_On_Cycle_Ch1[(int)invcalibBAOfactors_On_Cycle_Ch1(ir,0)] 689 = invcalibBAOfactors_On_Cycle_Ch1(ir,1); 690 }//eo loop on coef On 691 692 //Loop on cyles 693 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){ 694 695 //Look if the cycle has been calibrated... 696 bool isCycleCalibrated = 697 ( calibBAO_On_Run_Ch0.count(*ic) * 698 calibBAO_On_Run_Ch1.count(*ic) * 699 calibBAO_Off_Run_Ch0.count(*ic) * 700 calibBAO_Off_Run_Ch1.count(*ic) * 701 calibBAO_On_Cycle_Ch0.count(*ic) * 702 calibBAO_On_Cycle_Ch1.count(*ic) * 703 calibBAO_Off_Cycle_Ch0.count(*ic) * 704 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false; 705 706 if(para.debuglev_>9) { 707 cout << "Calibration coefficients for cycle "<<*ic << endl; 708 if (isCycleCalibrated) { 709 cout << "Cycle calibrated " << endl; 710 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " " 711 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n" 712 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " " 713 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n" 714 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " " 715 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n" 716 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " " 717 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl; 718 } else { 719 cout << "Cycle NOT calibrated " << endl; 720 } 721 }//debug 722 723 if ( ! isCycleCalibrated ) continue; 724 725 totalNumberCycles++; 726 //Fill NTuple 727 xnt[0] = rdateOfRun; 728 xnt[1] = totalNumberRuns; 729 xnt[2] = totalNumberCycles; 730 xnt[3] = calibBAO_Off_Run_Ch0[*ic]; 731 xnt[4] = calibBAO_Off_Run_Ch1[*ic]; 732 xnt[5] = calibBAO_On_Run_Ch0[*ic]; 733 xnt[6] = calibBAO_On_Run_Ch1[*ic]; 734 xnt[7] = calibBAO_Off_Cycle_Ch0[*ic]; 735 xnt[8] = calibBAO_Off_Cycle_Ch1[*ic]; 736 xnt[9] = calibBAO_On_Cycle_Ch0[*ic]; 737 xnt[10] = calibBAO_On_Cycle_Ch1[*ic]; 738 xnt[11] = invcalibBAO_Off_Run_Ch0[*ic]; 739 xnt[12] = invcalibBAO_Off_Run_Ch1[*ic]; 740 xnt[13] = invcalibBAO_On_Run_Ch0[*ic]; 741 xnt[14] = invcalibBAO_On_Run_Ch1[*ic]; 742 xnt[15] = invcalibBAO_Off_Cycle_Ch0[*ic]; 743 xnt[16] = invcalibBAO_Off_Cycle_Ch1[*ic]; 744 xnt[17] = invcalibBAO_On_Cycle_Ch0[*ic]; 745 xnt[18] = invcalibBAO_On_Cycle_Ch1[*ic]; 746 747 calibcoeffs.Fill(xnt); 748 749 //Quit if enough 750 if (totalNumberCycles >= para.maxNumberCycles_) break; 751 752 }//eo loop on cycles 753 if (totalNumberCycles >= para.maxNumberCycles_) break; 754 totalNumberRuns++; 755 }//eo loop on spectra in a file 756 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl; 757 758 {//save the results 759 stringstream tmp; 760 tmp << totalNumberCycles; 761 string fileName = para.outPath_+"/calibcoeffTuple_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf"; 762 763 POutPersist fos(fileName); 764 cout << "Save output in " << fileName << endl; 765 766 string tag = "calib"; 767 fos << PPFNameTag(tag) << calibcoeffs; 768 769 }//end of save 770 771 } 316 772 //------------------------------------------------------- 317 773 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra … … 356 812 sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5)); 357 813 //Lower and Higher freq. just arround 1420.4MHz Freq. bin to perform mean follow up 358 sa_size_t ch1420Low = freqToChan(1420.4-0.2); 359 sa_size_t ch1420High = freqToChan(1420.4+0.2); 360 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 } 361 829 //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up 362 830 sa_size_t ch1420aLow = freqToChan(1418); … … 397 865 string srcLower = tokens2[0]; 398 866 867 868 399 869 PInPersist fin(*iFile); 400 870 vector<string> vec = fin.GetNameTags(); … … 612 1082 nr = calibBAOfactors_Off_Run_Ch0.NRows(); 613 1083 for (sa_size_t ir=0; ir<nr; ++ir){ 614 cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "615 << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;1084 // cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val " 1085 // << calibBAOfactors_Off_Run_Ch0(ir,1) << endl; 616 1086 617 1087 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)] … … 664 1134 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl; 665 1135 } else { 666 cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile<< endl;1136 cout << "Cycle NOT calibrated " << endl; 667 1137 } 668 1138 }//debug … … 708 1178 //JEC 29/10/11 add ON-OFF/OFF 709 1179 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data 710 TMatrix<r_4> aSpecOffFi ltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);1180 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 711 1181 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window 712 medianFiltering(aSpecOff,halfWidth,aSpecOffFi ltered);1182 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered); 713 1183 714 diffOnOffOvOff_noCalib.Div(aSpecOffFi ltered); //division in place1184 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place 715 1185 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib; 716 1186 … … 721 1191 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib; 722 1192 1193 // 723 1194 totalNumberCycles++; 724 725 1195 //Fill NTuple 726 1196 xnt[0] = totalNumberCycles; … … 796 1266 //JEC 18/11/11 follow up the 1400-1420MHz OFF only 797 1267 TMatrix<r_4> aSpecOffovOff(aSpecOff,false); 798 aSpecOffovOff.Div(aSpecOffFi ltered);1268 aSpecOffovOff.Div(aSpecOffFitltered); 799 1269 800 1270 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS); 801 // meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);802 1271 meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib); 803 1272 … … 806 1275 807 1276 TMatrix<r_4> aSpecOnovOff(aSpecOn,false); 808 aSpecOnovOff.Div(aSpecOffFi ltered);1277 aSpecOnovOff.Div(aSpecOffFitltered); 809 1278 810 1279 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS); 811 //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);812 1280 meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib); 813 1281 … … 884 1352 885 1353 tag = "onoffevol"; 886 fos << PPFNameTag(tag) << onoffevolution; 887 1354 fos << PPFNameTag(tag) << onoffevolution; 888 1355 }//end of save 889 890 891 1356 } 892 1357 //JEC 14/11/11 New meanRawDiffOnOffCycles START … … 911 1376 TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // ON/Filtered_OFF 912 1377 TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF 1378 913 1379 //Tuple only for RAW things to follow 914 1380 static const int NINFO=11; … … 951 1417 for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { 952 1418 if (para.debuglev_>90){ 953 cout << "load file <" << *iFile << ">" << endl; 954 } 955 956 vector<string> tokens; 957 split(*iFile,"_",tokens); 958 string dateOfRun = tokens[1]; 959 if (para.debuglev_>90){ 960 cout << "date <" << dateOfRun << ">" << endl; 961 } 962 vector<string> tokens2; 963 split(tokens[2],".",tokens2); 964 string srcLower = tokens2[0]; 965 1419 } 966 1420 PInPersist fin(*iFile); 967 1421 vector<string> vec = fin.GetNameTags(); … … 993 1447 regexp(iSpec->c_str(),matchstr.c_str(),&b,&e); 994 1448 if (para.debuglev_>90){ 995 cout << " sp ectra <" << *iSpec << ">" << endl;1449 cout << " spactra <" << *iSpec << ">" << endl; 996 1450 cout << " cycle " << iSpec->substr(b) << endl; 997 1451 } … … 1021 1475 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){ 1022 1476 1023 // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell24401024 if ( *ic == 1 && srcLower == "abell1205" ) {1025 if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {1026 cout << "Skipping non-calibrated cycle " << *ic << endl;1027 continue;1028 }1029 } else if ( *ic == 1 && srcLower == "abell2440" ) {1030 if ( dateOfRun == "20110516" ) {1031 cout << "Skipping non-calibrated cycle " << *ic << endl;1032 continue;1033 }1034 } else if ( *ic == 3 && srcLower == "abell1205" ) {1035 if ( dateOfRun == "20110810" ) {1036 cout << "Skipping non-calibrated cycle " << *ic << endl;1037 continue;1038 }1039 }1040 1041 1477 string ppftag; 1042 1478 //load ON phase … … 1054 1490 //Perform the difference ON-OFF 1055 1491 TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff; 1056 1057 1492 meanDiffONOFF_noCalib += diffOnOff_noCalib; 1058 1493 1059 1494 //JEC 29/10/11 add ON-OFF/OFF 1060 1495 TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data 1061 TMatrix<r_4> aSpecOffFi ltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);1496 TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 1062 1497 sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window 1063 medianFiltering(aSpecOff,halfWidth,aSpecOffFi ltered);1498 medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered); 1064 1499 1065 diffOnOffOvOff_noCalib.Div(aSpecOffFi ltered); //division in place1500 diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place 1066 1501 meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib; 1502 1067 1503 1068 1504 //JEC 15/11/11 add ON/OFF and OFF/OFF 1069 1505 TMatrix<r_4> onOvOff(aSpecOn,false); 1070 onOvOff.Div(aSpecOffFi ltered);1506 onOvOff.Div(aSpecOffFitltered); 1071 1507 meanONovOFF_noCalib += onOvOff; 1072 1508 1073 1509 TMatrix<r_4> offOvOff(aSpecOff,false); 1074 offOvOff.Div(aSpecOffFi ltered);1510 offOvOff.Div(aSpecOffFitltered); 1075 1511 meanOFFovOFF_noCalib += offOvOff; 1076 1512 1077 1513 totalNumberCycles++; 1514 1078 1515 1079 1516 //Fill NTuple 1080 1517 xnt[0] = totalNumberCycles; 1081 1518 1519 1082 1520 //Follow up arround the 1420.4MHz Freq. 1083 1521 TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS); … … 1097 1535 1098 1536 1099 //JEC 25/10/11 follow 14 10-1415MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie1100 TVector<r_4> meanInRange_14 10a1415Freq_noCalib(NUMBER_OF_CHANNELS);1101 meanInRange(diffOnOffOvOff_noCalib,ch14 10,ch1415,meanInRange_1410a1415Freq_noCalib);1102 xnt[5] = meanInRange_14 10a1415Freq_noCalib(0);1103 xnt[6] = meanInRange_14 10a1415Freq_noCalib(1);1104 1105 //JEC 18/11/11 follow up the 14 10-1415MHz OFF only1537 //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie 1538 TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS); 1539 meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib); 1540 xnt[5] = meanInRange_1400a1420Freq_noCalib(0); 1541 xnt[6] = meanInRange_1400a1420Freq_noCalib(1); 1542 1543 //JEC 18/11/11 follow up the 1400-1420MHz OFF only 1106 1544 TMatrix<r_4> aSpecOffovOff(aSpecOff,false); 1107 aSpecOffovOff.Div(aSpecOffFi ltered);1545 aSpecOffovOff.Div(aSpecOffFitltered); 1108 1546 1109 1547 TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS); … … 1114 1552 1115 1553 TMatrix<r_4> aSpecOnovOff(aSpecOn,false); 1116 aSpecOnovOff.Div(aSpecOffFi ltered);1554 aSpecOnovOff.Div(aSpecOffFitltered); 1117 1555 1118 1556 TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS); … … 1162 1600 1163 1601 tag = "onoffevol"; 1164 fos << PPFNameTag(tag) << onoffevolution; 1602 fos << PPFNameTag(tag) << onoffevolution; 1165 1603 1166 1604 }//end save 1167 1168 1169 1605 } 1170 1606 //JEC 14/11/11 New meanRawDiffOnOffCycles END … … 1526 1962 } else if (action == "meanCalibBAODiffOnOff") { 1527 1963 meanCalibBAODiffOnOffCycles(); 1964 } else if (action == "calibCoeffNtp") { 1965 calibCoeffNtp(); 1528 1966 } else { 1529 1967 msg = "Unknown action " + action;
Note: See TracChangeset
for help on using the changeset viewer.