/[PAMELA software]/quicklook/OrbitalRate/src/OrbitalRate.cpp
ViewVC logotype

Diff of /quicklook/OrbitalRate/src/OrbitalRate.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by pam-rm2, Wed Dec 6 15:52:13 2006 UTC revision 1.4 by pam-rm2, Thu Dec 7 13:27:30 2006 UTC
# Line 346  void Rate(TString *filename, TString out Line 346  void Rate(TString *filename, TString out
346    
347    /********** Magnetic Field **************/    /********** Magnetic Field **************/
348    // Check that all this is correct!    // Check that all this is correct!
349    float dimo = 0.0; // dipole moment (computed from dat files)    float br, btheta, bphi;
   float bnorth, beast, bdown, babs;  
   float xl; // L value  
   float icode; // code value for L accuracy (see fortran code)  
   float bab1; // What's  the difference with babs?  
   float stps = 0.005; // step size for field line tracing  
   float bdel = 0.01; // required accuracy  
   float bequ;  // equatorial b value (also called b_0)  
   bool value = 0; // false if bequ is not the minimum b value  
   float rr0; // equatorial radius normalized to earth radius  
   
   // Initialize fortran routines!!!  
   initize_();  
350    
351    // I can now compute the magnetic dipole moment at the actual date,    // I can now compute the magnetic dipole moment at the actual date,
352    // using the cJulian date.  I don't to recompute it for every event    // using the cJulian date.  I don't to recompute it for every event
# Line 368  void Rate(TString *filename, TString out Line 356  void Rate(TString *filename, TString out
356    Int_t d = tledate.GetDay();    Int_t d = tledate.GetDay();
357    float year = (float) y + (m*31+d)/365;    float year = (float) y + (m*31+d)/365;
358    
359    // Compute the magnetic dipole moment.    // Initialize common data for geopack
360    feldcof_(&year, &dimo);    recalc_(y, m*31+d, 0, 0, 0);
361    /********** Magnetic Field **************/    /********** Magnetic Field **************/
362    
363    tr = (TTree*)rootFile->Get("Physics");    tr = (TTree*)rootFile->Get("Physics");
# Line 469  void Rate(TString *filename, TString out Line 457  void Rate(TString *filename, TString out
457        alt = coo.m_Alt;        alt = coo.m_Alt;
458    
459        /********** Magnetic Field **************/        /********** Magnetic Field **************/
460        feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);        igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi);
461        shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);        //      cout<<"("<<(coo.m_Alt+6371.2)/6371.2<<", "<<M_PI/2.-coo.m_Lat<<", "<<coo.m_Lon<<")"<<endl;
       findb0_(&stps, &bdel, &value, &bequ, &rr0);  
462        /********** Magnetic Field **************/        /********** Magnetic Field **************/
463    
464        // serve fare il controllo deltatime < 1?        // serve fare il controllo deltatime < 1?
465        if (deltaTime > 1) cout << endl << "******** deltaTime<1 ********" << endl;        if (deltaTime > 1) cout << endl << "******** deltaTime<1 ********" << endl;
466        // Does nothing for the first two events or if acquisition time if more        // Does nothing for the first two events or if acquisition time if more
467        // than 1s.        // than 1s.
468        if(i<2 || (deltaTime > 1)) continue;        if(i<1 || (deltaTime > 1)) continue;
469    
470        // CAS3 and CAS4 are not rates but only counters.  So I fill        // CAS3 and CAS4 are not rates but only counters.  So I fill
471        // with the bin with the difference beetween the actual counter        // with the bin with the difference beetween the actual counter
# Line 498  void Rate(TString *filename, TString out Line 485  void Rate(TString *filename, TString out
485        // this values but I need to count how many times I fill        // this values but I need to count how many times I fill
486        // each bin.  This is done by the histogram event_counter.        // each bin.  This is done by the histogram event_counter.
487        // I will normalize later.        // I will normalize later.
488        hbabs_counter->Fill(lon, lat, babs);        hbabs_counter->Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5);
489        hbnorth_counter->Fill(lon, lat, bnorth);        hbnorth_counter->Fill(lon, lat, -btheta*1e-5);
490        hbdown_counter->Fill(lon, lat, bdown);        hbdown_counter->Fill(lon, lat, -br*1e-5);
491        hbeast_counter->Fill(lon, lat, beast);        hbeast_counter->Fill(lon, lat, bphi*1e-5);
       hb0_counter->Fill(lon, lat, bequ);  
       hl_counter->Fill(lon, lat, xl);  
