EIC-IR/macros_plots/plot_compare_fun4all_eicroot.C

0001 #include "TraceBox.h"
0002
0003 int
0004 plot_compare_fun4all_eicroot()
0005 {
0006   /**** Chose input files ****/
0007   /* File with IR mangets configuration*/
0008   TString fname_irmag("example/proton-magnets-250GeV-opt2.dat");
0009   /* File with Fun4All output*/
0010   TString fname_fun4all("example/eRHIC_proton-magnets-250GeV-opt2_250GeV_0mrad.root");
0011   /* File with EICROOT output*/
0012   TString fname_eicroot("example/eicroot-track_proton-magnets-250GeV-opt2_250GeV_0mrad.txt");
0013
0014   /* Open iput file with trajectories from GEANT4 */
0015   TFile *fin = new TFile(fname_fun4all);
0016
0017   /* Get tree from file */
0018   TTree *tin = (TTree*)fin->Get("T");
0019
0020   int nhits = 0;
0021   tin->SetBranchAddress("n_G4HIT_FWDDISC",&nhits);
0022
0023   /* create graph of particle trajectory */
0024   /* Use only first event (for now) */
0025   tin->GetEntry(0);
0026   cout << "hits: " << nhits << endl;
0027
0028   tin->Draw("G4HIT_FWDDISC.x:G4HIT_FWDDISC.z","Entry$==0","");
0029
0030   TGraph* g1 = new TGraph(nhits*2, &(tin->GetV2()[0]), &(tin->GetV1()[0]));
0031   g1->SetMarkerStyle(7);
0032   g1->SetMarkerSize(1);
0033   g1->SetMarkerColor(kRed);
0034
0035   /* Get tree from file */
0036   TTree *tin2 = new TTree();
0037   tin2->ReadFile(fname_eicroot,"x/F:y:z");
0038
0039   int nhits = tin2->GetEntries();
0040   tin2->Draw("x:z","","");
0041
0042   TGraph* g2 = new TGraph(nhits, &(tin2->GetV2()[0]), &(tin2->GetV1()[0]));
0043   g2->SetMarkerStyle(7);
0044   g2->SetMarkerSize(1);
0045   g2->SetMarkerColor(kGreen+1);
0046
0047   /* Create frame histogram for plot */
0048   TH1F *h1 = new TH1F("h1","",10,0,15000);
0049   h1->GetXaxis()->SetRangeUser(0,5000);
0050   h1->GetYaxis()->SetRangeUser(-50,70);
0051   h1->GetXaxis()->SetTitle("Z(cm)");
0052   h1->GetYaxis()->SetTitle("X(cm)");
0053
0054   /* Plot frame histogram */
0055   TCanvas *c1 = new TCanvas();
0056   h1->Draw("AXIS");
0057
0058   /* Read IR configuration file- this needs to go somewhere else using parameters and a .root file to store them */
0059   ifstream irstream(fname_irmag);
0060
0061   if(!irstream.is_open())
0062     {
0063       cout << "ERROR: Could not open IR configuration file " << fname_irmag << endl;
0064       return -1;
0065     }
0066
0067   while(!irstream.eof()){
0068     string str;
0069     getline(irstream, str);
0070     if(str[0] == '#') continue; //for comments
0071     stringstream ss(str);
0072
0073     string name;
0074     double center_z, center_x, center_y, aperture_radius, length, angle, B, gradient;
0075
0076     ss >> name >> center_z >> center_x >> center_y >> aperture_radius >> length >> angle >> B >> gradient;
0077
0078     if ( name == "" ) continue; //for empty lines
0079
0080     // convert units from m to cm
0081     center_x *= 100;
0082     center_y *= 100;
0083     center_z *= 100;
0084     aperture_radius *= 100;
0085     length *= 100;
0086
0087     // define magnet outer radius
0088     float outer_radius = 30.0; // cm
0089
0090     //flip sign of dipole field component- positive y axis in Geant4 is defined as 'up',
0091     // positive z axis  as the hadron-going direction
0092     // in a right-handed coordinate system x,y,z
0093     B *= -1;
0094
0095     // convert angle from millirad to rad
0096     angle = (angle / 1000.);
0097
0098     // Place IR component
0099     cout << "New IR component: " << name << " at z = " << center_z << endl;
0100
0101     string volname = "IRMAGNET_";
0102     volname.append(name);
0103
0104     /* Draw box for magnet position on canvas */
0105     TPolyLine *b1 = TraceBox( angle, center_z, center_x, length, aperture_radius, outer_radius ); //upper box
0106     TPolyLine *b2 = TraceBox( angle, center_z, center_x, length, -1 * aperture_radius, -1 * outer_radius ); //lower box
0107
0108     if(B != 0 && gradient == 0.0){
0109       //dipole magnet
0110       b1->SetFillColor(kOrange+1);
0111       b2->SetFillColor(kOrange+1);
0112     }
0113     else if( B == 0 && gradient != 0.0){
0114       //quad magnet
0115       b1->SetFillColor(kBlue+1);
0116       b2->SetFillColor(kBlue+1);
0117     }
0118     else{
0119       //placeholder magnet
0120       b1->SetFillColor(kGray+1);
0121       b2->SetFillColor(kGray+1);
0122     }
0123     b1->Draw("Fsame");
0124     b2->Draw("Fsame");
0125   }
0126
0127   /* draw particle trajectory */
0128   g2->Draw("LPsame");
0129   g1->Draw("LPsame");
0130
0131   return 0;
0132 }