// Macro created A.Sanchez // It creates a geant simulation file for hyp { TStopwatch timer; timer.Start(); gDebug=0; // Load basic libraries // If it does not work, please check the path of the libs and put it by hands gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); basiclibs(); // Load this example libraries gSystem->Load("libGeoBase"); gSystem->Load("libParBase"); gSystem->Load("libBase"); gSystem->Load("libMCStack"); gSystem->Load("libField"); gSystem->Load("libPassive"); gSystem->Load("libGen"); gSystem->Load("libgenfit"); gSystem->Load("libtrackrep"); gSystem->Load("libtpc"); gSystem->Load("libtpcreco"); gSystem->Load("librecotasks"); //gSystem->Load("libHypGe"); gSystem->Load("libDpmEvtGen"); gSystem->Load("libPGen");// add other detector's lib if you need them gSystem->Load("libHyp"); gSystem->Load("libTof"); gSystem->Load("libtpc"); //gSystem->Load("libHypGe"); CbmRunSim *fRun = new CbmRunSim(); // Set the number of events Int_t nEvents = 5;//50000;//---first step XXb //Int_t nEvents = 7396;//---second step hyperfragments //TString inFile= "/d/panda01/asanchez/PANDA/ximinsol1Ascii.dat"; //TString inFile= "/d/panda01/asanchez/PANDA/ximin5dmsol1Ascii.dat"; //--Second xi minus momentum solution. //TString inFile= "/d/panda01/asanchez/PANDA/ximinAscii.dat"; //--Second xi minus momentum solution only Ximin P<0.5GeV. //TString inFile= "ximinAsciiStpRate.dat"; //---root file hyperfragment TString inFile= "hypBupVDecay.root";//absorption veretx is included //TString inFile= "hypBupDecay2.root";//hypBupDecay3.root gam without corresp.frag //TString inFile= "../../hyp/xibximin2AStpRate.dat"; // ------ anti xi-minus annihilation products for second solution. TString parOutput="BupG3checkparams.root"; //TString inFile2= "xibar2Ascii_t1.dat"; // set the MC version used // ------------------------ fRun->SetName("TGeant3"); //fRun->SetOutputFile("Buphypcheck.root"); //-------second step, production of hyperfragments //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/BuphypAllXXb.root"); //-------first step, determ. vertex //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hypAllXXb.root"); fRun->SetOutputFile("sim_BupGhypcheck.root"); //fRun->SetOutputFile("simnewgeo.root"); //hypST_new.geo new local coordinates x-y //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/hypRes/sim_hypStp5.root"); //old local coordinates x-z, y-z. //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/hypRes/sim_hypStp4.root"); //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hyplaystp.root"); //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hyptessts.root"); //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hypPipe.root"); //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hypPipeVac.root"); //fRun->SetOutputFile("/d/panda01/asanchez/PANDA/sim_hypXB1.root"); // Set Material file Name //----------------------- fRun->SetMaterials("media_pnd.geo"); // Create and add detectors //------------------------- CbmModule *Cave= new PndCave("CAVE"); Cave->SetGeometryFileName("cave.geo"); fRun->AddModule(Cave); /* CbmModule *Magnet= new CbmMagnet("MAGNET"); Magnet->SetGeometryFileName("magnet.geo"); fRun->AddModule(Magnet); */ CbmDetector *Hyp = new PndHyp("HYP",kTRUE); //Hyp->SetGeometryFileName("HypST_newxy.geo"); Hyp->SetGeometryFileName("HypST_newxy3C.geo"); fRun->AddModule(Hyp); CbmDetector *Tof = new PndTof("TOF",kTRUE); Tof->SetGeometryFileName("tofSciF.geo"); fRun->AddModule(Tof); // CbmDetector *Tpc = new PndTpcDetector("TPC", kTRUE); Tpc->SetGeometryFileName("tpc.geo"); fRun->AddModule(Tpc); /* CbmDetector *HypGe= new PndHypGe("HYPGE",kTRUE); //HypGe->SetGeometryFileName("GeClusterPosfull.dat"); fRun->AddModule(HypGe); */ // Create and Set Event Generatorcczz //------------------------------- CbmPrimaryGenerator* primGen = new CbmPrimaryGenerator(); fRun->SetGenerator(primGen); fRun->SetUserDecay(kTRUE); /* CbmParticle *Ion_3_2 = new CbmParticle( "He3" ,2 , 3 , 0,2.80841,2,kTRUE,-1); fRun->AddNewParticle(Ion_3_2);*/ CbmParticle *Ion_6_3 = new CbmParticle( "Ion_6_3" ,3 , 6 , 0,5.6015,3,kTRUE,-1); fRun->AddNewParticle(Ion_6_3); CbmParticle *Ion_11_5 = new CbmParticle( "Ion_11_5" ,5 , 11 ,0,10.2525,5,kTRUE,-1); fRun->AddNewParticle(Ion_11_5); CbmParticle *Ion_7_3 = new CbmParticle( "Ion_7_3" ,3 , 7 , 0,6.5338,3,kTRUE,-1); fRun->AddNewParticle(Ion_7_3); CbmParticle *Ion_8_3 = new CbmParticle( "Ion_8_3" ,3 , 8 , 0,7.4714,3,kTRUE,-1); fRun->AddNewParticle(Ion_8_3); CbmParticle *Ion_6_4 = new CbmParticle( "Ion_6_4" ,4 , 6 , 0,5.60533,4,kFALSE,7.17e-12); fRun->AddNewParticle(Ion_6_4); CbmParticle *Ion_7_4 = new CbmParticle( "Ion_7_4" ,4 , 7 ,0,6.53423,4,kTRUE,-1); fRun->AddNewParticle(Ion_7_4); CbmParticle *Ion_9_4 = new CbmParticle( "Ion_9_4" ,4 , 9 , 0,8.392,4,kTRUE,-1); fRun->AddNewParticle(Ion_9_4); CbmParticle *Ion_8_5 = new CbmParticle( "Ion_8_5" ,5 , 8 ,0,7.472,5,kTRUE,-1); fRun->AddNewParticle(Ion_8_5); CbmParticle *Ion_10_6 = new CbmParticle( "Ion_10_6" ,6 , 10 ,0,9.327,6,kTRUE,-1); fRun->AddNewParticle(Ion_10_6); CbmParticle *Ion_11_6 = new CbmParticle( "Ion_11_6" ,6 , 11 ,0,10.254,6,kTRUE,-1); fRun->AddNewParticle(Ion_11_6); CbmParticle *Ion_12_6 = new CbmParticle( "Ion_12_6" ,6 , 12 , 0,11.174,6,kTRUE,-1); fRun->AddNewParticle(Ion_12_6); CbmParticle *Ion_12_7 = new CbmParticle( "Ion_12_7" ,7 , 12 ,0,11.191,7,kTRUE,-1); fRun->AddNewParticle(Ion_12_7); CbmParticle *Ion_13_7 = new CbmParticle( "Ion_13_7" ,7 , 13 ,0,12.111,7,kTRUE,-1); fRun->AddNewParticle(Ion_13_7); CbmParticle *Ion_1_1 = new CbmParticle( "Ion_1_1" ,1 , 1 , 0,0.9383,1,kTRUE,-1); fRun->AddNewParticle(Ion_1_1); CbmParticle *Ion_8_4 = new CbmParticle( "Ion_8_4" ,4 , 8 , 0,7.454,4,kTRUE,-1); fRun->AddNewParticle(Ion_8_4); CbmParticle *Ion_9_5 = new CbmParticle( "Ion_9_5" ,5 , 9 , 0,8.393,5,kFALSE,1.26e-11); fRun->AddNewParticle(Ion_9_5); CbmParticle *Be9L = new CbmParticle("Be9L", 4,9, 1,8.56449,4,kFALSE,0.200); fRun->AddNewParticle(Be9L); CbmParticle *Li9LL = new CbmParticle("Li9LL", 3,9, 2,8.74885,3,kFALSE,0.125);//kTRUE,-1); fRun->AddNewParticle(Li9LL); CbmParticle *he5 = new CbmParticle("He5", 2,5, 1,4.8399,2,kFALSE,0.200); fRun->AddNewParticle(he5); CbmParticle *H3L = new CbmParticle("H3L", 1,3, 1,2.99114,1,kFALSE,0.200); fRun->AddNewParticle(H3L); CbmParticle *H4L = new CbmParticle("H4L", 1,4, 1,3.9225,1,kFALSE,0.200); fRun->AddNewParticle(H4L); CbmParticle *He4L = new CbmParticle("He4L", 2,4, 1,3.92119,2,kTRUE,-1);//kFALSE,0.200); fRun->AddNewParticle(He4L); // CbmParticle *He5L = new CbmParticle("He5L", 2,5, 1,4.8399,2,kTRUE,-1); // fRun->AddNewParticle(He5L); CbmParticle *He6L = new CbmParticle("He6L", 2,6, 1,5.779,2,kFALSE,0.200); fRun->AddNewParticle(He6L); CbmParticle *He8L = new CbmParticle("He8L", 2,8, 1,7.6542,2,kFALSE,0.200); fRun->AddNewParticle(He8L); CbmParticle *Li6L = new CbmParticle("Li6L", 3,6, 1,5.7787,3,kFALSE,0.200); fRun->AddNewParticle(Li6L); CbmParticle *Li7L = new CbmParticle("Li7L", 3,7, 1,6.71153,3,kFALSE,0.200); fRun->AddNewParticle(Li7L); CbmParticle *Li8L = new CbmParticle("Li8L", 3,8, 1,7.6418,3,kFALSE,0.200); fRun->AddNewParticle(Li8L); CbmParticle *Li9L = new CbmParticle("Li9L", 3,9, 1,8.5787,3,kFALSE,0.200); fRun->AddNewParticle(Li9L); CbmParticle *Be7L = new CbmParticle("Be7L", 4,7, 1,6.7159,4,kFALSE,0.200); fRun->AddNewParticle(Be7L); CbmParticle *Be8L = new CbmParticle("Be8L", 4,8, 1,7.64315,4,kFALSE,0.200); fRun->AddNewParticle(Be8L); CbmParticle *B10L = new CbmParticle("B10L", 5,10, 1,9.49969,5,kFALSE,0.200); fRun->AddNewParticle(B10L); CbmParticle *B11L = new CbmParticle("B11L", 5,11, 1,10.41176,5,kFALSE,0.200); fRun->AddNewParticle(B11L); CbmParticle *B12L = new CbmParticle("B12L", 5,12, 1,11.35712,5,kFALSE,0.200); fRun->AddNewParticle(B12L); CbmParticle *C12L = new CbmParticle("C12L", 6,12, 1,11.35871,6,kFALSE,0.200); fRun->AddNewParticle(C12L); CbmParticle *C13L = new CbmParticle("C13L", 6,13, 1,12.27871,6,kFALSE,0.200); fRun->AddNewParticle(C13L); CbmParticle *He6LL = new CbmParticle("He6LL", 2,6, 2,5.95137,2,kFALSE,0.200); fRun->AddNewParticle(He6LL); CbmParticle *He7LL = new CbmParticle("He7LL", 2,7, 2,6.88982,2,kFALSE,0.200); fRun->AddNewParticle(He7LL); CbmParticle *Li7LL = new CbmParticle("Li7LL", 3,7, 2,6.88935,3,kFALSE,0.200); fRun->AddNewParticle(Li7LL); CbmParticle *Li8LL = new CbmParticle("Li8LL", 3,8, 2,7.82127,3,kFALSE,0.200); fRun->AddNewParticle(Li8LL); CbmParticle *Be9LL = new CbmParticle("Be9LL", 4,9, 2,8.75104,4,kFALSE,0.200); fRun->AddNewParticle(Be9LL); CbmParticle *Be10LL = new CbmParticle("Be10LL", 4,10, 2,9.67108,4,kFALSE,0.200); fRun->AddNewParticle(Be10LL); CbmParticle *Be11LL = new CbmParticle("Be11LL", 4,11, 2,10.60335,4,kFALSE,0.200); fRun->AddNewParticle(Be11LL); CbmParticle *Be12LL = new CbmParticle("Be12LL", 4,12, 2,11.53405,4,kFALSE,0.200); fRun->AddNewParticle(Be12LL); CbmParticle *B12LL = new CbmParticle("B12LL", 5,12, 2,11.53252,5,kFALSE,0.200); fRun->AddNewParticle(B12LL); CbmParticle *B11LL = new CbmParticle("B11LL", 5,11, 2,10.6037,5,kFALSE,0.200); fRun->AddNewParticle(B11LL); CbmParticle *B13LL = new CbmParticle("B13LL", 5,13, 2,12.45659,5,kFALSE,0.200); fRun->AddNewParticle(B13LL); /*CbmParticle *alpha = new CbmParticle("Alpha", 2,4,0,3.727417,2,kTRUE,-1); //lifetime was 1E20 fRun->AddNewParticle(alpha);*/ fRun->SetUserDecay(kTRUE); //Int_t pdg; // TDatabasePDG *db= TDatabasePDG::Instance(); // TParticlePDG *p=0; // p=db->GetParticle(geant3->GetIonPdg(1,2));//He5 // if(p) pdg = p->PdgCode(); // cout<<" Lambda-He4 "<< pdg <GetPDG()<< endl; // cout<<" Alpha geant4 "<GetPDG()<< endl; // Box Generator: /* PndBoxGenerator* boxGen = new PndBoxGenerator(-211, 1); // 13 = muon; 1 = multipl. // 211 = pi+ // first number: PDG particle code: 2nd number: particle multiplicity per event boxGen->SetPRange(0.7,0.7); // GeV/c // boxGen->SetPtRange(1.,1.); // GeV/c boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree] boxGen->SetThetaRange(40., 180.); // Polar angle in lab system range [degree] boxGen->SetCosTheta(); // Set uniform ditribution in cos(theta) boxGen->SetXYZ(0., 0., -76.); // vertex coordinates [cm] primGen->AddGenerator(boxGen); */ //CbmParticleGenerator* partGen = new CbmParticleGenerator(3312, 1, -0.8,0.8,0.8, 0., 0., -76); //primGen->AddGenerator(partGen); //*** with Ascii inFile ***first step //CbmAsciiGenerator* AsciiGen = new CbmAsciiGenerator(inFile); //primGen->AddGenerator(AsciiGen); //*** with root inFile ***second step PndHypBupGenerator* partGen = new PndHypBupGenerator(inFile); primGen->AddGenerator(partGen); // CbmAsciiGenerator* AsciiGen2 = new CbmAsciiGenerator(inFile2); // primGen->AddGenerator(AsciiGen2); fRun->SetStoreTraj(kTRUE); // to store particle trajectories //magnetic field: no field when commented put //CbmFieldConst *fMagField=new CbmFieldConst(); //fMagField->SetField(0.,0.,20.); // values are in kG //fMagField->SetFieldRegion(-50, 50,-50, 50, -100, 100); // values are in cm (xmin,xmax,ymin,ymax,zmin,zmax) //fRun->SetField(fMagField); PndMultiField *fField= new PndMultiField(); PndTransMap *map= new PndTransMap("TransMap", "R"); PndDipoleMap *map1= new PndDipoleMap("DipoleMap", "R"); PndSolenoidMap *map2= new PndSolenoidMap("SolenoidMap", "R"); fField->AddField(map); fField->AddField(map1); fField->AddField(map2); fRun->SetField(fField); /* PndConstField *fMagField=new PndConstField(); fMagField->SetField(0, 0 ,20. ); // values are in kG // MinX=-75, MinY=-40,MinZ=-12 ,MaxX=75, MaxY=40 ,MaxZ=124 ); // values are in cm fMagField->SetFieldRegion(-50, 50,-50, 50, -200, 200); fRun->SetField(fMagField);*/ /*CbmTrajFilter* trajFilter = CbmTrajFilter::Instance(); trajFilter->SetStepSizeCut(0.001); // 1 cm // trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.); // trajFilter->SetMomentumCutP(10e-3); // p_lab > 10 MeV // trajFilter->SetEnergyCut(0., 1.02); // 0 < Etot < 1.04 GeV trajFilter->SetStorePrimaries(kTRUE); trajFilter->SetStoreSecondaries(kTRUE);*/ // not used for the others.???? fRun->Init(); // Fill the Parameter containers for this run //------------------------------------------- CbmRuntimeDb *rtdb=fRun->GetRuntimeDb(); Bool_t kParameterMerged=kTRUE; CbmParRootFileIo* output=new CbmParRootFileIo(kParameterMerged); //output->open("simparams.root"); output->open(parOutput.Data()); rtdb->setOutput(output); rtdb->saveOutput(); rtdb->print(); // Transport nEvents // ----------------- fRun->Run(nEvents); timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime); }