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

Diff of /quicklook/dataToXML/OrbitalRate.cpp

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

revision 1.2 by kusanagi, Tue May 2 14:32:23 2006 UTC revision 1.6 by kusanagi, Thu Jul 6 07:31:57 2006 UTC
# Line 5  Line 5 
5  *  *
6  * Description: .  * Description: .
7  *  *
 * Parameters:  
 *       base - the path where to find the PAMELA unpacked root file.  
 *       outDirectory - the path where to put the output file.  
 *       orbPath - the path where to find an orb format for the output.  
8  *  *
9  * version 1.0 - 27 April 2006  * 27 April 2006
10  * First implementation  * First implementation
11    *
12    * 6 July 2006
13    * Released version adding as an "emergency" Trigger and Anticounter plot. A future release should
14    * foreseen a general class called OrbitalRate, to be put into the "Utility" for plotting whatever
15    * data related to the orbit. Under a preliminary analysis this Class will contain
16    *
17    * OrbitalRate(String TLEFile) - file containing one or more TLE
18    *
19    * void    populateHisto(double data, double timeFromTLE) - fill the histogram
20    * TCanvas getCanvas() - return the populated TCanvas
21    *
22  */  */
23    #include <physics/anticounter/AnticounterEvent.h>
24    #include <physics/trigger/TriggerEvent.h>
25  #include <mcmd/McmdEvent.h>  #include <mcmd/McmdEvent.h>
26  #include <mcmd/McmdRecord.h>  #include <mcmd/McmdRecord.h>
27  #include <EventHeader.h>  #include <EventHeader.h>
# Line 22  Line 30 
30  #include "cTle.h"  #include "cTle.h"
31  #include "cEci.h"  #include "cEci.h"
32  #include "cOrbit.h"  #include "cOrbit.h"
33    #include "cJulian.h"
34  #include "TH2F.h"  #include "TH2F.h"
35  #include "TFrame.h"  #include "TFrame.h"
36  #include "TGraph.h"  #include "TGraph.h"
37  #include "TCanvas.h"  #include "TCanvas.h"
38  #include "TASImage.h"  #include "TASImage.h"
39  #include "TMarker.h"  #include "TMarker.h"
40    #include <TDatime.h>
41    
42  #include "TString.h"  #include "TString.h"
43  #include "TObjString.h"  #include "TObjString.h"
# Line 35  Line 45 
45  #include "TPaletteAxis.h"  #include "TPaletteAxis.h"
46  #include "TROOT.h"  #include "TROOT.h"
47  #include <sys/stat.h>  #include <sys/stat.h>
48    #include <fstream>
49    
50  #define TRUE  1  #define TRUE  1
51  #define FALSE 0  #define FALSE 0
# Line 104  void InitStyle() { Line 115  void InitStyle() {
115    
116  }  }
117    
118  void Rate(TString base, TString outDirectory = "", TString format = "jpg", TString mapFile = ""){  void Rate(TString base, TString outDirectory = "", TString format = "png", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000){
119          TTree                  *tr         = 0;          TTree                  *tr         = 0;
120          pamela::McmdEvent      *mcmdev     = 0;          pamela::McmdEvent      *mcmdev     = 0;
121          pamela::McmdRecord     *mcmdrc     = 0;          pamela::McmdRecord     *mcmdrc     = 0;
122          TArrayC                *mcmddata   = 0;          pamela::EventHeader    *eh         = 0;
123            pamela::PscuHeader     *ph         = 0;
124            TArrayC                *mcmddata;
125          ULong64_t               nevents    = 0;          ULong64_t               nevents    = 0;
126          ULong64_t               timesync   = 0;          double                  timesync   = 0;
         pamela::EventHeader    *eh         = 0;  
         pamela::PscuHeader     *ph         = 0;  
127          stringstream    oss;          stringstream    oss;
128          double offsetTime = 0;          double offsetTime = 0;
129          double absTime;          double absTime = 0;
130          UInt_t i = 0;          UInt_t i = 0;
131          UInt_t j = 0;          UInt_t j = 0;
132          struct stat buf;          struct stat buf;
# Line 129  void Rate(TString base, TString outDirec Line 140  void Rate(TString base, TString outDirec
140                    
141          TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();          TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
142          filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();          filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
143          // Test SGP4  
144          string str1 = "SGP4 Test";          string str1 = "RESURS-DK 1";
145          string str2 = "1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131";          string str2 = "1 29228U 06021A   06170.19643714  .00009962  00000-0  21000-3 0   196";
146          string str3 = "2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143";          string str3 = "2 29228 069.9363 054.7893 0167576 127.4359 017.0674 15.31839265   604";
147            if (tleFile != ""){
148                    fstream fileTle(tleFile.Data(),ios::in);
149                    if (fileTle.is_open()) {
150                            getline (fileTle,str1);
151                            getline (fileTle,str2);
152                            getline (fileTle,str3);
153                    }
154                    fileTle.close();
155            }
156    
157          cTle tle1(str1, str2, str3);          cTle tle1(str1, str2, str3);
158          cOrbit       orbit(tle1);          cOrbit       orbit(tle1);
# Line 155  void Rate(TString base, TString outDirec Line 175  void Rate(TString base, TString outDirec
175          nevents = tr->GetEntries();          nevents = tr->GetEntries();
176                    
177          bool timeFound = FALSE;          bool timeFound = FALSE;
178          while ((i < nevents) | !timeFound){          while (i < nevents) {
179                  tr->GetEntry(i);                  tr->GetEntry(i);
180                  recEntries = mcmdev->Records->GetEntries();                  recEntries = mcmdev->Records->GetEntries();
181                  while ((j < recEntries) | !timeFound){                  while (j < recEntries){
182                          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);                          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
183                          if (mcmdrc->ID1 == 0xE0){                          //mcmddata = mcmdrc->McmdData;
184                          mcmddata = mcmdrc->McmdData;                          //printf(" timesync TimeSync %i \n", (unsigned int)mcmddata->At(0));
185                          timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) +                          //It is a TimeSync?
186                                          (((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) +                          if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)){
187                                          (((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00)  +                                  mcmddata = mcmdrc->McmdData;
188                                          (((ULong64_t)mcmddata->At(3))&0x000000FF);                                  timesync = (((unsigned int)mcmddata->At(0)<<24)&0xFF000000) + (((unsigned int)mcmddata->At(1)<<16)&0x00FF0000) + (((unsigned int)mcmddata->At(2)<<8)&0x0000FF00) + (((unsigned int)mcmddata->At(3))&0x000000FF);
189                          offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);                                  timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
190                          timeFound = TRUE;                                  timeFound = TRUE;
191                                    //printf(" timesync TimeSync %i \n", timesync);
192                            }
193                            
194                            //It is an Inclination Mcmd?
195                            if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE2)){
196                                    mcmddata = mcmdrc->McmdData;
197                                    timesync = (((mcmddata->At(0) << 24) & 0xFF000000) + ((mcmddata->At(1) << 16) & 0x00FF0000) + ((mcmddata->At(2) << 8) & 0x0000FF00) + (mcmddata->At(3) & 0x000000FF))/128.0;
198                                    timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
199                                    timeFound = TRUE;
200                                    //printf(" timesync Inclination %16.8f \n", timesync);
201                          }                          }
202                  j++;                          
203                            if (timeFound) break;
204                            j++;
205                  }                  }
206            if (timeFound) break;
207          i++;          i++;
208          }          }
209                    
210            if (!timeFound) {
211                    printf("No timesync info have been found in the file %s", base.Data());
212                    exit(0);
213            }
214            
215            //Get the Julian date of the TLE Epoch
216            cJulian offTLE = cJulian(((int)tle1.getField(cTle::FLD_EPOCHYEAR) + 2000), tle1.getField(cTle::FLD_EPOCHDAY));
217            
218            //Get the Julian date of the Resours offset
219            TDatime offRes = TDatime(offDate, offTime);
220            cJulian offResours = cJulian(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond());
221            //Add to the Resours Offset the seconds differece beetwen the timesynch from MCMD and the relative OBT
222            offResours.addSec(timesync);
223            
224            //Get the MINUTES past since the TLE Epoch
225            offsetTime = (offResours.getDate() - offTLE.getDate()) * 24.0 * 60.0;
226          tr = (TTree*)rootFile->Get("Physics");          tr = (TTree*)rootFile->Get("Physics");
227          TBranch *headBr = tr->GetBranch("Header");          TBranch *headBr = tr->GetBranch("Header");
228          tr->SetBranchAddress("Header", &eh);          tr->SetBranchAddress("Header", &eh);
229            
230            /********** Anticounter **************/
231            pamela::anticounter::AnticounterEvent *antiev   = 0;
232            TH2F *antiRate40    = new TH2F("rate40", base, 360, -180, 180, 161, -80.5, 80.5);
233            TH2F *antiRate30    = new TH2F("rate30", base, 360, -180, 180, 161, -80.5, 80.5);
234            tr->SetBranchAddress("Anticounter", &antiev);
235            Int_t oldAntiRate40 = 0;
236            Int_t newAntiRate40 = 0;
237            Int_t oldAntiRate30 = 0;
238            Int_t newAntiRate30 = 0;
239            Int_t rateAnti40 = 0;
240            Int_t rateAnti30 = 0;
241            /********** Anticounter **************/
242            
243            /********** Trigger **************/
244            pamela::trigger::TriggerEvent *trigger   = 0;
245            TH2F *trigAnd   = new TH2F("trigAnd", base, 360, -180, 180, 161, -80.5, 80.5); // (S11+S12)*(S21+S22)*(S31+S32)
246            TH2F *trigS1    = new TH2F("trigS1", base, 360, -180, 180, 161, -80.5, 80.5); //(S11+S12)
247            TH2F *trigS111A = new TH2F("trigS111A", base, 360, -180, 180, 161, -80.5, 80.5); //(S111A)
248            tr->SetBranchAddress("Trigger", &trigger);
249            /********** Trigger **************/
250            
251          nevents = tr->GetEntries();          nevents = tr->GetEntries();
252          //Fill variables from root-ple          //Fill variables from root-ple
253          for (i = 0; i < nevents; i++){          for (i = 0; i < nevents; i++){
254            //for (i = 0; i < 1000; i++){
255                  tr->GetEntry(i);                  tr->GetEntry(i);
256                  ph = eh->GetPscuHeader();                  ph = eh->GetPscuHeader();
257                  absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60;                  absTime = offsetTime + (ph->GetOrbitalTime()*(1./60000.));
258                  orbit.getPosition(absTime, &eci);                  orbit.getPosition(absTime, &eci);
259                  coo = eci.toGeo();                  coo = eci.toGeo();
260                  rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat));                  /********** Anticounter **************/
261                  /*                  oldAntiRate40 = newAntiRate40;
262                  printf("             %16.8f %16.8f %16.8f\n",                  oldAntiRate30 = newAntiRate30;
263                  rad2deg(coo.m_Lat),                  newAntiRate40 = antiev->counters[0][6];
264                  rad2deg(coo.m_Lon),                  newAntiRate30 = antiev->counters[0][10];
265                  coo.m_Alt);                      if (newAntiRate40 < oldAntiRate40){
266                  */                          newAntiRate40 += 0xFFFF;
267                    }
268                    if (newAntiRate30 < oldAntiRate30){
269                            newAntiRate30 += 0xFFFF;
270                    }
271                    rateAnti40 = newAntiRate40 - oldAntiRate40;
272                    rateAnti30 = newAntiRate30 - oldAntiRate30;
273                    /********** Anticounter **************/
274                    
275                    if (coo.m_Lon > PI) {
276                            rate->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat));
277                            antiRate30->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), rateAnti30);
278                            antiRate40->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), rateAnti40);
279                            trigAnd->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), (1/4.)*trigger->trigrate[1]);
280                            trigS1->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), 1.*trigger->pmtpl[0]);
281                            trigS111A->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), 1.*trigger->pmtcount1[0]);
282                    } else {
283                            rate->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat));
284                            antiRate30->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), rateAnti30);
285                            antiRate40->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), rateAnti40);
286                            trigAnd->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), (1/4.)*trigger->trigrate[1]);
287                            trigS1->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), 1.*trigger->pmtpl[0]);
288                            trigS111A->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), 1.*trigger->pmtcount1[0]);
289                            //printf(" Time: %16.8f Lat: %16.8f  Long: %16.8f  Alt: %16.8f\n", absTime, rad2deg(coo.m_Lat), rad2deg(coo.m_Lon), coo.m_Alt);
290                    }
291          }          }
292          double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;          double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;
293    
         int ptcnt=0, color=0;  
