diff --git a/macro/mvd/ana_radlength_mvd.C b/macro/mvd/ana_radlength_mvd.C deleted file mode 100644 index 8792d3b32442844fa14f0ba66c36744caa31e3a6..0000000000000000000000000000000000000000 --- a/macro/mvd/ana_radlength_mvd.C +++ /dev/null @@ -1,153 +0,0 @@ -#if !defined(__CINT__) || defined(__MAKECINT__) - - -#include "TString.h" -#include "TFile.h" -#include "TTree.h" -#include "TClonesArray.h" -#include "TH1.h" -#include "TVector3.h" -#include "TCanvas.h" -#include "TSystem.h" -#include "TStyle.h" -#include "TROOT.h" -#include "TProfile2D.h" - -#include "/cvmfs/fairroot.gsi.de/fairroot/v-15.03a_fairsoft-mar15p2/include/FairRadLenPoint.h" - -#include <iostream> -#include <vector> -using std::cout; -using std::endl; -using std::vector; -#endif - -Int_t ana_radlength_mvd() -{ - - Int_t nEvents=10000000; - TString inFile = "data/radlength.mvd_v18a.mc.root"; - TFile* f = new TFile(inFile); - TTree *t=(TTree *) f->Get("cbmsim") ; - - TString dir = gSystem->Getenv("VMCWORKDIR"); - TString file = dir + "/gconfig/basiclibs.C"; - gROOT->LoadMacro(file); - gROOT->ProcessLine("basiclibs()"); - gSystem->Load("libGeoBase"); - gSystem->Load("libParBase"); - gSystem->Load("libBase"); - gSystem->Load("libCbmBase"); - - - TClonesArray* radlen_array = new TClonesArray("FairRadLenPoint"); - t->SetBranchAddress("RadLen", &radlen_array); - - TH1D* hisRadLen = new TH1D("hisRadLen","Radiation Length", 1000,0,100); - - const int NStations = 4; - const int NBins = 1500; - int RMax = 11; - TProfile2D* hStaRadLen[NStations]; - - for ( int i = 0; i < NStations; ++i ) { - TString name = "Radiation Thickness [%],"; - name += " Station"; - name += i; - if ( i > 1) RMax = 15; - hStaRadLen[i] = new TProfile2D(name, name, NBins,-RMax,RMax, NBins,-RMax,RMax); - } - - int startEvent = 0; - - TVector3 vecs,veco; - std::map<int,int> trackHitMap; - - for (Int_t j=startEvent; j<(nEvents+startEvent) && j<t->GetEntriesFast(); j++) { - t->GetEntry(j); - if ( 0 == j%10000 ) { - cout<<">>>> Event No "<<j<<endl; - } - - - const int NTracks = 1; - std::vector<double> RadLengthOnTrack (NTracks,0.0); //trackID, vector with points on track - std::vector<double> TrackLength (NTracks,0.0); //trackID, vector with points on track - - vector<double> RadThick(NStations,0); - double x,y; - for (Int_t i=0; i<radlen_array->GetEntriesFast(); i++) { -// cout<<"Point No "<<i<<endl; - - FairRadLenPoint *point=(FairRadLenPoint*)radlen_array->At(i); - -// cout << "Track ID: " << point->GetTrackID() << std::endl; - - TVector3 pos, posOut, res; - pos = point->GetPosition(); - posOut = point->GetPositionOut(); - res = posOut - pos; - - const TVector3 middle2 = (posOut + pos); - x = middle2.X()/2; - y = middle2.Y()/2; - const double z = middle2.Z()/2; - - const double radThick = res.Mag()/point->GetRadLength(); - RadLengthOnTrack[point->GetTrackID()] += radThick; - TrackLength[point->GetTrackID()] += res.Mag(); - - // cout << radThick << endl; - // cout << pos.X() << " " << pos.Y() << " " << pos.Z() << " " << endl; - // cout << posOut.X() << " " << posOut.Y() << " " << posOut.Z() << " " << endl; - - - int iSta = pos.Z()/5 - 1 + 0.5; - int iStaOut = posOut.Z()/5 - 1 + 0.5; - // int iSta = pos.Z()/10 - 3 + 0.5; // suppose equdistant stations at 30-100 cm - // int iStaOut = posOut.Z()/10 - 3 + 0.5; - if ( iSta != iStaOut ) continue; - if ( iSta >= NStations || iSta < 0 ) continue; - RadThick[iSta] += radThick; - - } - - for ( int i = 0; i < NStations; ++i ) { - // cout << i << " " << x << " " << y << " " << RadThick[i] << endl; - hStaRadLen[i]->Fill( x, y, RadThick[i]*100 ); - } - - for (int k = 0; k < RadLengthOnTrack.size(); k++) { - if (RadLengthOnTrack[k] > 0) { -// std::cout << "Full TrackLength: " << TrackLength[k] << std::endl; -// std::cout << "Full RadiationLength: " << RadLengthOnTrack[k] << std::endl; - hisRadLen->Fill(RadLengthOnTrack[k]); //% - } - } - } - - TCanvas* can1 = new TCanvas(); - can1->Divide(NStations/2,2); - gStyle->SetPalette(1); - gStyle->SetOptStat(0); - - TFile *rtFile = new TFile("mvd_matbudget_v18a.root","RECREATE"); - rtFile->cd(); - - for ( int i = 0; i < NStations; ++i ) { - can1->cd(i+1); - hStaRadLen[i]->GetXaxis()->SetTitle("x [cm]"); - hStaRadLen[i]->GetYaxis()->SetTitle("y [cm]"); - //hStaRadLen[i]->GetZaxis()->SetTitle("radiation thickness [%]"); - hStaRadLen[i]->SetAxisRange(0, 2, "Z"); - hStaRadLen[i]->Draw("colz"); - hStaRadLen[i]->Write(); - } - can1->SaveAs("mvd_matbudget_v18a.pdf"); - can1->SaveAs("mvd_matbudget_v18a.eps"); - can1->SaveAs("mvd_matbudget_v18a.png"); - - cout << endl << "i am out, pease!" << endl; - return 0; -} - diff --git a/macro/mvd/run_radlength_mvd.C b/macro/mvd/run_radlength_mvd.C deleted file mode 100644 index 4ed86397caeeb946744aa5d1abdf32c6963d3fa8..0000000000000000000000000000000000000000 --- a/macro/mvd/run_radlength_mvd.C +++ /dev/null @@ -1,154 +0,0 @@ -// -------------------------------------------------------------------------- -// -// Macro for standard transport simulation using UrQMD input and GEANT3 -// Standard CBM setup with MVD, STS, RICH, TRD, TOF and ECAL -// -// V. Friese 22/02/2007 -// -// -------------------------------------------------------------------------- - -void run_radlength_mvd(Int_t nEvents = 10000000) -{ - - // ======================================================================== - // Adjust this part according to your requirements - - // ----- Paths and file names -------------------------------------------- - TString inDir = gSystem->Getenv("VMCWORKDIR"); - TString inFile = inDir + "/input/urqmd.ftn14"; - TString outDir = "data"; - TString outFile = outDir + "/radlength.mvd_v15a_vacuum.mc.root"; - TString parFile = outDir + "/params.mvd_v15a_vacuum.root"; - - // ----- Geometries ----------------------------------------------------- - TString caveGeom = "cave.geo"; - - TString pipeGeom = "pipe/pipe_v14l.root"; - TString mvdGeom = "mvd/MVD_v15a_vacuum.root"; - - - // Magnet geometry and field map - TString magnetGeom = "magnet/magnet_v12b.geo.root"; - TString fieldMap = "field_v12b"; // name of field map - Int_t fieldZ = 40.; // field centre z position - Int_t fieldScale = 1.; // field scaling factor - Int_t fieldSymType = 3; - - - // In general, the following parts need not be touched - // ======================================================================== - - - - - // ---- Debug option ------------------------------------------------- - gDebug = 0; - // ------------------------------------------------------------------------ - - - - // ----- Timer -------------------------------------------------------- - TStopwatch timer; - timer.Start(); - // ------------------------------------------------------------------------ - - - // ----- Create simulation run ---------------------------------------- - FairRunSim* fRun = new FairRunSim(); - fRun->SetName("TGeant3"); // Transport engine - fRun->SetOutputFile(outFile); // Output file - FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); - // ------------------------------------------------------------------------ - - - // ----- Create media ------------------------------------------------- - fRun->SetMaterials("media.geo"); // Materials - // ------------------------------------------------------------------------ - - - // ----- Create detectors and passive volumes ------------------------- - if ( caveGeom != "" ) - FairModule* cave = new CbmCave("CAVE"); - cave->SetGeometryFileName(caveGeom); - fRun->AddModule(cave); - - FairModule* pipe= new CbmPipe("PIPE"); - pipe->SetGeometryFileName(pipeGeom); - fRun->AddModule(pipe); - - FairDetector* mvd = new CbmMvd("MVD", kTRUE); - mvd->SetGeometryFileName(mvdGeom); - mvd->SetMotherVolume("pipevac1"); - fRun->AddModule(mvd); - - // ------------------------------------------------------------------------ - - - - // ----- Create magnetic field --------------------------------------- - if ( 2 == fieldSymType ) { - CbmFieldMap* magField = new CbmFieldMapSym2(fieldMap); - } else if ( 3 == fieldSymType ) { - CbmFieldMap* magField = new CbmFieldMapSym3(fieldMap); - } - magField->SetPosition(0., 0., fieldZ); - magField->SetScale(fieldScale); - fRun->SetField(magField); - // ------------------------------------------------------------------------ - - // ----- Create PrimaryGenerator -------------------------------------- - FairPrimaryGenerator* primGen = new FairPrimaryGenerator(); - fRun->SetGenerator(primGen); - - //ROOTino - FairBoxGenerator *fBox1 = new FairBoxGenerator(0, 1); - fBox1->SetBoxXYZ(-15.,-15.,15.,15.,0.); - fBox1->SetPRange(0.1,0.5); - fBox1->SetThetaRange(0.,0.); - fBox1->SetPhiRange(0.,360.); - primGen->AddGenerator(fBox1); - - fRun->SetStoreTraj(kFALSE); - fRun->SetRadLenRegister(kTRUE); - // ------------------------------------------------------------------------ - - // ----- Run initialisation ------------------------------------------- - fRun->Init(); - // ------------------------------------------------------------------------ - - // ----- Runtime database --------------------------------------------- - CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar"); - fieldPar->SetParameters(magField); - fieldPar->setChanged(); - fieldPar->setInputVersion(fRun->GetRunId(),1); - Bool_t kParameterMerged = kTRUE; - FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged); - parOut->open(parFile.Data()); - rtdb->setOutput(parOut); - rtdb->saveOutput(); - rtdb->print(); - // ------------------------------------------------------------------------ - - - // ----- Start run ---------------------------------------------------- - fRun->Run(nEvents); - // ------------------------------------------------------------------------ - fRun->CreateGeometryFile("data/geo_rad_mvd15a_vacuum.root"); - - - // ----- 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 << endl; - // ------------------------------------------------------------------------ - - cout << " Test passed" << endl; - cout << " All ok " << endl; -} -