/[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.3 by pam-rm2, Wed Dec 6 16:25:52 2006 UTC revision 1.5 by pam-rm2, Sat Dec 23 21:34:58 2006 UTC
# Line 210  void InitStyle() { Line 210  void InitStyle() {
210  void Rate(TString *filename, TString outDirectory = "", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000)  void Rate(TString *filename, TString outDirectory = "", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000)
211  {  {
212    // **** Offset to temporarily correct the TDatime bug ****/    // **** Offset to temporarily correct the TDatime bug ****/
213    offTime += 10000;    //  offTime += 10000;
214    //********************************************************/    //********************************************************/
215    
216    TTree                  *tr         = 0;    TTree                  *tr         = 0;
# 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?
# 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 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  }  }
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23