294          TMarker *tma=NULL;          TMarker *tma=NULL;
         TLine *tli=NULL;  
         double step=0;  
   
   
295          TImage *tImage=TImage::Open(mapFile);          TImage *tImage=TImage::Open(mapFile);
296          int width=(int)(tImage->GetWidth()*0.80);          int width=(int)(tImage->GetWidth()*0.80);
297          int height=(int)(tImage->GetHeight()*0.80);          int height=(int)(tImage->GetHeight()*0.80);
298          InitStyle();          InitStyle();
299          TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars          TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars
300                            
301          TH1F *hframe=NULL;          TH1F *hframe=NULL;
302          hframe=gPad->DrawFrame(-180,-90,180,90);          hframe=gPad->DrawFrame(-180,-90,180,90);
303            
304          oss.str("");          oss.str("");
305          oss << filename << " - Event Rate (Hz)";          oss << filename << " - Event Rate (Hz)";
306          hframe->SetTitle(oss.str().c_str());          hframe->SetTitle(oss.str().c_str());
307          hframe->SetXTitle("Longitude (deg)");          hframe->SetXTitle("Longitude (deg)");
308          hframe->SetYTitle("Latitude (deg)");          hframe->SetYTitle("Latitude (deg)");
         c1->cd();  
309                    
310            /********** Anticounter **************/
311            TCanvas *anti30 = new TCanvas("anti30","rate/orbit",-width, height); // - : removes the menu bars
312            TH1F *hframe30=NULL;
313            hframe30=gPad->DrawFrame(-180,-90,180,90);
314            
315            oss.str("");
316            oss << filename << " - CAS3 rate/degree";
317            hframe30->SetTitle(oss.str().c_str());
318            hframe30->SetXTitle("Longitude (deg)");
319            hframe30->SetYTitle("Latitude (deg)");
320            
321            TCanvas *anti40 = new TCanvas("anti40","rate/orbit",-width, height); // - : removes the menu bars
322            TH1F *hframe40=NULL;
323            hframe40=gPad->DrawFrame(-180,-90,180,90);
324            
325            oss.str("");
326            oss << filename << " - CAS4 rate/degree";
327            hframe40->SetTitle(oss.str().c_str());
328            hframe40->SetXTitle("Longitude (deg)");
329            hframe40->SetYTitle("Latitude (deg)");
330            /********** Anticounter **************/
331                    
332            /********** Trigger **************/
333            TCanvas *trigAndC = new TCanvas("AND","rate/orbit",-width, height); // - : removes the menu bars
334            TH1F *trigAndF=NULL;
335            trigAndF=gPad->DrawFrame(-180,-90,180,90);
336            
337            oss.str("");
338            oss << filename << " - (S11+S12)*(S21+S22)*(S31+S32) rate/degree";
339            trigAndF->SetTitle(oss.str().c_str());
340            trigAndF->SetXTitle("Longitude (deg)");
341            trigAndF->SetYTitle("Latitude (deg)");
342            
343            TCanvas *trigS1C = new TCanvas("trigS1","rate/orbit",-width, height); // - : removes the menu bars
344            TH1F *trigS1F=NULL;
345            trigS1F=gPad->DrawFrame(-180,-90,180,90);
346            
347            oss.str("");
348            oss << filename << " - (S11+S12) rate/degree";
349            trigS1F->SetTitle(oss.str().c_str());
350            trigS1F->SetXTitle("Longitude (deg)");
351            trigS1F->SetYTitle("Latitude (deg)");
352                    
353            TCanvas *trigS111AC = new TCanvas("trigS111A","rate/orbit",-width, height); // - : removes the menu bars
354            TH1F *trigS111AF=NULL;
355            trigS111AF=gPad->DrawFrame(-180,-90,180,90);
356                    
357            oss.str("");
358            oss << filename << " - S111A rate/degree";
359            trigS111AF->SetTitle(oss.str().c_str());
360            trigS111AF->SetXTitle("Longitude (deg)");
361            trigS111AF->SetYTitle("Latitude (deg)");
362            /********** Trigger **************/
363                    
364            c1->cd();
365          TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);          TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);
366          p2->Draw();  
367          p2->cd();          p2->cd();
368          TPaletteAxis *tpa=NULL;          TPaletteAxis *tpa=NULL;
369          TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);          TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);
370          forpal->SetAxisColor(kBlack); //Delete the stat box          forpal->SetAxisColor(kBlack); //Delete the stat box
371          forpal->SetStats(0); //Delete the stat box          forpal->SetStats(0); //Delete the stat box
372          forpal->SetMinimum(0.1);          forpal->SetMinimum(0.1);
373          forpal->SetMaximum(200);          forpal->SetMaximum(rate->GetMaximum());
374            p2->SetLogy();
375          forpal->SetBinContent(5,1);    // just to initialize the histo          forpal->SetBinContent(5,1);    // just to initialize the histo
376          forpal->SetContour(50);          forpal->SetContour(50);
377          TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));          TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
         step=forpal->GetContourLevel(1)-forpal->GetContourLevel(0);  
