hadr_examp(){ //================================================================ // Toy Monte-Carlo to simulate hadronic cascade in calorimeter // taking into account energy conservation law // for all components of hadronic shower // V.L. Morgunov 9-March-2007 //================================================================ Float_t enr = 10.0; TH1F *h1 = new TH1F("vis","Visible",500,0.0,0.25*enr); TH1F *h2 = new TH1F("abs","Absorber",500,0.0,1.5*enr); TH1F *h3 = new TH1F("bind","Binding",500,0.0,1.5*enr); TH1F *h4 = new TH1F("vis_fact","Visible*fact",300,0.0,1.5*enr); TH1F *h5 = new TH1F("had","Hadronic",500,0.0,1.5*enr); TH1F *h6 = new TH1F("em","Electromagnetic",500,0.0,1.5*enr); TH2F *h7 = new TH2F("em_had_fact","Hadr vs EM fact", 100,0.0,1.5*enr,100,0.0,1.5*enr); h7->GetXaxis()->SetTitle("Ef "); h7->GetYaxis()->SetTitle("Hf"); TH2F *h8 = new TH2F("em_had","Hadr vs EM", 100,0.0,1.5*enr,100,0.0,1.5*enr); h8->GetXaxis()->SetTitle("EM "); h8->GetYaxis()->SetTitle("HAD"); //================================================================ // Factor to convert from visible to physical energy scale Float_t factor = 29.3; // Another factor, if binding losses can be measured // Float_t factor = 24.0; // One tile crossing by one particles = 836 keV Float_t mip_cost = 836.e-6; Float_t h_e_ratio = 0.75; Float_t binding_fraction = 0.25; //================================================================ for (int j=0;j<10000;j++){ // Event loop // Float_t rh = gRandom->Rndm(); //Uniform distribution Float_t rh = -1.0; while(!(rh > 0.0 && rh < 1.0)) // Gaussian-like distribution rh = gRandom->Gaus(0.5,0.3); // Hadron fraction in cascade Float_t re = 1.0-rh; // EM fraction in cascade // Binding energy is proportional to hadron fraction Float_t bind_g = binding_fraction*enr*rh; Float_t rest_bind = -1.0; while(rest_bind<=0.0) rest_bind = gRandom->Gaus(bind_g,0.2*bind_g); // but spreaded Float_t rest_he = enr - rest_bind; // All other goes to Had+EM Float_t rest_had = rest_he*rh; // Hadr still proportional itself Float_t rest_em = rest_he - rest_had; // EM fraction of energy // Visible parts of energy are about in MIPs Int_t nmips_had = rest_had/30./mip_cost; Int_t nmips_em = rest_em/30./mip_cost; Int_t nmips_bind = rest_bind/30./mip_cost; Float_t vis_fe = 0.0; Float_t vis_fh = 0.0; Float_t vis_fb = 0.0; // Randomize visible energy in pseudo-MIPs for(int i = 0;iGaus(2.0,3.0); for(int i = 0;iLandau(1.0,0.07); for(int i = 0;iGaus(1.2,1.2); // Convert back to absolute energy from pseudo--MIPs vis_fh = h_e_ratio*vis_fh*mip_cost; vis_fe = vis_fe*mip_cost; vis_fb = vis_fb*mip_cost; Float_t mrest_bind = rest_bind - vis_fb; // subtract visible from binding Float_t mrest_had = rest_had - vis_fh; // subtarct visible hard from hard Float_t mrest_em = rest_em - vis_fe; // subtarct visible em from em Float_t e_absorber = mrest_had + mrest_em; // If we can measure binding signal by neutrons -- whole sum //Float_t e_visible = fe+fh+fb; //Float_t factor = 24.0; // If we can NOT measure binding signal by neutrons -- partial sum Float_t e_visible = vis_fe+vis_fh; Float_t efact = e_visible*factor; Float_t ee = vis_fe*factor; Float_t eh = vis_fh*factor; //================================================================ h5->Fill(rest_had,1.0); h6->Fill(rest_em,1.0); h8->Fill(mrest_em,mrest_had,1.0); h7->Fill(ee,eh,1.0); h1->Fill(e_visible,1.0); h2->Fill(e_absorber,1.0); h3->Fill(mrest_bind,1.0); h4->Fill(efact,1.0); } // end of event loop //================================================================ gStyle->SetPadGridX(true); gStyle->SetPadGridY(true); gStyle->SetOptFit(1); gStyle->SetOptStat(0); gStyle->SetPalette(1); //-------------------------------------------------------- TCanvas* c1 = new TCanvas("c1","1",300,300); gPad->SetLogy(); h1->SetLineColor(107); h1->SetLineWidth(2); h1->Draw(""); //-------------------------------------------------------- TCanvas* c2 = new TCanvas("c2","2",300,300); gPad->SetLogy(); h2->SetLineColor(4); h2->SetLineWidth(2); h2->Draw(""); //-------------------------------------------------------- TCanvas* c3 = new TCanvas("c3","3",300,300); gPad->SetLogy(); h3->SetLineColor(1); h3->SetLineWidth(2); h3->Draw(""); //-------------------------------------------------------- TCanvas* c4 = new TCanvas("c4","4",300,300); // gPad->SetLogy(); h4->SetLineColor(2); h4->SetLineWidth(2); h4->Draw(""); //-------------------------------------------------------- // TCanvas* c5= new TCanvas("c5","5",300,300); // gPad->SetLogy(); // h5->SetLineColor(104); // h5->SetLineWidth(2); // h5->Draw(""); //-------------------------------------------------------- TCanvas* c6= new TCanvas("c6","6",300,300); // gPad->SetLogy(); h6->SetLineColor(2); h6->SetLineWidth(2); h6->Draw(""); h5->SetLineColor(104); h5->SetLineWidth(2); h5->Draw("Same"); //-------------------------------------------------------- TCanvas* ca = new TCanvas("ca","All",800,600); // gPad->SetLogy(); h4->SetMaximum(300.); h4->Fit("gaus"); h1->Draw("SAME"); h2->Draw("SAME"); h3->Draw("SAME"); //-------------------------------------------------------- TCanvas* c7 = new TCanvas("c7","HAD vs EM",300,300); gStyle->SetPadGridX(true); gStyle->SetPadGridY(true); gPad->SetLogz(); h8->Draw("colz"); lm1 = new TF1("lm1","10.0-x",0.,10.); lm1->Draw("SAME"); //-------------------------------------------------------- TCanvas* c8 = new TCanvas("c8","HAD vs EM fact",300,300); gStyle->SetPadGridX(true); gStyle->SetPadGridY(true); gPad->SetLogz(); h7->Draw("colz"); // lm1 = new TF1("lm1","10.0-x",0.,10.); lm1->Draw("SAME"); //-------------------------------------------------------- TFile *f = new TFile("hadr_exampl.root","RECREATE"); //-------------------------------------------------------- h1->Write(); h2->Write(); h3->Write(); h4->Write(); h5->Write(); h6->Write(); h7->Write(); h8->Write(); f->Close(); }