| 1 | #define PI 3.1415
|
|---|
| 2 | #define D2R PI/180.
|
|---|
| 3 | #define R2D 180./PI
|
|---|
| 4 |
|
|---|
| 5 | void createTables( string sInFile, string sOutFile )
|
|---|
| 6 | {
|
|---|
| 7 |
|
|---|
| 8 | gStyle->SetOptStat(0);
|
|---|
| 9 |
|
|---|
| 10 | // you can run without an output file by leaving the sOutFile string blank ""
|
|---|
| 11 | bool bOut = false;
|
|---|
| 12 |
|
|---|
| 13 | // Open output file
|
|---|
| 14 | if( sOutFile.c_str() != "" )
|
|---|
| 15 | {
|
|---|
| 16 | bOut = true;
|
|---|
| 17 | TFile *fOutFile = new TFile( sOutFile.c_str(), "CREATE" );
|
|---|
| 18 | }
|
|---|
| 19 |
|
|---|
| 20 | // open input file
|
|---|
| 21 | TFile *fInFile = new TFile( sInFile.c_str() );
|
|---|
| 22 |
|
|---|
| 23 | if( fInFile->IsZombie() )
|
|---|
| 24 | {
|
|---|
| 25 | cout << "File Error " << sInFile.c_str() << endl;
|
|---|
| 26 | exit(-1);
|
|---|
| 27 | }
|
|---|
| 28 |
|
|---|
| 29 | // Get the tree
|
|---|
| 30 | TTree *hTree = (TTree*)fInFile->Get("hillastree");
|
|---|
| 31 |
|
|---|
| 32 | // Variables needed
|
|---|
| 33 | float fMCEnergy;
|
|---|
| 34 | float fMCZd;
|
|---|
| 35 | double fLength;
|
|---|
| 36 | double fWidth;
|
|---|
| 37 | double fDistance;
|
|---|
| 38 | double fSize;
|
|---|
| 39 |
|
|---|
| 40 | // get the branches
|
|---|
| 41 | hTree->SetBranchAddress("MCEnergy",&fMCEnergy);
|
|---|
| 42 | hTree->SetBranchAddress("MCZd",&fMCZd);
|
|---|
| 43 | hTree->SetBranchAddress("Length",&fLength);
|
|---|
| 44 | hTree->SetBranchAddress("Width",&fWidth);
|
|---|
| 45 | hTree->SetBranchAddress("Size",&fSize);
|
|---|
| 46 | hTree->SetBranchAddress("Distance",&fDistance);
|
|---|
| 47 |
|
|---|
| 48 | // Setup the histograms
|
|---|
| 49 | TH1F *hMCEnergy = new TH1F("hMCEnergy","True Energy",100,-1.,2.);
|
|---|
| 50 | TH1F *hMCZd = new TH1F("hMCZd","MC Zenith",100,0.,40.);
|
|---|
| 51 | TH1F *hLength = new TH1F("hLength","Length",100,0.,150.);
|
|---|
| 52 | TH1F *hWidth = new TH1F("hWidth","Width",100,0.,50.);
|
|---|
| 53 | TH1F *hDistance = new TH1F("hDistance","Distance",100,0.,190.);
|
|---|
| 54 | TH1F *hSize = new TH1F("hSize","Size",100,0.,6.);
|
|---|
| 55 |
|
|---|
| 56 | TH2F *hTableEnergy = new TH2F("hTableEnergy","Energy",30,1.,4.,30,0.,170.);
|
|---|
| 57 | TH2F *hTableLength = new TH2F("hTableLength","Length",30,1.,4.,30,0.,170.);
|
|---|
| 58 | TH2F *hTableWidth = new TH2F("hTableWidth","Width",30,1.,4.,30,0.,170.);
|
|---|
| 59 |
|
|---|
| 60 | // This is the normalisation histogram
|
|---|
| 61 | TH2F *hTableNorm = new TH2F("hNTableNorm","Norm",30,1.,4.,30,0.,170.);
|
|---|
| 62 |
|
|---|
| 63 | // Titles
|
|---|
| 64 | hTableEnergy->GetXaxis()->SetTitle("log10(Size)");
|
|---|
| 65 | hTableLength->GetXaxis()->SetTitle("log10(Size)");
|
|---|
| 66 | hTableWidth->GetXaxis()->SetTitle("log10(Size)");
|
|---|
| 67 |
|
|---|
| 68 | hTableEnergy->GetYaxis()->SetTitle("Distance");
|
|---|
| 69 | hTableLength->GetYaxis()->SetTitle("Distance");
|
|---|
| 70 | hTableWidth->GetYaxis()->SetTitle("Distance");
|
|---|
| 71 |
|
|---|
| 72 | // Half the entries
|
|---|
| 73 | int BIG = int(hTree->GetEntries()/2.);
|
|---|
| 74 |
|
|---|
| 75 | // Loop over tree
|
|---|
| 76 | for( int i = 0; i < BIG; i++ )
|
|---|
| 77 | {
|
|---|
| 78 |
|
|---|
| 79 |
|
|---|
| 80 | hTree->GetEntry(i);
|
|---|
| 81 | hMCEnergy->Fill(log10(fMCEnergy/1000.));
|
|---|
| 82 | hSize->Fill(log10(fSize));
|
|---|
| 83 | hDistance->Fill(fDistance);
|
|---|
| 84 | hLength->Fill(fLength);
|
|---|
| 85 | hWidth->Fill(fWidth);
|
|---|
| 86 | hMCZd->Fill(fMCZd*R2D);
|
|---|
| 87 |
|
|---|
| 88 | if( fabs(fMCZd*R2D - 10.) < 2.0 )
|
|---|
| 89 | {
|
|---|
| 90 | hTableNorm->Fill(log10(fSize),fDistance);
|
|---|
| 91 | hTableEnergy->Fill(log10(fSize),fDistance,log10(fMCEnergy));
|
|---|
| 92 | hTableLength->Fill(log10(fSize),fDistance,fLength);
|
|---|
| 93 | hTableWidth->Fill(log10(fSize),fDistance,fWidth);
|
|---|
| 94 | }
|
|---|
| 95 |
|
|---|
| 96 | }
|
|---|
| 97 |
|
|---|
| 98 | // Average the histogram
|
|---|
| 99 | hTableEnergy->Divide(hTableNorm);
|
|---|
| 100 | hTableLength->Divide(hTableNorm);
|
|---|
| 101 | hTableWidth->Divide(hTableNorm);
|
|---|
| 102 |
|
|---|
| 103 | // We need to smooth somehow
|
|---|
| 104 | //hTableEnergy->Smooth();
|
|---|
| 105 | //hTableLength->Smooth();
|
|---|
| 106 | //hTableWidth->Smooth();
|
|---|
| 107 |
|
|---|
| 108 | // Plot stuff
|
|---|
| 109 | TCanvas *c1 = new TCanvas("c1","c1",0,0,500,500);
|
|---|
| 110 | c1->SetLogy();
|
|---|
| 111 | hMCEnergy->Draw();
|
|---|
| 112 |
|
|---|
| 113 | TCanvas *c2 = new TCanvas("c2","c2",0,0,500,500);
|
|---|
| 114 | c2->SetLogy();
|
|---|
| 115 | hLength->Draw();
|
|---|
| 116 |
|
|---|
| 117 | TCanvas *c3 = new TCanvas("c3","c3",0,0,500,500);
|
|---|
| 118 | c3->SetLogy();
|
|---|
| 119 | hWidth->Draw();
|
|---|
| 120 |
|
|---|
| 121 | TCanvas *c4 = new TCanvas("c4","c4",0,0,500,500);
|
|---|
| 122 | hMCZd->Draw();
|
|---|
| 123 |
|
|---|
| 124 | TCanvas *c5 = new TCanvas("c5","c5",0,0,500,500);
|
|---|
| 125 | hDistance->Draw();
|
|---|
| 126 |
|
|---|
| 127 | TCanvas *c6 = new TCanvas("c6","c6",0,0,500,500);
|
|---|
| 128 | c6->SetLogy();
|
|---|
| 129 | hSize->Draw();
|
|---|
| 130 |
|
|---|
| 131 | TCanvas *c7 = new TCanvas("c7","c7",0,0,500,500);
|
|---|
| 132 | hTableEnergy->Draw("zcol");
|
|---|
| 133 |
|
|---|
| 134 | TCanvas *c8 = new TCanvas("c8","c8",0,0,500,500);
|
|---|
| 135 | hTableLength->Draw("zcol");
|
|---|
| 136 |
|
|---|
| 137 | TCanvas *c9 = new TCanvas("c9","c9",0,0,500,500);
|
|---|
| 138 | hTableWidth->Draw("zcol");
|
|---|
| 139 |
|
|---|
| 140 |
|
|---|
| 141 | //Output
|
|---|
| 142 | if( bOut )
|
|---|
| 143 | {
|
|---|
| 144 | fOutFile->cd();
|
|---|
| 145 | hTableNorm->Write();
|
|---|
| 146 | hTableEnergy->Write();
|
|---|
| 147 | hTableLength->Write();
|
|---|
| 148 | hTableWidth->Write();
|
|---|
| 149 |
|
|---|
| 150 | fOutFile->Close();
|
|---|
| 151 | }
|
|---|
| 152 |
|
|---|
| 153 | }
|
|---|