378          tpa=new TPaletteAxis(184,-90,195,90,forpal);          tpa=new TPaletteAxis(184,-90,195,90,forpal);
379          tpa->SetLabelColor(kWhite);          tpa->SetLabelColor(kWhite);
380          forpal->Draw("colz");          forpal->Draw("colz");
381            
382            c1->cd();
383            p2->DrawClone();
384            
385            /********** Anticounter **************/
386            Double_t stepPP30   = (antiRate30->GetMaximum() / 50.0);
387            forpal->SetMaximum(antiRate30->GetMaximum());
388            anti30->cd();
389            p2->DrawClone();
390            
391            Double_t stepPP40   = (antiRate40->GetMaximum() / 50.0);        
392            forpal->SetMaximum(antiRate40->GetMaximum());
393            anti40->cd();
394            p2->DrawClone();
395            /********** Anticounter **************/
396            
397            /********** Trigger **************/
398            Double_t trigAndPP   = (trigAnd->GetMaximum() / 50.0);
399            forpal->SetMaximum(trigAnd->GetMaximum());
400            trigAndC->cd();
401            p2->DrawClone();
402            
403            Double_t trigS1PP   = (trigS1->GetMaximum() / 50.0);
404            forpal->SetMaximum(trigS1->GetMaximum());
405            trigS1C->cd();
406            p2->DrawClone();
407            
408            Double_t trigS111APP   = (trigS111A->GetMaximum() / 50.0);
409            forpal->SetMaximum(trigS111A->GetMaximum());
410            trigS111AC->cd();
411            p2->DrawClone();
412            /********** Trigger **************/
413    
414          c1->cd();          c1->cd();
415          TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);          TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);
416          p1->Draw();          p1->Draw();
417            
418          p1->cd();          p1->cd();
419          tImage->Draw("X");          tImage->Draw("X");
420            
421            /********** Anticounter **************/
422            anti30->cd();
423            p1->DrawClone();
424            anti40->cd();
425            p1->DrawClone();
426            /********** Anticounter **************/
427            
428            /********** Trigger **************/
429            trigAndC->cd();
430            p1->DrawClone();
431            trigS1C->cd();
432            p1->DrawClone();
433            trigS111AC->cd();
434            p1->DrawClone();
435            /********** Trigger **************/
436            
437          c1->cd();          c1->cd();
438          gPad->RedrawAxis();          gPad->RedrawAxis();
439          c1->Update();          c1->Update();
440                    
441            /********** Anticounter **************/
442            anti30->cd();
443            gPad->RedrawAxis();
444            anti30->Update();
445            anti40->cd();
446            gPad->RedrawAxis();
447            anti40->Update();
448            /********** Anticounter **************/
449            
450            /********** Trigger **************/
451            trigAndC->cd();
452            gPad->RedrawAxis();
453            trigAndC->Update();
454            
455            trigS1C->cd();
456            gPad->RedrawAxis();
457            trigS1C->Update();
458            
459            trigS111AC->cd();
460            gPad->RedrawAxis();
461            trigS111AC->Update();
462            /********** Trigger **************/
463            
464    c1->cd();    c1->cd();
   ptcnt=0;  
