//// root macro to analyze the simulation output //#include //#include //#include //#include // //// #include //#include //#include //#include //#include //#include //#include // for global root variables //#include //#include //#include //#include //#include //#include //#include //#include //#include // ////#include "$VMCWORKDIR/base/FairRadLenPoint.h" //#include "/Users/ralfk/Panda/pandaroot/base/FairRadLenPoint.h" void materialana(int nEvents = 1000, bool verbose = false) { // ----- Load libraries ------------------------------------------------ gROOT->Macro("$VMCWORKDIR/example/gconfig/basiclibs.C"); gSystem->Load("libProof"); gSystem->Load("libFairTools"); gSystem->Load("libFairDB"); gSystem->Load("libGeoBase"); gSystem->Load("libParBase"); gSystem->Load("libBase"); gSystem->Load("libMCStack"); //gSystem->Load("libGen"); //gSystem->Load("libfield"); // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // Fetch data arrays std::string inFile = "eic_rad_length.mc.root"; TFile* f = new TFile(inFile.c_str()); // the sim file you want to analyse TTree *t=(TTree *) f->Get("cbmsim"); TClonesArray* mc_array=new TClonesArray("FairMCTrack"); t->SetBranchAddress("MCTrack",&mc_array);//Branch names TClonesArray* rad_array=new TClonesArray("FairRadLenPoint"); t->SetBranchAddress("RadLen",&rad_array); // Setup Histograms Double_t boxwidth=300; TH1D* hisx = new TH1D("hX","X",100,-boxwidth,boxwidth); TH1D* hisy = new TH1D("hY","Y",100,-boxwidth,boxwidth); TH1D* hisz = new TH1D("hZ","Z",100,-boxwidth,boxwidth); TProfile* thetaprofile = new TProfile("thetaprof","Radiaion length (#theta)",180,0.,TMath::Pi()); TProfile* phiprofile = new TProfile("phiprof","Radiaion length (#phi)",180,-1.*TMath::Pi(),TMath::Pi()); TProfile2D* histhephi = new TProfile2D("hThePhi","",90,0.2,2.8,90,-1.*TMath::Pi(),TMath::Pi()); histhephi->SetTitle("Radiation length map #theta & #phi;#theta;#phi"); //TString detname,volname; double radlen=0.,effradl=0.,effradlsum=0.,theta=0.,phi=0.; TVector3 in, out, dist, point; TVector3 mom; for (Int_t event=0; eventGetEntriesFast(); event++) { t->GetEntry(event); if(verbose) cout<<"Event No "<GetEntriesFast();trackno++){ FairMCTrack* aTrack = (FairMCTrack*)mc_array->At(trackno); aTrack->GetMomentum(mom); // write momentum to mom theta = mom.Theta(); phi = mom.Phi(); effradl=0.; effradlsum=0.; for (Int_t k=0; kGetEntriesFast(); k++){ FairRadLenPoint* radpoint = (FairRadLenPoint*)rad_array->At(k); if (radpoint->GetTrackID() != trackno) continue; radlen = radpoint->GetRadLength(); in = radpoint->GetPosition(); out = radpoint->GetPositionOut(); dist = in - out; point.SetXYZ(0.5*(in.x()+out.x()),0.5*(in.y()+out.y()),0.5*(in.z()+out.z())); if(in.Mag() < 0.02) continue; // cut target effradl = dist.Mag()/radlen; effradlsum += effradl; if(event < 1) // testing output printf("Test: phi(%f), theta(%f), radlen(%g), effradl(%g)\n",phi,theta,radlen,effradl); hisx->Fill(in.X(),effradl); hisy->Fill(in.Y(),effradl); hisz->Fill(in.Z(),effradl); }//radpoints thetaprofile->Fill(theta,effradlsum); phiprofile->Fill(phi,effradlsum); histhephi->Fill(theta,phi,effradlsum); }//tracks }// end for event // Draw the Histos Int_t a=3,b=2; int resol = 250; TCanvas* can1 = new TCanvas("can1","MCHit view in MVD",0,0,a*resol,b*resol); can1->Divide(a,b); can1->cd(1); hisx->SetFillColor(38); hisx->DrawCopy(); can1->cd(2); hisy->SetFillColor(45); hisy->DrawCopy(); can1->cd(3); hisz->SetFillColor(30); hisz->DrawCopy(); can1->cd(4); gPad->SetLogy(); thetaprofile->SetFillColor(30); thetaprofile->SetLineColor(30); thetaprofile->Draw("hist"); can1->cd(5); gPad->SetLogy(); phiprofile->SetFillColor(45); phiprofile->SetLineColor(45); phiprofile->Draw("hist"); can1->cd(6); histhephi->SetStats(false); histhephi->DrawCopy("colz"); // ----- Finish ------------------------------------------------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << endl << endl; cout << "Macro finished succesfully." << endl; //cout << "Output file is " << outFile << endl; //cout << "Parameter file is " << parFile << endl; cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl; cout << endl; // ------------------------------------------------------------------------ }