TPC/SCDCorrections/simRho.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include "sHelix.h"
0004 #include "sChargeMap.h"
0005
0006 #include "TStyle.h"
0007 #include "TCanvas.h"
0008 #include "TMath.h"
0009 #include "TRandom3.h"
0010 #include "TFile.h"
0011 #include "TH1D.h"
0012 #include "TH2D.h"
0013 #include "TH3D.h"
0014
0015 const int npid = 12+1;
0016 int apid[npid] =     {11, 13, // 2
0017                       211, // 1
0018                       321, // 1
0019                       411, 431, // 2
0020                       511, // 1
0021                       2212, // 1
0022                       3222, 3112, 3312, 3334, // 6
0023                       0};
0024 int acha[npid] = {-1,-1,+1,+1,+1,+1,+1,+1,+1,-1,-1,-1};
0025 TString spid[npid] = {"e^{-}",
0026                       "#mu^{-}",
0027                       "#pi^{+}",
0028                       "K^{+}",
0029                       "D^{+}",
0030                       "D^{+}_{s}",
0031                       "B^{+}",
0032                       "p",
0033                       "#Sigma^{+}",
0034                       "#Sigma^{-}",
0035                       "#Xi^{-}",
0036                       "#Omega^{-}",
0037
0038                       "Neutral"};
0039 TH1F *hpid      = new TH1F("hpid","hpid",npid,-0.5,npid-0.5);
0040 TH1F *hrpass    = new TH1F("hrpass","hrpass; r [cm]",100,-15,+15);
0041 TH1F *hzpass    = new TH1F("hzpass","hzpass; z [cm]",100,-15,+15);
0042 TH1F *hpidpass  = new TH1F("hpidpass","hpidpass",npid,-0.5,npid-0.5);
0043 TH1F *hptpass   = new TH1F("hptpass","hptpass; pt [GeV]",100,0,10);
0044 TH1F *hetapass  = new TH1F("hetapass","hetapass; eta",100,-3,+3);
0045 TH1F *hlengthpass=new TH1F("hlengthpass","hlengthpass; length [cm]",100,0,150);
0046 TH1F *hbxing    = new TH1F("hbxing","hbxing",150,-0.5,300-0.5);
0047 TH3F *tracking  = new TH3F("tracking","Generated tracks: tracking;xi;yi;zi",120,-100,+100,120,-100,+100,120,-100,+100);
0048
0049 int main() {
0050   for(int n=0; n!=npid; ++n) hpid->GetXaxis()->SetBinLabel(1+n,spid[n].Data());
0051   for(int n=0; n!=npid; ++n) hpidpass->GetXaxis()->SetBinLabel(1+n,spid[n].Data());
0052
0053   gStyle->SetOptStat(111111111);
0054   //================
0055   // Gas parameters
0056   //float gs_wf_me = 35; // [eV] mean energy for ionisation
0057   //float gs_density = 1e-3; // [gr/cm^3]
0058   //float gs_np = 14.35; // [cm^-1]
0059   float gs_nt = 47.80; // [cm^-1]
0060   float gs_ion_mobility = 4; // [cm^2/(V.s)]
0061
0062   //================
0063   // TPC running conditions
0064   float tpc_magnetic_field = 1.5; // [Tesla]
0065   float tpc_electric_field = 400; // [V/cm]
0066   float tpc_drift_velocity_e = 8; // [cm/us]
0067   float tpc_drift_velocity_i = tpc_electric_field*gs_ion_mobility*1e-3; // [cm/ms]
0068   int lumi = 189; //  1/(189*106e-9) => 49915.1 Hz
0069
0070   // input file
0071   std::ifstream inputf("/Users/cperez//ampt/ana/ampt.dat");
0072
0073   sChargeMap *map = new sChargeMap(100,30,80,
0074                    1,0,TMath::TwoPi(),
0075                    4*int(80/(tpc_drift_velocity_i*(lumi*106e-6)))+1,80,
0076                    tpc_drift_velocity_e,tpc_drift_velocity_i);
0077   // main routine
0078   int ncoll;
0079   for(ncoll=0; ncoll!=3; ++ncoll) {
0080     // =========
0081     // COLLISION
0082     if(ncoll%100==0) std::cout << ncoll << std::endl;
0083     // event wise
0084     int eventno, testno, nhadrons, npartA, npartB, npartElaA, npartIneA, npartElaB, npartIneB;
0085     float impactPar;
0086     // hadron wise
0087     int pid;
0088     float px, py, pz, mass, x, y, z, t;
0089     inputf >> eventno;
0090     if(!inputf.good()) {
0091       inputf.close();
0092       inputf.open("ampt.dat");
0093       inputf >> eventno;
0094     }
0095     inputf >> testno >> nhadrons >> impactPar >> npartA >> npartB >> npartElaA >> npartIneA >> npartElaB >> npartIneB;
0096     for(int np=0; np!=nhadrons; ++np) {
0097       inputf >> pid >> px >> py >> pz >> mass >> x >> y >> z >> t;
0098       if( TMath::AreEqualAbs(t,0,1e-3) ) continue;
0099       x *= 1e-13;
0100       y *= 1e-13;
0101       z *= 1e-13;
0102       int fill = npid-1;
0103       for(int n=0; n!=npid-1; ++n) if(apid[n]==TMath::Abs(pid)) fill = n;
0104       hpid->Fill( fill );
0105       if(fill==npid-1) continue;
0106       int q;
0107       for(int n=0; n!=npid-1; ++n) if(apid[n]==TMath::Abs(pid)) q = acha[n];
0108       sHelix a(x,y,z,px,py,pz,q,tpc_magnetic_field);
0109       float t1 = a.findFirstInterceptTo(30,80);
0110       if(t1>999) continue;
0111       float t2 = a.findFirstInterceptTo(80,80);
0112       if(t2>999) continue;
0113       float length = a.s(t1,t2);
0114       if(length<1) continue;
0115       hrpass->Fill( sqrt(x*x+y*y) );
0116       hzpass->Fill( z );
0117       hpidpass->Fill( fill );
0118       hptpass->Fill( sqrt(px*px+py*py) );
0119       float pmom = sqrt(px*px+py*py+pz*pz);
0120       float eta = 1.e30;
0121       if (pmom != TMath::Abs(pz)) eta = 0.5*TMath::Log((pmom+pz)/(pmom-pz));
0122       hetapass->Fill( eta );
0123       hlengthpass->Fill( length );
0124       float nprim = gRandom->Poisson(gs_nt*length);
0125       float track[100][3];
0126       a.breakIntoPieces(t1,t2,track);
0127       for(int i=0; i!=100; ++i) {
0128     float ri = TMath::Sqrt( track[i][0]*track[i][0] + track[i][1]*track[i][1] );
0129     if(ri<30-0.1) break;
0130     if(ri>80+0.1) break;
0131     float pi = 1;
0132         map->Fill(ri, pi, track[i][2], -nprim/100);
0133         map->Fill(ri, pi, track[i][2], nprim/100);
0134     //if(ncoll==1&&
0135     //   (
0136     //    np==622||
0137     //    np==935
0138     //    )) {
0139     tracking->Fill(track[i][0],track[i][1],track[i][2]);
0140     //  if(i==0) {
0141     //  std::cout << "pmom " << pmom << "| eta " << eta;
0142     //  std::cout << "| x y z " << x << " " << y << " " << z << " " << px << " " << py << " " << pz << std::endl;
0143     //  std::cout << "t " << t1 << "|" << t2 << std::endl;
0144     //}
0145     //std::cout << ri << " ";
0146     //if(i==99)
0147     //  std::cout << std::endl;
0148     //}
0149       }
0150     }
0151
0152     // ==========
0153     // DRIFT TIME
0154     int steps = gRandom->Poisson( lumi );
0155     hbxing->Fill(steps);
0156     float lapse = 106*steps*1e-6; // [ns] -> [ms]
0157     map->Propagate(lapse);
0158   }
0159   map->ScreenShot("map","root",ncoll);
0160   inputf.close();
0161   map->SaveIonMap("ionmap.root");
0162   //map->SaveRho("mark1_0.root",100,1,1600);
0163
0164   TFile *qa = new TFile("qa.root","RECREATE");
0165   qa->WriteTObject(hpid,hpid->GetName());
0166   qa->WriteTObject(hrpass,hrpass->GetName());
0167   qa->WriteTObject(hzpass,hzpass->GetName());
0168   qa->WriteTObject(hpidpass,hpidpass->GetName());
0169   qa->WriteTObject(hptpass,hptpass->GetName());
0170   qa->WriteTObject(hetapass,hetapass->GetName());
0171   qa->WriteTObject(hlengthpass,hlengthpass->GetName());
0172   qa->Close();
0173   tracking->SaveAs("tracking.root","root");
0174
0175   return 0;
0176 }