Changeset 6493 for trunk/MagicSoft/Mars/macros
- Timestamp:
- 02/15/05 17:33:43 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/macros/collarea.C
r2098 r6493 17 17 ! 18 18 ! Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de> 19 ! Author(s): Abelardo Moralejo, 2/2005 <mailto:moralejo@pd.infn.it> 19 20 ! 20 ! Copyright: MAGIC Software Development, 2000-200 321 ! Copyright: MAGIC Software Development, 2000-2005 21 22 ! 22 23 ! 23 24 \* ======================================================================== */ 24 25 26 // 27 // Example macro on how to calculate effective areas. The input is a file 28 // containing Hillas parameters of the MC events surviving a certain kind 29 // of analysis. It must also contain a tree called "OriginalMC" with a branch 30 // "MMcEvtBasic" and one entry per original simulated shower used to produce 31 // the events in that input file. The current (2/2005) camera simulation in 32 // the cvs already writes out all events with the MMcEvtBasic container to 33 // the Events tree. In Mars/macros/mccalibrate.C and starmc2.C you can see 34 // how this branch is put to a separate tree "OriginalMC" and written out 35 // to the "star" file. This will not work with star files not containing this 36 // extra tree. Please contact moralejo@pd.infn.it in case of questions 37 // about this new approach to the calculation of effective areas. 38 // 25 39 26 void collarea(TString filename=" camera.root", TString outname="")40 void collarea(TString filename="star_gamma_test.root", TString outname="area.root") 27 41 { 28 42 MStatusDisplay *d = new MStatusDisplay; 29 43 // redirect logging output to GUI, kFALSE to disable stream to stdout 30 44 d->SetLogStream(&gLog, kTRUE); 31 45 32 // 33 // first we have to create our empty lists 34 // 35 MParList parlist; 36 MTaskList tasklist; 46 // 47 // First loop: fill the MHMcCollectionArea::fHistAll histogram containing 48 // all the events produced at the Corsika level: 49 // 50 MParList parlist; 51 MTaskList tasklist; 37 52 38 53 parlist.AddToList(&tasklist); 39 54 40 // 41 // Setup out tasks: 42 // - First we have to read the events 43 // - Then we can fill the efficiency histograms 44 // 45 MReadMarsFile reader("Events", filename); 46 MMcCollectionAreaCalc effi; 55 // 56 // Setup out tasks: 57 // - First we have to read the events 58 // 59 MReadMarsFile reader("OriginalMC", filename); 60 reader.DisableAutoScheme(); 47 61 48 tasklist.AddToList(&reader); 49 tasklist.AddToList(&effi); 62 // 63 // Create collection area object and add it to the task list 64 // 65 MHMcCollectionArea collarea; 66 parlist.AddToList(&collarea); 50 67 51 // 52 // set up the loop for the processing 53 // 54 MEvtLoop magic; 55 magic.SetParList(&parlist); 68 // 69 // Set the COARSE binnings in which you will divide your data, to get 70 // the collection area in these same binnings. You can set also this 71 // binnings on an existing and filled MHMcCollectionAreaObject, and 72 // call Calc again to get the effective area in other coarse bins. 73 // 74 MBinning binsTheta("binsTheta"); 75 MBinning binsEnergy("binsEnergy"); // name must be standardized! 76 Int_t nbins = 8; 77 TArrayD edges(nbins+1); 78 edges[0] = 0; 79 edges[1] = 10.; 80 edges[2] = 20.; 81 edges[3] = 30.; 82 edges[4] = 40.; 83 edges[5] = 50.; 84 edges[6] = 60.; 85 edges[7] = 70.; 86 edges[8] = 80.; 87 binsTheta.SetEdges(edges); // Theta [deg] 88 binsEnergy.SetEdgesLog(20, 2., 20000); // Energy [GeV] 89 parlist.AddToList(&binsTheta); 90 parlist.AddToList(&binsEnergy); 56 91 57 // 58 // Start to loop over all events 59 // 60 magic.SetDisplay(d); 61 if (!magic.Eventloop()) 62 return; 92 // Task which fills the necessary histograms in MHMcCollectionArea 93 MMcCollectionAreaCalc effi; 63 94 64 tasklist.PrintStatistics(); 95 // Tentative energy spectrum. This may slightly modify the calculation 96 // of effective areas. If not added to the parameter list, a flat spectrum 97 // in dN/dE will be assumed. "x" stands for the energy. The name of the 98 // function must be standardized! 65 99 66 // 67 // Now the histogram we wanted to get out of the data is 68 // filled and can be displayed 69 // 70 if ((d = magic.GetDisplay())) 71 d->AddTab("Collection Area"); 72 else 73 new TCanvas("CollArea", "Result of the collection area calculation"); 100 TF1 Spectrum("Spectrum", "pow(x,-5.)"); 101 effi.SetSpectrum(&Spectrum); 74 102 75 parlist.FindObject("MHMcCollectionArea")->DrawClone("nonew"); 103 tasklist.AddToList(&reader); 104 tasklist.AddToList(&effi); 105 106 // 107 // set up the loop for the processing 108 // 109 MEvtLoop magic; 110 magic.SetParList(&parlist); 111 112 // 113 // Start to loop over all events 114 // 115 magic.SetDisplay(d); 116 117 if (!magic.Eventloop()) 118 return; 119 120 tasklist.PrintStatistics(); 76 121 77 122 78 // 79 // Write histogram to a file in case an output 80 // filename has been supplied 81 // 82 if (outname.IsNull()) 83 return; 123 // 124 // Second loop: now fill the histogram MHMcCollectionArea::fHistSel 125 // 84 126 85 TFile f(outname, "recreate"); 86 if (!f) 87 return; 127 MParList parlist2; 128 MTaskList tasklist2; 129 parlist2.AddToList(&tasklist2); 130 parlist2.AddToList((MHMcCollectionArea*)parlist->FindObject("MHMcCollectionArea")); 88 131 89 parlist.FindObject("MHMcCollectionArea")->Write(); 132 MReadMarsFile reader2("Events", filename); 133 reader2.DisableAutoScheme(); 134 135 tasklist2.AddToList(&reader2); 136 137 // The same task as before, now will fill MHMcCollectionArea::fHistSel 138 // and call MHMcCollectionArea::Calc in the PostProcess 139 tasklist2.AddToList(&effi); 140 141 // 142 // set up the loop for the processing 143 // 144 MEvtLoop magic2; 145 magic2.SetParList(&parlist2); 146 147 // 148 // Start to loop over all events 149 // 150 magic2.SetDisplay(d); 151 152 if (!magic2.Eventloop()) 153 return; 154 155 tasklist2.PrintStatistics(); 156 157 // 158 // Now the histogram we wanted to get out of the data is 159 // filled and can be displayed 160 // 161 if ((d = magic2.GetDisplay())) 162 d->AddTab("Collection Area"); 163 else 164 new TCanvas("CollArea", "Result of the collection area calculation"); 165 166 gStyle->SetPalette(1,0); 167 gPad->SetLogx(); 168 TH2D* areaplot = ((MHMcCollectionArea*)parlist2.FindObject 169 ("MHMcCollectionArea"))->GetHistCoarse(); 170 171 areaplot->SetStats(0); 172 areaplot->SetTitle("Effective area (m^{2})"); 173 areaplot->GetXaxis()->SetTitle("Energy (GeV)"); 174 areaplot->GetYaxis()->SetTitle("\\theta (deg)"); 175 areaplot->DrawCopy("zcol"); 176 177 // 178 // Write histogram to a file in case an output 179 // filename has been supplied 180 // 181 182 if (outname.IsNull()) 183 return; 184 185 TFile f(outname, "recreate"); 186 if (!f) 187 return; 188 189 // Write to the output file: 190 parlist2.FindObject("MHMcCollectionArea")->Write(); 90 191 }
Note:
See TracChangeset
for help on using the changeset viewer.