ForwardTracking/macros/DrawFieldMap.C

0001 // $Id: $
0002
0003 /*!
0004  * \file DrawFieldMap.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "SaveCanvas.C"
0012
0013 void DrawFieldMap(const TString sfield_map = "/phenix/upgrades/decadal/fieldmaps/fsPHENIX.2d.root")
0014 //void DrawFieldMap(const TString sfield_map = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root")
0015 {
0016   gROOT->LoadMacro("sPhenixStyle.C");
0017   SetsPhenixStyle();
0018
0019   TFile *_file0 = TFile::Open(sfield_map);
0020   TTree *fieldmap = (TTree *) _file0->GetObjectChecked("fieldmap", "TTree");
0021
0022   TH2F *hfield = new TH2F("hfield", "hfield", 480 / 2, -30, 450, 280 / 2, 0, 280);
0023   TH2F *hbendingfield = hfield->Clone("hbendingfield");
0024
0025   fieldmap->SetAlias("perb_angle", "atan2(r,z) + pi/2");
0026   fieldmap->SetAlias("bendingfield", "br*cos(perb_angle) + bz*sin(perb_angle)");
0027
0028   fieldmap->Draw("r:z>>hfield", "sqrt(br**2 + bz**2)", "goff");
0029   fieldmap->Draw("r:z>>hbendingfield", "bendingfield", "goff");
0030
0031   TCanvas *c1 = new TCanvas(TString("DrawFieldMap_hfield_") + gSystem->BaseName(sfield_map),
0032                             TString("DrawFieldMap_hfield_") + gSystem->BaseName(sfield_map), 1800,
0033                             900);
0034   c1->Divide(1, 1);
0035   int idx = 1;
0036   TPad *p;
0037   p = (TPad *) c1->cd(idx++);
0038   c1->Update();
0039
0040   p->SetRightMargin(0.2);
0041
0042   hfield->GetZaxis()->SetRangeUser(0., 2.5);
0043   hfield->SetTitle("#it{#bf{sPHENIX}} Field Map;z [cm];R [cm];Field strength [T]");
0044   hfield->GetYaxis()->SetTitleOffset(0.9);
0045   hfield->Draw("colz");
0046
0047   //  %h-going ECAL
0048   //  \newcommand{ \hecalfsphenix }
0049   //  {
0050   //
0051   //    % towerMap_FEMC_v005.txt
0052   //    \draw[fill=blizzardblue] ( 3.28250, 1.82655 ) -- ( 3.28250, 0.1 ) -- ( 2.91750, 0.1 ) -- ( 2.91750, 1.82655 ) -- ( 3.28250, 1.82655 );
0053   //
0054   //  }
0055   TBox *b = new TBox(2.91750 * 100, 0.1 * 100, 3.2 * 100, 1.82655 * 100);
0056   b->SetFillColor(0);
0057   b->SetFillStyle(0);
0058   b->SetLineColor(kGray);
0059   b->SetLineWidth(4);
0060   b->Draw();
0061   //
0062   //  %h-going HCAL
0063   //  \newcommand{ \hhcalfsphenix }
0064   //  {
0065   //
0066   //    % towerMap_FHCAL_v004.txt
0067   //  \draw[fill=bluegray] ( 4.5, 2.62 ) -- ( 4.5, 0.05 ) -- ( 3.5, 0.05 ) -- ( 3.5, 2.62 ) -- ( 4.5, 2.62 );
0068   //
0069   //  }
0070   TBox *b = new TBox(3.5 * 100, 0.05 * 100, 4.5 * 100, 2.62 * 100);
0071   b->SetFillColor(0);
0072   b->SetFillStyle(0);
0073   b->SetLineColor(kGray);
0074   b->SetLineWidth(4);
0075   b->Draw();
0076
0077   //  %h-going GEM fsPhenix
0078   //  \newcommand{ \hgemfsphenix }
0079   //  {
0080   //
0081   //  \draw[ultra thick,bazaar] ( 1.2, 0.1 ) -- ( 1.2, 0.8 );
0082   //  \draw[ultra thick,bazaar] ( 1.5, 0.1 ) -- ( 1.5, 1 );
0083   //  \draw[ultra thick,bazaar] ( 2.7, 0.1 ) -- ( 2.7, 0.8 ) -- ( 2.5, 1.3);
0084   //
0085   //  }
0086   TLine *l = new TLine(1.2 * 100, 0.1 * 100, 1.2 * 100, 0.8 * 100);
0087   l->SetLineColor(kGray);
0088   l->SetLineWidth(4);
0089   l->Draw();
0090   TLine *l = new TLine(1.5 * 100, 0.1 * 100, 1.5 * 100, 1 * 100);
0091   l->SetLineColor(kGray);
0092   l->SetLineWidth(4);
0093   l->Draw();
0094   TLine *l = new TLine(2.7 * 100, 0.1 * 100, 2.7 * 100, 0.8 * 100);
0095   l->SetLineColor(kGray);
0096   l->SetLineWidth(4);
0097   l->Draw();
0098   TLine *l = new TLine(2.7 * 100, 0.8 * 100, 2.5 * 100, 1.3 * 100);
0099   l->SetLineColor(kGray);
0100   l->SetLineWidth(4);
0101   l->Draw();
0102
0103   TBox *b = new TBox(-30, 0.02 * 100, 100, 78);
0104   b->SetFillColor(0);
0105   b->SetFillStyle(0);
0106   b->SetLineColor(kGray);
0107   b->SetLineWidth(4);
0108   b->Draw();
0109
0110   SaveCanvas(c1, TString(c1->GetName()), kTRUE);
0111   //  c1->Print(TString(c1->GetName())+".pdf");
0112   c1->Print(TString(c1->GetName()) + ".eps");
0113   //  c1->Print(TString(c1->GetName())+".svg");
0114
0115   //  TCanvas *c1 = new TCanvas("DrawFieldMap_hbendingfield", "DrawFieldMap_hbendingfield", 1800,
0116   //                            900);
0117   //  c1->Divide(1, 1);
0118   //  int idx = 1;
0119   //  TPad *p;
0120   //  p = (TPad *) c1->cd(idx++);
0121   //  c1->Update();
0122   //
0123   //  p->SetRightMargin(0.1);
0124   //
0125   //  hbendingfield->Draw("colz");
0126
0127   //  SaveCanvas(c1, TString(c1->GetName()), kTRUE);
0128
0129   // output field map
0130
0131   const int n_data = fieldmap->Draw("r:z:br:bz", "", "goff");
0132
0133   string file_name = string(c1->GetName()) + ".data";
0134   fstream foutput(file_name.c_str(), ios_base::out);
0135
0136   for (int i = 0; i < n_data; ++i)
0137   {
0138     foutput << (fieldmap->GetV1())[i] << "\t";
0139     foutput << (fieldmap->GetV2())[i] << "\t";
0140     foutput << (fieldmap->GetV3())[i] << "\t";
0141     foutput << (fieldmap->GetV4())[i] << "\n";
0142   }
0143
0144   foutput.close();
0145
0146   TFile *_file1 = TFile::Open(TString("ReverseBr_") + gSystem->BaseName(sfield_map), "recreate");
0147
0148   TNtuple *ReverseBr_fieldmap =
0149       new TNtuple("fieldmap",
0150                   "Field Map for fsPHENIX with corrected Br",
0151                   "z:r:bz:br");
0152
0153   for (int i = 0; i < n_data; ++i)
0154   {
0155     ReverseBr_fieldmap->Fill((fieldmap->GetV2())[i], (fieldmap->GetV1())[i], (fieldmap->GetV4())[i], -(fieldmap->GetV3())[i]);
0156   }
0157   ReverseBr_fieldmap->Write();
0158   _file1->Close();
0159 }