492        
493        // This histograms is now filled with the number of entries.        // This histograms is now filled with the number of entries.
494        // Below we will divide with the time (in seconds) to get        // Below we will divide with the time (in seconds) to get
# Line 574  void Rate(TString *filename, TString out Line 559  void Rate(TString *filename, TString out
559    TH2F *hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");    TH2F *hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");
560    TH2F *hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");    TH2F *hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");
561    TH2F *hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");    TH2F *hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");
   TH2F *hb0_norm = (TH2F*) hb0_counter->Clone("hb0_norm");  
   TH2F *hl_norm = (TH2F*) hl_counter->Clone("hl_norm");  
562    
563    // Now we divide each histogram _counter with the time histogram    // Now we divide each histogram _counter with the time histogram
564    // obtBinTime to have an histogram _rate.  Note that, when a second    // obtBinTime to have an histogram _rate.  Note that, when a second
# Line 591  void Rate(TString *filename, TString out Line 574  void Rate(TString *filename, TString out
574    trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, "");    trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, "");
575    oss.str("");    oss.str("");
576    oss << basename.Data() << "_orbit_trigS111A.png";    oss << basename.Data() << "_orbit_trigS111A.png";
577    trigS111A_rate->SetMinimum(10);    trigS111A_rate->SetMinimum(9);
578    printHist(trigS111A_rate, mapFile, outDirectory, oss.str().c_str(), "S111A (Hz)", -width, height, true, 0);    printHist(trigS111A_rate, mapFile, outDirectory, oss.str().c_str(), "S111A (Hz)", -width, height, true, 0);
579    
580    antiCAS4_rate->Divide(antiCAS4_counter, obtBinTime, 1, 1, "");    antiCAS4_rate->Divide(antiCAS4_counter, obtBinTime, 1, 1, "");
581    oss.str("");    oss.str("");
582    oss << basename.Data() << "_orbit_CAS4.png";    oss << basename.Data() << "_orbit_CAS4.png";
583    antiCAS4_rate->SetMinimum(100);    antiCAS4_rate->SetMinimum(99);
584    printHist(antiCAS4_rate, mapFile, outDirectory, oss.str().c_str(), "CAS4 (Hz)", -width, height, true, 0);    printHist(antiCAS4_rate, mapFile, outDirectory, oss.str().c_str(), "CAS4 (Hz)", -width, height, true, 0);
585    
586    antiCAS3_rate->Divide(antiCAS3_counter, obtBinTime, 1, 1, "");    antiCAS3_rate->Divide(antiCAS3_counter, obtBinTime, 1, 1, "");
587    oss.str("");    oss.str("");
588    oss << basename.Data() << "_orbit_CAS3.png";    oss << basename.Data() << "_orbit_CAS3.png";
589    antiCAS3_rate->SetMinimum(100);    antiCAS3_rate->SetMinimum(99);
590    printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0);    printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0);
591    
592    event_rate->Divide(event_counter, obtBinTime, 1, 1, "");    event_rate->Divide(event_counter, obtBinTime, 1, 1, "");
# Line 614  void Rate(TString *filename, TString out Line 597  void Rate(TString *filename, TString out
597    trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, "");    trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, "");
598    oss.str("");    oss.str("");
599    oss << basename.Data() << "_orbit_trigS11andS12.png";    oss << basename.Data() << "_orbit_trigS11andS12.png";
600    antiCAS3_rate->SetMinimum(100);    trigS11andS12_rate->SetMinimum(99);
601    printHist(trigS11andS12_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12) (Hz)", -width, height, 1, 0);    printHist(trigS11andS12_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12) (Hz)", -width, height, 1, 0);
602    
603    trigS12andS21andS22_rate->Divide(trigS12andS21andS22_counter, obtBinTime, 1, 1, "");    trigS12andS21andS22_rate->Divide(trigS12andS21andS22_counter, obtBinTime, 1, 1, "");
604    oss.str("");    oss.str("");
605    oss << basename.Data() << "_orbit_trigS12andS21andS22.png";    oss << basename.Data() << "_orbit_trigS12andS21andS22.png";
606    antiCAS3_rate->SetMinimum(10);    trigS12andS21andS22_rate->SetMinimum(9);
607    printHist(trigS12andS21andS22_rate, mapFile, outDirectory, oss.str().c_str(), "(S12*S12*S21) (Hz)", -width, height, true, 0);    printHist(trigS12andS21andS22_rate, mapFile, outDirectory, oss.str().c_str(), "(S12*S12*S21) (Hz)", -width, height, true, 0);
608    
609    trigAndOr_rate->Divide(trigAndOr_counter, obtBinTime, 1, 1, "");    trigAndOr_rate->Divide(trigAndOr_counter, obtBinTime, 1, 1, "");
# Line 662  void Rate(TString *filename, TString out Line 645  void Rate(TString *filename, TString out
645    oss << basename.Data() << "_orbit_Beast.png";    oss << basename.Data() << "_orbit_Beast.png";
646    printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);    printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);
647    
   hb0_norm->Divide(hb0_counter, event_counter, 1, 1, "");  
   oss.str("");  
   oss << basename.Data() << "_orbit_B0.png";  
   printHist(hb0_norm, mapFile, outDirectory, oss.str().c_str(), "B_0 (G)", -width, height, 0, 0);  
   
   hl_norm->Divide(hl_counter, event_counter, 1, 1, "");  
   oss.str("");  
   oss << basename.Data() << "_orbit_L.png";  
   printHist(hl_norm, mapFile, outDirectory, oss.str().c_str(), "L shell", -width, height, 0, 0);  
   
