drich-dev/scripts/test_tracks.C

0001 // draw track points of projected trajectories on dRICH planes
0002 //
0003
0004
0005 class radiator_config {
0006   public:
0007     radiator_config(int id_, std::string name_, Color_t color_, Style_t style_, TTreeReader& reader) :
0008       id(id_),
0009       name(name_),
0010       color(color_),
0011       style(style_),
0012       collName(Form("DRICH%s_0",name.c_str())),
0013       x_arr({reader, collName+".position.x"}),
0014       y_arr({reader, collName+".position.y"}),
0015       z_arr({reader, collName+".position.z"}) {}
0016     ~radiator_config() {}
0017     int         id;
0018     std::string name;
0019     TString     collName;
0020     Color_t     color;
0021     Style_t     style;
0022     TTreeReaderArray<Float_t> x_arr, y_arr, z_arr;
0023 };
0024
0025
0026
0027 void test_tracks(TString infileN = "out/rec.edm4hep.root") {
0028   auto infile = new TFile(infileN);
0029   auto tr = (TTree*) infile->Get("events");
0030   TTreeReader tr_reader("events", infile);
0031
0032   std::vector<radiator_config*> radiators = {
0033     new radiator_config( 0, "AerogelTracks", kAzure+7, kFullCircle, tr_reader ),
0034     new radiator_config( 1, "GasTracks",     kRed,     kFullCircle, tr_reader ),
0035     // new radiator_config( 2, "MergedTracks",  kGreen+1, kFullCircle, tr_reader ),
0036   };
0037
0038   double rmax = 1900;
0039   double xmax = rmax;
0040   double ymax = rmax;
0041   double zmin = 1800;
0042   double zmax = 3300;
0043
0044   // ymax = 200; // zoom to horizontal (y=0) plane
0045
0046   // 2D plots
0047   enum views { zx, zy, xy, Nview };
0048   TString viewN[Nview], viewT[Nview];
0049   viewN[zx] = "zx";
0050   viewN[zy] = "zy";
0051   viewN[xy] = "xy";
0052   viewT[zx] = "Top View (toward -y);z [mm];x [mm]";
0053   viewT[zy] = "Side View (toward +x);z [mm];y [mm]";
0054   viewT[xy] = "Front View (toward +z);x [mm];y [mm]";
0055   std::vector<TGraph*> points2[Nview];
0056   TMultiGraph *points2mgr[Nview];
0057   for(int view=0; view<Nview; view++) {
0058     points2mgr[view] = new TMultiGraph();
0059     points2mgr[view]->SetName("points_"+viewN[view]);
0060     points2mgr[view]->SetTitle(viewT[view]);
0061   }
0062   for(auto& radiator : radiators) {
0063     for(int view=0; view<Nview; view++) {
0064       auto gr = new TGraph();
0065       points2[view].push_back(gr);
0066       points2mgr[view]->Add(gr);
0067       gr->SetMarkerColor(radiator->color);
0068       gr->SetMarkerStyle(radiator->style);
0069       gr->SetName(Form("points_%s_%s",viewN[view].Data(),radiator->name.c_str()));
0070       gr->SetTitle(viewT[view]);
0071       tr_reader.Restart();
0072       while(tr_reader.Next()) {
0073         auto i_offset = gr->GetN();
0074         for(int i=0; i<radiator->x_arr.GetSize(); i++) {
0075           auto x = radiator->x_arr[i];
0076           auto y = radiator->y_arr[i];
0077           auto z = radiator->z_arr[i];
0078           switch(view) {
0079             case zx: gr->SetPoint(i+i_offset,z,x); break;
0080             case zy: gr->SetPoint(i+i_offset,z,y); break;
0081             case xy: gr->SetPoint(i+i_offset,x,y); break;
0082           }
0083         }
0084       }
0085     }
0086   }
0087   auto canv2 = new TCanvas("canv2","canv2",Nview*800,600);
0088   canv2->Divide(Nview,1);
0089   for(int view=0; view<Nview; view++) {
0090     canv2->cd(view+1);
0091     canv2->GetPad(view+1)->SetGrid(0,1);
0092     canv2->GetPad(view+1)->SetLeftMargin(0.15);
0093     TString mode = "AP";
0094     if(view==xy) mode+=" RX";
0095     points2mgr[view]->Draw(mode);
0096   }
0097   points2mgr[zx]->GetXaxis()->SetLimits(     zmin, zmax );
0098   points2mgr[zx]->GetYaxis()->SetRangeUser( -xmax, xmax );
0099   points2mgr[zy]->GetXaxis()->SetLimits(     zmin, zmax );
0100   points2mgr[zy]->GetYaxis()->SetRangeUser( -ymax, ymax );
0101   points2mgr[xy]->GetXaxis()->SetLimits(    -xmax, xmax );
0102   points2mgr[xy]->GetYaxis()->SetRangeUser( -ymax, ymax );
0103   canv2->SaveAs("out/tracks.png");
0104
0105   // 3D plots (NOTE: point positions inaccurate, since 3D histogram bins)
0106   std::vector<TH3D*> points3;
0107   for(auto& radiator : radiators) {
0108     TString histName = Form("%s_hist",radiator->name.c_str());
0109     points3.push_back(
0110         new TH3D(
0111           histName,
0112           "Propagated track points",
0113           100, zmin,  zmax,
0114           100, -xmax, xmax,
0115           100, -ymax, ymax
0116           )
0117         );
0118     points3.back()->SetMarkerColor(radiator->color);
0119     points3.back()->SetMarkerStyle(kFullCircle);
0120     tr->Project(
0121         histName,
0122         radiator->collName+".position.y:"+radiator->collName+".position.x:"+radiator->collName+".position.z"
0123         );
0124   }
0125   auto canv3 = new TCanvas();
0126   for(auto plot : points3) plot->Draw("same");
0127 }
0128