// ---------------------------------------------------------------------------- // $Id: root_GoGmGmBkg.C,v 1.1 2001/12/07 19:22:05 toshi Exp $ // ---------------------------------------------------------------------------- // // Run Fast MC with overlayed gamma-gamma backgrounds, interactively with Root // // June 8 20001 T. Abe // // :: How to run :: // // 1) to call the root, type // // root // // 2) // root [0] .x root_GoGmGmBkg.C // // .. that's it. // // For example, if you like to run 300 events, type like // root [0] .x root_GoGmGmBkg.C(300) // Double_t partchg(); Int_t partchg3(); Int_t DoInit(); int root_GoGmGmBkg(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"); if (strcmp(gSystem->GetName(),"WinNT")){ gSystem->Load("libPythia6"); } gSystem->Load("libEG"); gSystem->Load("libEGPythia6"); LCDEvent event; //Create Physics event container first. LCDEvent event0; //Create Physics event container first. // Input root data Char_t* rootFile = "test.root"; LCDreadRootFile source(rootFile,&event0); // 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(); //-------------------------------------------------------------------------- // gamma-gamma background stuff. (generated by Pythia) //-------------------------------------------------------------------------- Double_t ecm = 500.; //center of mass energy. You have to choose same value of input MC events. TPythia6 gamgam; // interface to Pythia6.1 Event Generator gamgam.SetPMAS(6,1,175.); // set top mass. //...Sample program for minumum bias events by gamma*gamma*, //...internally convoluted with gamma* flux in e. //...WARNING: matching to DIS processes when one photon //...is very virtual is still missing in PYTHIA 6.143. //...Furthermore, so far the treatment of elastic and diffractive //...GVMD states defaults back to VMD ones, i.e. to lower masses //...than ought to be the case. Plus some other simplifications. //...Three (main) options of event classes: //...= 1 : jet events only (set pTmin in CKIN(3)). //...= 2 : jet plus low-pT events (pTmin irrelevant). //...= 3 : jet plus low-pT plus elastic and diffractive events. gamgam.SetCKIN(3,0.0); gamgam.SetMSEL(2); //...Set minimal and maximal W of gamma*gamma* system. gamgam.SetCKIN(77,5.0); gamgam.SetCKIN(78,ecm/2); //...Use minimal W to set minimum photon energy fractions. gamgam.SetCKIN(61,TMath::Power(gamgam.GetCKIN(77)/ecm,2)); gamgam.SetCKIN(63,gamgam.GetCKIN(61)); gamgam.SetCKIN(73,gamgam.GetCKIN(61)); gamgam.SetCKIN(75,gamgam.GetCKIN(61)); //...Set maximum photon virtuality. gamgam.SetCKIN(66,1.0); gamgam.SetCKIN(68,1.0); //...Initialize. gamgam.Pyinit("CMS","GAMMA/E+","GAMMA/E-",ecm); gamgam.Pystat(2); //-------------------------------------------------------------------------- // 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 Gamma-Gamma 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 iEvent; for (iEvent = 0; iEvent < nEvent; iEvent++) { if (!source.GetEvent()) break; // Read an event from the root file. Int_t nmc00=event0.MCparticles()->GetEntries(); event.Clean(); for (Int_t igen=0 ; igen < nmc00 ; igen++) { new((*event.MCparticles())[igen]) LCDMcPart((LCDMcPart*)(event0.MCparticles()->At(igen))); } // Poisson. Generate gamma-gamma backgounds per 0.2 event. Int_t ngen_gamgam=ran_gamgam.Poisson(0.2); NGmGm->Fill(ngen_gamgam); // Generate gamma-gamma background. for (Int_t igen_gamgam=0; igen_gamgam < ngen_gamgam; igen_gamgam++) { gamgam.GenerateEvent(); Double_t zoffset= ran_gamgam.Gaus(0.0,0.010);//Gaus mean=0 sigma=0.01cm ZGmGm->Fill(zoffset); Int_t nmc0 = event.MCparticles()->GetEntries(); Int_t nmcpy = gamgam.GetN(); for (Int_t imcpy=0 ; imcpy < nmcpy ; imcpy++) { Int_t statpy = gamgam.GetK(imcpy,1); Int_t idpy = gamgam.GetK(imcpy,2); Double_t chrgpy = partchg(idpy); LCDMcPart* part = new((*event.MCparticles())[nmc0+imcpy]) LCDMcPart(); part->SetStatus(statpy); part->SetParticleID(idpy); part->SetCharge(chrgpy); part->SetPosition(TVector3(gamgam.GetV(imcpy,1)/10.0, gamgam.GetV(imcpy,2)/10.0, gamgam.GetV(imcpy,3)/10.0+zoffset)); part->Set4Momentum(TLorentzVector(gamgam.GetP(imcpy,1), gamgam.GetP(imcpy,2), gamgam.GetP(imcpy,3), gamgam.GetP(imcpy,4))); // Fill parent and daughter information Int_t iprnt=gamgam.GetK(imcpy,3)-1+nmc0; part->SetParentIdx(iprnt); if (iprnt>=0) { LCDMcPart* prnt = (LCDMcPart*)(event->MCparticles()->At(iprnt)); prnt->SetTermPosition(part->GetPosition()); } } } Int_t nmcgmgm=event.MCparticles()->GetEntries() - nmc00 ; 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; } //----------------------------------------------------------------------------- // belows are copied from Util/src/partchg.c //----------------------------------------------------------------------------- /* Original table from stdhep file hepchg.F data ichg/-1,2,-1,2,-1,2,-1,2,2*0, 1 - 10 1 -3,0,-3,0,-3,0,-3,3*0, 11 - 20 2 3*0,3,6*0, 21 - 30 3 3*0,3,0,0,3,3*0, 31 - 40 4 10*0 41 - 50 5 0,0,6,3,6,5*0, 51 - 60 6 40*0,10*0/ */ /* Initialize table. This routine is called by PartChg3 first time only. */ static void DoInit(Int_t *tbl) { /* most of the entries are 0 */ Int_t nLeft = 110; Int_t *pEntry = tbl; while (nLeft--) { *pEntry++ = 0; } /* Now put in the non-zero entries */ tbl[1] = -1; tbl[2] = 2; tbl[3] = -1; tbl[4] = 2; tbl[5] = -1; tbl[6] = 2; tbl[7] = -1; tbl[8] = 2; tbl[11] = -3; tbl[13] = -3; tbl[15] = -3; tbl[17] = -3; tbl[24] = 3; tbl[34] = 3; tbl[37] = 3; tbl[53] = 6; tbl[54] = 3; tbl[55] = 6; return; } Int_t partchg3(Int_t id) { UInt_t kqa, kqn, kq1, kq2, kq3, kqj, irt, kqx; Int_t hepchg; // // To lessen chances of Fortran-to-c conversion problems, let "real" entries // start at index = 1 // Int_t ichg[110]; DoInit(&ichg[0]); /* Initial values. Simple case of direct readout. */ hepchg = 0; kqa = abs(id); kqn = kqa/1000000000 % 10; kqx = kqa/1000000 % 10; kq3 = kqa/1000 % 10; kq2 = kqa/100 % 10; kq1 = kqa/10 % 10; kqj = kqa % 10; irt = kqa % 10000; /* illegal or ion */ if ((kqa == 0) || (kqa >= 10000000)) { /* set ion charge to zero - not enough information */ if (kqn == 1) hepchg = 0; } /* direct translation; i.e., can use abs(id) as index into table */ else if (kqa <= 100) { hepchg = ichg[kqa]; } /* KS and KL (and undefined) */ else if (kqj == 0) { hepchg = 0; } /* direct translation -- use irt as index into table */ else if ((kqx > 0) && (irt <= 100)) { hepchg = ichg[irt]; } /* Construction from quark content for heavy meson, diquark, baryon. */ else if (kq3 == 0) { /* Mesons. */ hepchg = ichg[kq2] - ichg[kq1]; /* Strange or beauty mesons. */ if ((kq2 == 3) || (kq2 == 5)) hepchg = ichg[kq1] - ichg[kq2]; } else if (kq1 == 0) { /* Diquarks. */ hepchg = ichg[kq3] + ichg[kq2]; } else { /* Baryons */ hepchg = ichg[kq3] + ichg[kq2] + ichg[kq1]; } /* fix sign of charge */ if ((id < 0) && (hepchg != 0)) hepchg = -hepchg; return hepchg; } Double_t partchg(Int_t id) { return (partchg3(id) / 3.0); }