648    
649    delete obtBinTime;    delete obtBinTime;
650    delete event_counter;    delete event_counter;
# Line 699  void Rate(TString *filename, TString out Line 672  void Rate(TString *filename, TString out
672    delete hbnorth_counter;    delete hbnorth_counter;
673    delete hbdown_counter;    delete hbdown_counter;
674    delete hbeast_counter;    delete hbeast_counter;
   delete hb0_counter;  
   delete hl_counter;  
675    delete hbabs_norm;    delete hbabs_norm;
676    delete hbnorth_norm;    delete hbnorth_norm;
677    delete hbdown_norm;    delete hbdown_norm;
678    delete hbeast_norm;    delete hbeast_norm;
   delete hb0_norm;  
   delete hl_norm;  
679    
680    rootFile->Close();    rootFile->Close();
681  }  }
# Line 737  int printHist(TH2F *h, TString mapFile, Line 706  int printHist(TH2F *h, TString mapFile,
706    // Create a canvas and draw the TH2F with a nice colormap for z    // Create a canvas and draw the TH2F with a nice colormap for z
707    // values, using log scale for z values, if requested, and setting    // values, using log scale for z values, if requested, and setting
708    // some title.    // some title.
709    TCanvas *canvas = new TCanvas("h", "passed histogram", width*2, height*2);    TCanvas *canvas = new TCanvas("h", "h histogram", width*2, height*2);
710    
711    if(use_log) {    if(use_log) canvas->SetLogz();
     canvas->SetLogz();  
   }  
712    
713    h->SetTitle(title);    h->SetTitle(title);
714    h->SetXTitle("Longitude (deg)");    h->SetXTitle("Longitude (deg)");
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23