465    tma=new TMarker();    tma=new TMarker();
466    tma->SetMarkerStyle(4);    tma->SetMarkerStyle(21);
467    tli=new TLine();    tma->SetMarkerSize(0.5);
468    tli->SetLineColor(kMagenta);    Double_t color = 0.0;  
469    Stat_t freq;    Double_t freq  = 0.0;
470        for (Int_t kk = 0; kk < 360; kk++){    Double_t stepPP = (rate->GetMaximum() / 50.0);
471          for (Int_t jj = 0; jj < 161; jj++){    //Double_t cutoffFreq = stepPP/1000.0;
472      //Double_t alpha = log10(cutoffFreq/stepPP) * (-1.0); //Obviously is 3.0
473      
474            /********** Anticounter **************/
475            Double_t color30 = 0.0;  
476            Double_t freq30  = 0.0;
477            //Double_t cutoffFreq30 = stepPP30/1000.0;
478            //Double_t alpha30 = log10(cutoffFreq30/stepPP30) * (-1.0); //Obviously is 3.0
479            
480            Double_t color40 = 0.0;  
481            Double_t freq40  = 0.0;
482            //Double_t cutoffFreq40 = stepPP40/1000.0;
483            //Double_t alpha40 = log10(cutoffFreq40/stepPP40) * (-1.0); //Obviously is 3.0
484            /********** Anticounter **************/
485      
486            /********** Trigger **************/
487            Double_t trigAndCol = 0.0;  
488            Double_t trigAndFreq  = 0.0;
489            
490            Double_t trigS1Col = 0.0;  
491            Double_t trigS1Freq  = 0.0;
492            
493            Double_t trigS111ACol = 0.0;  
494            Double_t trigS111AFreq  = 0.0;
495            /********** Trigger **************/
496      
497      for (Int_t kk = 0; kk < 360; kk++){
498            for (Int_t jj = 0; jj < 161; jj++){
499                  freq = rate->GetBinContent(kk, jj);                  freq = rate->GetBinContent(kk, jj);
500                  if (freq>0) {                  
501                          posx=(kk - 180); posy= jj - 80.5;                  /********** Anticounter **************/
502                    freq30 = antiRate30->GetBinContent(kk, jj);
503                    freq40 = antiRate40->GetBinContent(kk, jj);
504                    /********** Anticounter **************/
505                    
506                    /********** Trigger **************/
507                    trigAndFreq = trigAnd->GetBinContent(kk, jj);
508                    trigS1Freq = trigS1->GetBinContent(kk, jj);
509                    trigS111AFreq = trigS111A->GetBinContent(kk, jj);
510                    /********** Trigger **************/
511                    
512                    //printf(" Freq:  %16.8f \n", freq);
513                    posx=kk - 180;
514                    posy= jj - 80.5;
515                    if (freq > 0) {
516                            c1->cd();
517                          // color: palette colors from 51 to 100 ie 50 levels                          // color: palette colors from 51 to 100 ie 50 levels
518                          color=51+(int) ((log10((rate==0)?0.1:freq)-log10(0.1))/step); // step defined by palette                          //color = ((log10(freq/stepPP))/alpha)*49.0 + 100.0; // step defined by palette
519                          if (color>100) color=100; // just in case if max rate is not max...                          color = 50.0 + (freq / stepPP); // step defined by palette
520                          tma->SetMarkerColor(color);                          //printf(" Color: %16.8f \n", color);
521                          if (!(posx<0 && oldposx>0) && oldposx!=-1000 && oldposy!=-1000) tli->DrawLine(oldposx,oldposy,posx,posy);                          //if (freq < stepPP) color = 51;
522                            //if (color < -100.0) color = 100.0; // just in case if max rate is not max...
523                            tma->SetMarkerColor((int)color);
524                          tma->DrawMarker(posx,posy);                          tma->DrawMarker(posx,posy);
                         oldposx=posx;  
                         oldposy=posy;    
525                  }                  }
526          }                  
527        }                  /********** Anticounter **************/
528                    if (freq30 > 0) {
529                            anti30->cd();
530                            color30 = 50.0 + (freq30 / stepPP30); // step defined by palette
531                            //printf(" Color: %16.8f \n", color30);
532                            //color30 = ((log10(freq30/stepPP30))/alpha30)*49.0 + 100.0; // step defined by palette
533                            //printf(" Color: %16.8f Step: %16.8f Freq:  %16.8f \n", color, stepPP, freq);
534                            //if (freq30 < cutoffFreq30) color30 = 51;
535                            //if (color30 > 100.0) color30 = 100.0; // just in case if max rate is not max...
536                            tma->SetMarkerColor((int)color30);
537                            tma->DrawMarker(posx,posy);
538                    }
539                    
540                    if (freq40 > 0) {
541                            anti40->cd();
542                            color40 = 50.0 + (freq40 / stepPP40); // step defined by palette
543                            //color40 = ((log10(freq40/stepPP40))/alpha40)*49.0 + 100.0; // step defined by palette
544                            //printf(" Color: %16.8f Step: %16.8f Freq:  %16.8f \n", color, stepPP, freq);
545                            //if (freq40 < cutoffFreq40) color40 = 51;
546                            //if (color40 > 100.0) color40 = 100.0; // just in case if max rate is not max...
547                            tma->SetMarkerColor((int)color40);
548                            tma->DrawMarker(posx,posy);
549                    }
550                    /********** Anticounter **************/
551                    
552                    /********** Trigger **************/
553                    if (trigAndFreq > 0) {
554                            trigAndC->cd();
555                            trigAndCol = 50.0 + (trigAndFreq / trigAndPP); // step defined by palette
556                            tma->SetMarkerColor((int)trigAndCol);
557                            tma->DrawMarker(posx,posy);
558                    }
559                    
560                    if (trigS1Freq > 0) {
561                            trigS1C->cd();
562                            trigS1Col = 50.0 + (trigS1Freq / trigS1PP); // step defined by palette
563                            tma->SetMarkerColor((int)trigS1Col);
564                            tma->DrawMarker(posx,posy);
565                    }
566                    
567                    if (trigS111AFreq > 0) {
568                            trigS111AC->cd();
569                            trigS111ACol = 50.0 + (trigS111AFreq / trigS111APP); // step defined by palette
570                            tma->SetMarkerColor((int)trigS111ACol);
571                            tma->DrawMarker(posx,posy);
572                    }
573                    /********** Anticounter **************/
574            }
575     }
576    
577          oss.str("");          oss.str("");
578          if (outDirectory == "") {          if (outDirectory == "") {
579                  oss << filename.Data() << "_OrbitRate." << format.Data();                  oss << filename.Data() << "_OrbitRate." << format.Data();
580          } else {          } else {
581                  oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();                  oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();
582          }          }
583          c1->SaveAs(oss.str().c_str());          c1->SaveAs(oss.str().c_str());
584            
585            oss.str("");
586            if (outDirectory == "") {
587                    oss << filename.Data() << "_CAS3." << format.Data();
588            } else {
589                    oss << outDirectory.Data() << filename.Data() << "_CAS3." << format.Data();
590            }
591            anti30->SaveAs(oss.str().c_str());
592            
593            oss.str("");
594            if (outDirectory == "") {
595                    oss << filename.Data() << "_CAS4." << format.Data();
596            } else {
597                    oss << outDirectory.Data() << filename.Data() << "_CAS4." << format.Data();
598            }
599            anti40->SaveAs(oss.str().c_str());
600            
601            oss.str("");
602            if (outDirectory == "") {
603                    oss << filename.Data() << "_TrigAND." << format.Data();
604            } else {
605                    oss << outDirectory.Data() << filename.Data() << "_TrigAND." << format.Data();
606            }
607            trigAndC->SaveAs(oss.str().c_str());
608            
609            oss.str("");
610            if (outDirectory == "") {
611                    oss << filename.Data() << "_TrigS1." << format.Data();
612            } else {
613                    oss << outDirectory.Data() << filename.Data() << "_TrigS1." << format.Data();
614            }
615            trigS1C->SaveAs(oss.str().c_str());
616                    
617            oss.str("");
618            if (outDirectory == "") {
619                    oss << filename.Data() << "_trigS111A." << format.Data();
620            } else {
621                    oss << outDirectory.Data() << filename.Data() << "_trigS111A." << format.Data();
622            }
623            trigS111AC->SaveAs(oss.str().c_str());
624    
625          rootFile->Close();          rootFile->Close();
626  }  }
627    
628    
629    
630    
631  int main(int argc, char* argv[]){  int main(int argc, char* argv[]){
632      TString outDir  = "./";      TString outDir  = "./";
633      TString mapFile = "";      TString mapFile = "";
634      TString format  = "jpg";      TString tleFile = "";
635        TString format  = "png";
636        int offDate = 20060614;
637        int offTime = 210000;
638            
639   if (argc < 2){   if (argc < 2){
640      printf("You have to insert at least the file to analyze and the mapFile \n");      printf("You have to insert at least the file to analyze and the mapFile \n");
# Line 292  int main(int argc, char* argv[]){ Line 646  int main(int argc, char* argv[]){
646          printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");          printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
647          printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");          printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
648          printf( "\t --help                 Print this help and exit \n");          printf( "\t --help                 Print this help and exit \n");
649            printf( "\t -tleFile[path]         Path where to find the tle infos [default dummyOrbit] \n");
650            printf( "\t\tThe tle file have to satisfy a 3-line structure like (this is the included dummyOrbit)\n ");
651            printf( "\t\t\tGP4 Test\n");
652            printf( "\t\t\t1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131\n");
653            printf( "\t\t\t2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143\n");
654          printf( "\t -outDir[path]          Path where to put the output [default ~/tmp] \n");          printf( "\t -outDir[path]          Path where to put the output [default ~/tmp] \n");
655          printf( "\t -format[jpg|gif|ps]    Format for output files [default 'jpg'] \n");          printf( "\t -format[jpg|gif|ps]    Format for output files [default 'png'] \n");
656            printf( "\t -offDate               Date of resetting of the Resource counter [format YYMMDD default 20060614] \n");
657            printf( "\t -offTime               Time of resetting of the Resource counter [format HHMMSS default 210000] \n");
658          exit(1);          exit(1);
659    }    }
660    
# Line 309  int main(int argc, char* argv[]){ Line 670  int main(int argc, char* argv[]){
670          }          }
671      }      }
672    
673       if (!strcmp(argv[i], "-tle")){
674            if (++i >= argc){
675                printf( "-tle needs arguments. \n");
676                printf( "Try '--help' for more information. \n");
677                exit(1);
678            } else {
679                tleFile = argv[i];
680                continue;
681            }
682        }
683    
684        if (!strcmp(argv[i], "-offTime")){
685          if (++i >= argc){
686            printf( "-offTime needs arguments. \n");
687            printf( "Try '--help' for more information. \n");
688            exit(1);
689          }
690          else{
691            offTime = atol(argv[i]);
692            continue;
693          }
694        }
695    
696        if (!strcmp(argv[i], "-offDate")){
697          if (++i >= argc){
698            printf( "-offDate needs arguments. \n");
699            printf( "Try '--help' for more information. \n");
700            exit(1);
701          }
702          else{
703            offDate = atol(argv[i]);
704            continue;
705          }
706        }
707    
708      if (!strcmp(argv[i], "-map")){      if (!strcmp(argv[i], "-map")){
709          if (++i >= argc){          if (++i >= argc){
710              printf( "-map needs arguments. \n");              printf( "-map needs arguments. \n");
# Line 332  int main(int argc, char* argv[]){ Line 728  int main(int argc, char* argv[]){
728      }      }
729  }  }
730  if (mapFile != ""){  if (mapFile != ""){
731          Rate(argv[1], outDir, format, mapFile);          Rate(argv[1], outDir, format, mapFile, tleFile, offDate, offTime);
732  } else {  } else {
733          printf("You have to insert at least the file to analyze and the mapFile \n");          printf("You have to insert at least the file to analyze and the mapFile \n");
734          printf("Try '--help' for more information. \n");          printf("Try '--help' for more information. \n");
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23