/** * Rate * author Nagni * version 1.0 - 27 April 2006 * * Description: . * * 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. * * version 1.0 - 27 April 2006 * First implementation */ #include #include #include #include #include #include "cTle.h" #include "cEci.h" #include "cOrbit.h" #include "TH2F.h" #include "TFrame.h" #include "TGraph.h" #include "TCanvas.h" #include "TASImage.h" #include "TMarker.h" #include "TString.h" #include "TObjString.h" #include "TStyle.h" #include "TPaletteAxis.h" #include "TROOT.h" #define TRUE 1 #define FALSE 0 /** * * @param base * @param outDirectory * @param xslPath */ void InitStyle() { gROOT->SetStyle("Plain"); TStyle *myStyle[2], *tempo; myStyle[0]=new TStyle("StyleWhite", "white"); myStyle[1]=new TStyle("StyleBlack", "black"); tempo=gStyle; Int_t linecol, bkgndcol, histcol; for(Int_t style=0; style<2; style++) { linecol=kWhite*style+kBlack*(1-style); bkgndcol=kBlack*style+kWhite*(1-style); histcol=kYellow*style+kBlack*(1-style); // was 95 myStyle[style]->Copy(*tempo); myStyle[style]->SetCanvasBorderMode(0); myStyle[style]->SetCanvasBorderSize(1); myStyle[style]->SetFrameBorderSize(1); myStyle[style]->SetFrameBorderMode(0); myStyle[style]->SetPadBorderSize(1); myStyle[style]->SetStatBorderSize(1); myStyle[style]->SetTitleBorderSize(1); myStyle[style]->SetPadBorderMode(0); myStyle[style]->SetPalette(1,0); myStyle[style]->SetPaperSize(20,27); myStyle[style]->SetFuncColor(kRed); myStyle[style]->SetFuncWidth(1); myStyle[style]->SetLineScalePS(1); myStyle[style]->SetCanvasColor(bkgndcol); myStyle[style]->SetAxisColor(linecol,"XYZ"); myStyle[style]->SetFrameFillColor(bkgndcol); myStyle[style]->SetFrameLineColor(linecol); myStyle[style]->SetLabelColor(linecol,"XYZ"); myStyle[style]->SetPadColor(bkgndcol); myStyle[style]->SetStatColor(bkgndcol); myStyle[style]->SetStatTextColor(linecol); myStyle[style]->SetTitleColor(linecol,"XYZ"); myStyle[style]->SetTitleFillColor(bkgndcol); myStyle[style]->SetTitleTextColor(linecol); myStyle[style]->SetLineColor(linecol); myStyle[style]->SetMarkerColor(histcol); myStyle[style]->SetTextColor(linecol); myStyle[style]->SetGridColor((style)?13:kBlack); myStyle[style]->SetHistFillStyle(1001*(1-style)); myStyle[style]->SetHistLineColor(histcol); myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow); } myStyle[1]->cd(); gROOT->ForceStyle(); } void Rate(TString base, TString outDirectory = "", TString format = "jpg"){ TTree *tr = 0; pamela::McmdEvent *mcmdev = 0; pamela::McmdRecord *mcmdrc = 0; TArrayC *mcmddata = 0; ULong64_t nevents = 0; ULong64_t timesync = 0; pamela::EventHeader *eh = 0; pamela::PscuHeader *ph = 0; stringstream oss; double offsetTime = 0; double absTime; UInt_t i = 0; UInt_t j = 0; TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString(); filename = ((TObjString*)filename.Tokenize('.')->First())->GetString(); // Test SGP4 string str1 = "SGP4 Test"; string str2 = "1 25544U 98067A 06117.32388940 .00009459 00000-0 64427-4 0 8131"; string str3 = "2 25544 51.6388 89.2928 0009043 155.3293 354.6512 15.75673697425143"; cTle tle1(str1, str2, str3); cOrbit orbit(tle1); cEci eci; cCoordGeo coo; TH2F *rate = new TH2F("rate", base, 360, -180, 180, 161, -80.5, 80.5); TFile *rootFile = new TFile(base); if (rootFile->IsZombie()) { printf("The %s file does not exist", base.Data()); exit(0); } /* * process Mcmd */ long int recEntries; tr = (TTree*)rootFile->Get("Mcmd"); tr->SetBranchAddress("Mcmd", &mcmdev); nevents = tr->GetEntries(); bool timeFound = FALSE; while ((i < nevents) | !timeFound){ tr->GetEntry(i); recEntries = mcmdev->Records->GetEntries(); while ((j < recEntries) | !timeFound){ mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j); if (mcmdrc->ID1 == 0xE0){ mcmddata = mcmdrc->McmdData; timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) + (((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) + (((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00) + (((ULong64_t)mcmddata->At(3))&0x000000FF); offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.); timeFound = TRUE; } j++; } i++; } tr = (TTree*)rootFile->Get("Physics"); TBranch *headBr = tr->GetBranch("Header"); tr->SetBranchAddress("Header", &eh); nevents = tr->GetEntries(); //Fill variables from root-ple for (i = 0; i < nevents; i++){ tr->GetEntry(i); ph = eh->GetPscuHeader(); absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60; orbit.getPosition(absTime, &eci); coo = eci.toGeo(); rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat)); /* printf(" %16.8f %16.8f %16.8f\n", rad2deg(coo.m_Lat), rad2deg(coo.m_Lon), coo.m_Alt); */ } double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000; int ptcnt=0, color=0; TMarker *tma=NULL; TLine *tli=NULL; double step=0; TImage *tImage=TImage::Open("Data/terra4.png"); int width=(int)(tImage->GetWidth()*0.80); int height=(int)(tImage->GetHeight()*0.80); InitStyle(); TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars TH1F *hframe=NULL; hframe=gPad->DrawFrame(-180,-90,180,90); oss.str(""); oss << filename << " - Event Rate (Hz)"; hframe->SetTitle(oss.str().c_str()); hframe->SetXTitle("Longitude (deg)"); hframe->SetYTitle("Latitude (deg)"); c1->cd(); TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1); p2->Draw(); p2->cd(); TPaletteAxis *tpa=NULL; TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2); forpal->SetAxisColor(kBlack); //Delete the stat box forpal->SetStats(0); //Delete the stat box forpal->SetMinimum(0.1); forpal->SetMaximum(200); forpal->SetBinContent(5,1); // just to initialize the histo forpal->SetContour(50); TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette")); step=forpal->GetContourLevel(1)-forpal->GetContourLevel(0); tpa=new TPaletteAxis(184,-90,195,90,forpal); tpa->SetLabelColor(kWhite); forpal->Draw("colz"); c1->cd(); TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92); p1->Draw(); p1->cd(); tImage->Draw("X"); c1->cd(); gPad->RedrawAxis(); c1->Update(); c1->cd(); ptcnt=0; tma=new TMarker(); tma->SetMarkerStyle(4); tli=new TLine(); tli->SetLineColor(kMagenta); Stat_t freq; for (Int_t kk = 0; kk < 360; kk++){ for (Int_t jj = 0; jj < 161; jj++){ freq = rate->GetBinContent(kk, jj); if (freq>0) { posx=(kk - 180); posy= jj - 80.5; // color: palette colors from 51 to 100 ie 50 levels color=51+(int) ((log10((rate==0)?0.1:freq)-log10(0.1))/step); // step defined by palette if (color>100) color=100; // just in case if max rate is not max... tma->SetMarkerColor(color); if (!(posx<0 && oldposx>0) && oldposx!=-1000 && oldposy!=-1000) tli->DrawLine(oldposx,oldposy,posx,posy); tma->DrawMarker(posx,posy); oldposx=posx; oldposy=posy; } } } oss.str(""); if (outDirectory == "") { oss << filename.Data() << "_OrbitRate." << format.Data(); } else { oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data(); } c1->SaveAs(oss.str().c_str()); rootFile->Close(); } int main(int argc, char* argv[]){ TString outDir = "./"; TString format = "jpg"; if (argc < 2){ printf("You have to insert at least the file to analyze \n"); printf("Try '--help' for more information. \n"); exit(1); } if (!strcmp(argv[1], "--help")){ printf( "Usage: LogToXML FILE [OPTION] \n"); printf( "\t --help Print this help and exit \n"); printf( "\t -outDir[path] Path where to put the output [default ~/tmp] \n"); printf( "\t -format[jpg|gif|ps] Format for output files [default 'jpg'] \n"); exit(1); } for (int i = 2; i < argc; i++){ if (!strcmp(argv[i], "-outDir")){ if (++i >= argc){ printf( "-outDir needs arguments. \n"); printf( "Try '--help' for more information. \n"); exit(1); } else { outDir = argv[i]; continue; } } if (!strcmp(argv[i], "-format")) if (++i >= argc){ printf( "-format needs arguments. \n"); printf( "Try '--help' for more information. \n"); exit(1); } else { format = argv[i]; continue; } } Rate(argv[1], outDir, format); }