EIC-IR/macros_plots/plot_track_multi.C

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