Calorimeter/FastHadronMC/fasthadv4.C

0001 #include "fasthadv4.h"
0002
0003 /*
0004   First shot at a fast hadron simulation tool.
0005   The "detector" is rectangular, nbin*nbin granularity, ddepth deep.
0006   Dimensions are arbitrary (make sense only in context with the
0007   electromagnetic and hadronic  shower width parameters dwidthem, dwidthhad.
0008   Impact point and theta, phi angles are arbitrary.
0009   Each hadron enters the detector as minimum ionizing; afetr each (fixed)
0010   MIP deposit throws a random number to decide whether to start a shower or not.
0011   Once it decided to form a shower, the remaining length in the detector is
0012   calculated, and serves as a weight for the energy it will deposit (the deposit
0013   in any single step is constant now, so one plays with the number of steps to change
0014   the total deposit).
0015   Also, the total energy it will actually deposited is weighted (inversely) with
0016   a formula, essentially with the log of the original hadron energy.
0017   The total energy it will deposit is "dshe".
0018   Next, it will be decided how much of "dshe" will be electromagnetic and
0019   how much is hadronic ("reh").    Since the electromagnetic part usually starts
0020   later, a spatial delay is assigned to the electromagnetic shower.
0021   Finally the two showers (em, had) are generated in three dimensions,
0022   gradually widening with depth, but randomized with a Gaussian.
0023  */
0024
0025 // Impact point and angle
0026 double dimpx,dimpy,dimpxtheta,dimpyphi;
0027 double dpath,dfullpath,ddepth,dtote,dtotemip,dtoteem,dtotehad;
0028 double x,y,z,xhad,yhad,zhad,reh;
0029 double dshe,dhade,deme;
0030 double ddzhad,ddzem;  // stepsize for had, em part (mip is fixed)
0031 double dfudge;
0032 double dcostheta,dcosphi,dsintheta,dsinphi;
0033 double demdelay;
0034
0035
0036 int nmipstep,nshstep;
0037 int i,j;
0038
0039
0040 int ntrysh,ntryem,ntryhad;
0041
0042
0043   TH2D *h2dhits = new TH2D("h2dhits","Hit pattern",nbin,xlo,xhi,nbin,ylo,yhi);
0044   TH2D *h2dlive = new TH2D("h2dlive","Live display",nbin,xlo,xhi,nbin,ylo,yhi);
0045   TH1D *htote = new TH1D("htote","Total energy",100,0.0,15.0);
0046   TH1D *htotemip = new TH1D("htotemip","Total energy",25,0.0,1.0);
0047   TH1D *htoteem = new TH1D("htoteem","Total energy",160,0.0,8.0);
0048   TH1D *htotehad = new TH1D("htotehad","Total energy",160,0.0,8.0);
0049   TH2D *h2dmipnonmip = new TH2D("h2dmipnonmip","MIP vs non MIP energy",
0050                 80,0.0,8.0,50,0.0,0.5);
0051   TH3D *h3dlivemip =  new TH3D("h3dlivemip","3D deposit",
0052                   nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
0053   TH3D *h3dliveem =  new TH3D("h3dliveem","3D em deposit",
0054                   nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
0055   TH3D *h3dlivehad =  new TH3D("h3dlivehad","3D had deposit",
0056                   nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
0057   TH3D *h3dliveall =  new TH3D("h3dliveall","3D all deposit",
0058                   nbin,xlo,xhi,nbin,ylo,yhi,nbin,zlo,zhi);
0059
0060
0061
0062 void fasthadv4(double dimpx=0.5, double dimpy=0.5,
0063            double dimpxtheta=0.0, double dimpyphi=0.0)
0064 {
0065
0066
0067   nmipstep = ndefmipstep;
0068   nshstep = ndefshstep;
0069
0070   dcostheta =  TMath::Cos(dimpxtheta);
0071   dcosphi =  TMath::Cos(dimpyphi);
0072   dsintheta =  TMath::Sin(dimpxtheta);
0073   dsinphi =  TMath::Sin(dimpyphi);
0074   if(dcostheta<=0.0 || dcosphi<=0.0)
0075     {
0076       std::cout << "Something dead wrong with impact angle " << std::endl;
0077       return;
0078     }
0079
0080   nmipstep = int(float(ndefmipstep)/(dcostheta*dcosphi));
0081   nshstep = int(float(ndefshstep)/(dcostheta*dcosphi));
0082   dfullpath = dlength/(dcostheta*dcosphi);
0083   std::cout << nmipstep << "    "  << nshstep << "   "  <<  dfullpath << std::endl;
0084
0085
0086
0087   dhade = dtothade / float(nshstep);
0088   deme = dtothade / float(nshstep);
0089   ddzem = dlength / float(ndefshstep);
0090   ddzhad = dlength / float(ndefshstep);
0091
0092 TCanvas  * c10 = new TCanvas("c10","Live display",1200,600);
0093   c10->Divide(2,1);
0094
0095
0096   for(i=1; i<=ntry;i++)
0097     {
0098
0099       dpath = ddepth = dtote = dtotemip = dtoteem = dtotehad = dshe = 0.0;
0100       h3dlivemip->Reset();
0101       h3dliveem->Reset();
0102       h3dlivehad->Reset();
0103       h2dlive->Reset();
0104
0105       for(j=1; j<=nmipstep; j++)
0106     {
0107       dpath = j*ddzmip;
0108       x = dimpx + dpath*dsintheta;
0109       y = dimpy + dpath*dsintheta;
0110       z = dpath * dcostheta*dcosphi;
0111       dtotemip  +=ddedx;
0112       h2dhits->Fill(x,y,ddedx);
0113       //      i%nlive == 0 ? h2dlive->Fill(x,y,ddedx) : ;
0114       if(i%nlive == 0)
0115         {
0116           h2dlive->Fill(x,y,ddedx);
0117           h3dlivemip->Fill(x,y,z);
0118         }
0119
0120       if(gRandom->Uniform()>dmipstop) break;
0121     }
0122       dtotemip *= 1.0 + dresol * gRandom->Gaus();
0123       htotemip->Fill(dtotemip);
0124       xhad = x;  yhad = y;
0125       zhad = dpath * dcostheta*dcosphi;
0126       ddepth = dpath * dcostheta*dcosphi;
0127       if(ddepth>=dlength)
0128     {
0129       htote->Fill(dtotemip);
0130       continue;  // done with this, generate the next particle
0131     }
0132
0133       // Get the remaining energy
0134       dshe = dtothade - dtotemip;
0135       // Estimate how much of it will be deposited (function of dtothade)
0136       dshe *= TMath::Sqrt(1.0/(1.0+0.5*TMath::Log(dshe+1.0)));
0137
0138       // Right now for simplicity the deposits per step are identical,
0139       // so you manipulate the number of steps to get the proper deposit
0140       ntrysh = int(nshstep * dshe/dtothade);
0141
0142       // Ratio of electromagnetic and hadronic part of the shower
0143       reh = 0.1 + 0.15*gRandom->Uniform();
0144       reh /=(1.0+0.6*TMath::Log(dshe+1.0));
0145
0146       ntryem = int(ntrysh * reh);
0147       ntryhad = int(ntrysh * (1.0 - reh));
0148       ntryem = int(ntryem * ((dlength-ddepth) / dlength));
0149       ntryhad = int(ntryhad * ((dlength-ddepth) / dlength));
0150
0151       // The em part of the hadron shower doesn't start immediately
0152       // it is
0153       demdelay = 2.0 + gRandom->Gaus();
0154       demdelay > 0.0 ?  : demdelay = -demdelay;
0155       if(ntryem>0) ddzem = (dlength-ddepth)/float(ntryem);
0156
0157
0158       dpath = 0.0;
0159       for(j=1; j<ntryem; j++)
0160     {
0161       dpath = j*ddzem + demdelay;
0162       if((ddepth + dpath * dcostheta*dcosphi)>dlength)
0163         {
0164           htote->Fill(dtoteem);
0165           break;
0166         }
0167       x = xhad + dpath*(dsintheta+dwidthem*gRandom->Gaus());
0168       y = yhad + dpath*(dsinphi+dwidthem*gRandom->Gaus());
0169       z = zhad+dpath * dcostheta*dcosphi;
0170       dtoteem += deme;
0171       h2dhits->Fill(x,y,deme);
0172       if(i%nlive == 0)
0173         {
0174           h2dlive->Fill(x,y,deme);
0175           h3dliveem->Fill(x,y,z);
0176         }
0177     }
0178       dtoteem *= 1.0 + dresol * gRandom->Gaus();
0179       htoteem->Fill(dtoteem);
0180
0181       if(ntryhad>0) ddzhad = (dlength-ddepth)/float(ntryhad);
0182
0183
0184       dpath = 0.0;
0185       for(j=1; j<ntryhad; j++)
0186     {
0187       dpath = j*ddzhad;
0188       if((ddepth + dpath * dcostheta*dcosphi)>dlength)
0189         {
0190           htote->Fill(dtotehad);
0191           break;
0192         }
0193       x = xhad + dpath*(dsintheta+dwidthhad*gRandom->Gaus());
0194       y = yhad + dpath*(dsinphi+dwidthhad*gRandom->Gaus());
0195       z = zhad+dpath * dcostheta*dcosphi;
0196       dtotehad += dhade;
0197       h2dhits->Fill(x,y,dhade);
0198       if(i%nlive == 0)
0199         {
0200           h2dlive->Fill(x,y,dhade);
0201           h3dlivehad->Fill(x,y,z);
0202         }
0203     }
0204
0205       dtotehad *= 1.0 + dresol * gRandom->Gaus();
0206       htotehad->Fill(dtotehad);
0207
0208       htote->Fill(dtotemip+dtoteem+dtotehad);
0209       h2dmipnonmip->Fill(dtoteem+dtotehad,dtotemip);
0210
0211       if(i%nlive==0)
0212     {
0213       c10->cd(1);
0214       h2dlive->Draw("colz");
0215       c10->cd(2);
0216       h3dlivemip->SetMarkerStyle(20);
0217       h3dlivemip->SetMarkerColor(1);
0218       h3dlivemip->SetMarkerSize(0.2);
0219       h3dlivemip->Draw();
0220       h3dliveem->SetMarkerStyle(20);
0221       h3dliveem->SetMarkerColor(2);
0222       h3dliveem->SetMarkerSize(0.4);
0223       h3dliveem->Draw("SAME");
0224       h3dlivehad->SetMarkerStyle(20);
0225       h3dlivehad->SetMarkerColor(3);
0226       h3dlivehad->SetMarkerSize(0.5);
0227       h3dlivehad->Draw("SAME");
0228       c10->Update();
0229       std::cout << i << std::endl;
0230       //      if(i==300) c10->Print("evt300.jpg","jpg");
0231
0232       h3dlivemip->Reset();
0233       h3dliveem->Reset();
0234       h3dlivehad->Reset();
0235       h2dlive->Reset();
0236
0237
0238
0239     }
0240
0241     }
0242
0243
0244 TFile  *fout = new TFile("fasthadv4_out.root","RECREATE");
0245
0246   h2dhits->Write();
0247   htote->Write();
0248   htotemip->Write();
0249   htoteem->Write();
0250   htotehad->Write();
0251   h2dmipnonmip->Write();
0252   h2dlive->Write();
0253   h3dlivemip->Write();
0254   h3dliveem->Write();
0255   h3dlivehad->Write();
0256   h3dliveall->Write();
0257
0258
0259
0260   fout->Close();
0261
0262 }
0263
0264
0265
0266
0267
0268
0269
0270