// ---------------------------------------------------------------------------- // $Id: stdhep_GoOverlayEvent.C,v 1.2 2001/09/20 01:47:42 toshi Exp $ // ---------------------------------------------------------------------------- // // Run Fast MC with overlayed events, interactively with Root // // :: How to run :: // // 1) to call the root, type // // root // // 2) // root [0] .x stdhep_GoOverlayEvent.C // // .. that's it. // // For example, if you like to run 300 events, type like // root [0] .x stdhep_GoOverlayEvent.C(300) // int stdhep_GoOverlayEvent(int nEvent=100) { gROOT->Reset(); TString parfile_dir; parfile_dir += gSystem->Getenv("LCDROOT") ;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("libLCDEvent"); gSystem->Load("libLCDRootAppsUtil"); gSystem->Load("libLCDFastMC"); gSystem->Load("libLCDStdHEPUtil"); LCDEvent event; //Create Physics event container first. // Input stdhep data Char_t* stdFile = "test.dat"; LCDreadStdFile source(stdFile,&event); // StdHep file index=0 printf("LCDreadStdFile initialization done\n"); // Overlayed event data Char_t* stdFile2= "test2.dat"; // StdHep file name for overlayed event source.AddStdFile(stdFile2); // StdHep file index=1 printf("LCDreadStdFile::AddStdFile initialization done\n"); // Open detector parameter file Char_t* geomFile = "SDMar01.par"; // SD //"LDMar01.par"; // LD 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"; //"lookup_large.nobmcon"; TString s_smearFile(parfile_dir); s_smearFile += smearFile; // Initialize FastMC LCDFastMC fmc(&gp,s_smearFile.Data(),&event); // If you don't want the Cluster merging //fmc.SetNoClusterMerge(); // Define a histogram or two TH1F *ClsMult = new TH1F("ClsMult","Cluster Multiplicity",50, 0.0, 200.0); TH1F *ClsE = new TH1F("ClsE","Cluster Energy",50, 0.0, 50.0); TH1F *TrkMult = new TH1F("TrkMult","Track Multiplicity",50, 0.0, 100.0); TH1F *TrkE = new TH1F("TrkE","Track Energy",50, 0.0, 50.0); TH1F *NGmGm = new TH1F("NGmGm","N of overlayed backgrounds",10,0.0,10.0); TH1F *ZGmGm = new TH1F("ZGmGm","Z offset",100,-0.1,+0.1); TH1F *MGmGm = new TH1F("MGmGm","N of bkg MC parts",100,0.0,100.0); Int_t iseed=7; TRandom ran_gamgam(iseed); // random number generator // Event loop Int_t nmcv[100]; //buffer to save # of events per read. Int_t iEvent; for (iEvent = 0; iEvent < nEvent; iEvent++) { if (!(nmcv[0]=source.ReadEvent(0))) break; // Read an event. // Make LCDMcPart source.MakeMcPart(); // Poisson. read overlayed events per 0.2 event. Int_t ngen_gamgam=ran_gamgam.Poisson(0.2); NGmGm->Fill(ngen_gamgam); // Read overlayed event. for (Int_t igen_gamgam=0; igen_gamgam < ngen_gamgam; igen_gamgam++) { nmcv[igen_gamgam+1]=source.ReadEvent(1); //1 means overlayed evt. file. // Make LCDMcPart source.MakeMcPart(0,1);//1 means overlay readevent printf("iEvent=%d\n",iEvent); } // Add zoffset, if you like. Int_t nmc=nmcv[0]; for (Int_t iov=0; iov < ngen_gamgam; iov++) { Double_t zoffset= ran_gamgam.Gaus(0.0,0.010);//Gaus mean=0 sigma=0.01cm ZGmGm->Fill(zoffset); for (Int_t imc=nmc ; imc < nmc+nmcv[iov+1] ; imc++) { LCDMcPart* part = (LCDMcPart*)event.MCparticles()->At(imc); Double_t newzi=part->GetPositionPtr()->Z() + zoffset; Double_t newzt=part->GetTermPositionPtr()->Z() + zoffset; part->GetPositionPtr()->SetZ(newzi); part->GetTermPositionPtr()->SetZ(newzt); } nmc += nmcv[iov+1]; } Int_t nmcgmgm=event.MCparticles()->GetEntries() - nmcv[0] ; MGmGm->Fill(nmcgmgm); fmc.Doit(); // Through FastMC TLorentzVector part_mom(0.,0.,0.,0.); // TObjArray* McPart = event.MCparticles(); Int_t nMCpart = McPart->GetEntries(); for (Int_t Imc=1 ; Imc < nMCpart ; Imc++) { LCDMcPart *mcpart = (LCDMcPart *) McPart->At(Imc); Int_t part_st = mcpart->GetStatus(); Int_t part_id = mcpart->GetParticleID(); part_mom = mcpart->Get4Momentum(); /* printf("%i ID = %i status = %i momentum ( %f, %f, %f ) E %f \n", Imc, part_id, part_st, part_mom.Px(), part_mom.Py(),part_mom.Pz(),part_mom.E()); */ } // Look at Clusters. Can access them by the [] operator. int nClusters = event.ClusterLst()->GetEntries(); ClsMult->Fill(nClusters); //printf("nClusters %i\n", nClusters); for (Int_t iClus = 0; iClus < nClusters; iClus++) { LCDCluster* cls = (LCDCluster*)event.ClusterLst()->At(iClus); if (cls) { Double_t eClus = cls->GetEnergy(); // Histogram energies //printf("energy = %f\n", eClus); ClsE->Fill(eClus); } } // Look at tracks. Can access them by the [] operator. int nTracks = event.Tracks()->GetEntries(); TrkMult->Fill(nTracks); // Track loop for (Int_t ntrk=0; ntrkAt(ntrk); if(trk) { Double_t* tP3 = trk->GetMomentum(); Double_t E_trk = 0.1396 * 0.1396; // Assume pion mass for (int i=0; i<3; i++) { E_trk += (*(tP3+i))*(*(tP3+i)); } E_trk = sqrt(E_trk); // Histogram energies TrkE->Fill(E_trk); } } } ClsMult->SetFillColor(4); ClsE->SetFillColor(4); TrkMult->SetFillColor(4); TrkE->SetFillColor(4); TCanvas *c1 = new TCanvas("c1", "c1", 600, 600); c1->Divide(2,2); c1->cd(1); ClsMult->Draw("B"); c1->cd(2); ClsE->Draw("B"); c1->cd(3); TrkMult->Draw("B"); c1->cd(4); TrkE->Draw("B"); NGmGm->SetFillColor(4); ZGmGm->SetFillColor(4); MGmGm->SetFillColor(4); TCanvas *c2 = new TCanvas("c2", "c2", 600, 600); c2->Divide(2,2); c2->cd(1); NGmGm->Draw("B"); c2->cd(2); ZGmGm->Draw("B"); c2->cd(3); MGmGm->Draw("B"); return iEvent; }