// ---------------------------------------------------------------------------- // $Id: pandpyth_GoLCDVTopl.C,v 1.6 2001/06/28 17:38:42 toshi Exp $ // ---------------------------------------------------------------------------- // $Log: pandpyth_GoLCDVTopl.C,v $ // Revision 1.6 2001/06/28 17:38:42 toshi // Changes for windows. // // Revision 1.5 2001/06/27 22:58:04 masako // Use LCDPandPyth2MCPart instead of LCDPyj2McPart. // // Revision 1.4 2001/06/19 16:51:57 toshi // Use gSystem->Getenv() to get Parfiles directory. // // Revision 1.3 2001/05/11 21:52:29 toshi // Remove .so suffix from gSystem->Load("...") in order to call macros from // Windows system. // // Revision 1.2 2001/05/09 21:56:28 toshi // Add CVS tag $Id: pandpyth_GoLCDVTopl.C,v 1.6 2001/06/28 17:38:42 toshi Exp $. // // // Run Pandora_Pythia and Fast MC interactively with Root // Sample program using LCDVTopl (translated from SLD ZVTOP) // 04/26/2001 Toshi // // :: How to run :: // // 1) to call the root, type // // root // // 2) // root [0] .x pandpyth_GoLCDVTopl.C // // .. that's it. // // For example, if you like to run 300 events, type like // root [0] .x pandpyth_GoLCDVTopl.C(300) // Int_t pandpyth_GoLCDVTopl(Int_t nEvent=300){ gROOT->Reset(); TString parfile_dir; parfile_dir += gSystem->Getenv("LCDROOT") ;parfile_dir += "/"; parfile_dir += gSystem->Getenv("LCDVERSION");parfile_dir += "/"; parfile_dir += "ParFiles/"; // In order to refer to shareable libraries in this simple // form, include their paths in a .rootrc file. gSystem->Load("libPandora"); gSystem->Load("libPandoraProc"); if (strcmp(gSystem->GetName(),"WinNT")) { gSystem->Load("libPythia6"); } gSystem->Load("libEG"); gSystem->Load("libEGPythia6"); if (strcmp(gSystem->GetName(),"WinNT")) { gSystem->Load("libTauola"); } gSystem->Load("libPandPyth"); gSystem->Load("libLCDEvent"); gSystem->Load("libLCDRootAppsUtil"); gSystem->Load("libLCDFastMC"); gSystem->Load("libLCDPhUtil"); gSystem->Load("libLCDGenUtil"); //------------< Setup Generator (Pandora-Pythia) >------------------------ printf("Start Generator Setup!\n"); #include "pandpythconsts.h" Int_t iseed_pan = 7; // default pandora initial seed /* define the pandora event selection in this space */ double ECM = 91.26; double PolE = 0.0; // Polarization Electron double PolP = 0.0; // Polarization Positron ebeam b1(ECM/2.0, PolE, electron, electron); b1.setup(NLC500); b1.ISRoff(); // Turn off initial-state radiation b1.beamstrahlungoff(); // Turn off beamstrahlung ebeam b2(ECM/2.0, PolP, positron, positron); b2.setup(NLC500); b2.ISRoff(); b2.beamstrahlungoff(); double ECMmin = 20.; // Higgs Mass eetoqqbar pr(ECMmin); //For example:: Z0 decays to anything //const Int_t invisibleOnly = 557; //const Int_t bOnly = 553; //const Int_t cOnly = 552; //const Int_t allStates = 546; //pr.onlyDecay(allStates); pandora P(b1,b2,pr); pandorarun PR(&P, epluseminus, ECM, nEvent, iseed_pan); //----- For PYTHIA ------ //Or you can also change PYTHIA seed here. int iseed_pyth = 22846; PR.newseed(iseed_pyth); // Change Seed //PR.partonshoweroff(); // Turn off parton shower //PR.fragmentationoff(); // Turn off fragmentation //PR.FSRoff(); // Turn off final-state QCD QED radiation /* end of definition */ PR.initialize(0); /* // If you want/don't want to decay some specific particle(s), // you can call PartDecayOn/PartDecayOff here. // For example: Do not allow to decay K0S and Lambda int ipart1 = 310; // K0S int ipart2 = 3122; // Lambda PR.PartDecayOff(ipart1); PR.PartDecayOff(ipart2); */ printf("End of Generator Setup!\n"); //Pyjet -> McPart converter //LCDPyj2McPart py2mc; LCDPandPyth2McPart py2mc; //------------< End Generator (Pandora-Pythia) setup>------------------- //Physics info.. LCDEvent event; //---------------- //---> Detector Geometry char* geomFile = "Silicon.par"; // Silicon TString s_geomFile(parfile_dir); s_geomFile += geomFile; FILE* m_geomfile = fopen(s_geomFile.Data(),"r"); LCDGetParameters gp(m_geomfile); fclose(m_geomfile); // If you want to change some detector parameter(s) for FastMC // do here // For example.. //gp.SetMagneticField(2.); // Change Magnetif field //gp.SetCoilinHAD(); // Set HAD inside the Coil (This case // you need to change some detector size also..) //---> Tracking smearing parameter file Char_t* smearFile = "lookup_silicon.nobmcon"; TString s_smearFile(parfile_dir); s_smearFile += smearFile; //---> FastMC LCDFastMC fmc(&gp,s_smearFile.Data(),&event); // Define histograms TH1F* h_nvtvrt = new TH1F("h_nvrt","N of LCDVToplVRT" ,10,0.0,10.0); TH1F* h_rwm = new TH1F("h_rwm" ,"Raw Vrt Mass" ,40,0.0,10.0); TH1F* h_ptm1 = new TH1F("h_ptm1" ,"Pt corr. Vertex Mass" ,28,0.0,10.0); TH1F* h_ptm2 = new TH1F("h_ptm2" ,"Pt corr. Vertex Mass" ,28,0.0,10.0); TH2F* h_mvsp = new TH2F("h_mvsp" ,"Pt Mass vs Vertex P uds" , 20,0.0,2.0,25,0.0,50.); h_nvrt ->SetFillColor(4); h_rwm ->SetFillColor(4); h_ptm1 ->SetFillColor(4); h_ptm2 ->SetFillColor(4); h_mvsp ->SetFillColor(4); TCanvas *c1 = new TCanvas("c1", "c1", 600,600); c1->Divide(2,2); TCanvas *c2 = new TCanvas("c2", "c2", 600,600); c2->Divide(2,2); // Setup of JetFinder LCDJetFinder jetfind(0.004); TClonesArray forjet("TLorentzVector",1000); // Track cut value Double_t cosCut=0.9, pCut=0.2; // Event loop LCDVTopl topl; Int_t iEvent; for (iEvent = 0; iEvent < nEvent; iEvent++) { PR.getevent();//Event Generate event.MCparticles()->Delete();//prepare for making McPart //py2mc.MakeMcPart(event.MCparticles());//Make McPart py2mc.MakeMcPart(&PR,event.MCparticles());//Make McPart fmc.Doit();//do fast detector simulation... // Look at tracks. Can access them by the [] operator. Int_t nTracks = event->Tracks()->GetEntries(); // Cut 1... if (nTracks < 5) goto canvasdraw; // Check number of good tracks TObjArray gtrk_list; LCDTrack* trk=0; TVector3 ptrk; Int_t itrk=0; for (itrk=0; itrkTracks()->At(itrk)); ptrk=trk->GetMomentumVector(0.0); if (TMath::Abs(ptrk.CosTheta())pCut) { gtrk_list.Add((TObject*)trk); } } //Cut 2... Int_t nGTracks=gtrk_list.GetEntries(); if (nGTracks < 5 ) goto canvasdraw; TVector3 vec3; // Do Jet Finding forjet.Delete(); // Clean up veclist for (itrk=0; itrkGetMomentumVector(0.0); Double_t Etrk = vec3.Mag2() + 0.1396 * 0.1396;// Assume pion mass new((forjet)[itrk]) TLorentzVector(vec3,Etrk); } jetfind.SetPartList(&forjet); // Input particle 4-vector list jetfind.SetDURHAM(); // Use DURHAM algorithm jetfind.SetYCut(0.004); // Set YCut value jetfind.doFindJets(); // JetFindig TVector3 ipv=((LCDBeam*)(event->BEAM()))->GetPos(); // IP position //TVector3 ipv(0.,0.,0.); TVector3 bfield(0.0,0.0,gp.GetMagneticField()); // B field TLorentzVector p4vrt; Double_t rawMass, ptMass, massTag; TLorentzVector jet4vec; TVector3 jet3vec; Int_t nJet = jetfind.njets(); //Loop Over # of jets Int_t* trkList=new Int_t[nGTracks]; for (Int_t iJet=0 ; iJet < nJet ; iJet++) { Int_t nTrks = jetfind.nParticlesPerJet(iJet); if (nTrks <= 2) continue; Int_t iTrk2=0; for (int ipart =0; ipartFill(nVRT); if (nVRT <= 1) continue; p4vrt=topl.GetP4Vertex(); rawMass = topl.GetRawMass(); if (rawMass <= 0.) continue; ptMass = topl.GetPtMass(); massTag = topl.GetMassTag(); h_rwm ->Fill(rawMass); h_ptm1->Fill(ptMass); h_ptm2->Fill(massTag); if (massTag<2.) h_mvsp->Fill(massTag,p4vrt.Rho()); } delete trkList; canvasdraw: if (iEvent%100 == 0 || iEvent == nEvent-1) { c1->cd(1); h_nvrt->Draw("B"); c1->cd(2); h_mvsp->Draw("BOX"); c1->Modified();c1->Update(); c1->cd(2); h_mvsp->Draw("BOX"); c1->Modified();c1->Update(); c2->cd(1); h_rwm ->Draw("B"); c2->cd(2); h_ptm1 ->Draw("B"); c2->cd(3); h_ptm2 ->Draw("B"); c2->Modified();c2->Update(); } if (iEvent%1000 == 999 || iEvent == nEvent-1) { printf("iEvent = %i was done...\n", iEvent+1); } } forjet.Delete(); topl.Clean(); event.Clean(); return iEvent; }