00001 #include "Validation/Geometry/interface/MaterialBudgetCastorHistos.h"
00002
00003 #include "DetectorDescription/Core/interface/DDFilter.h"
00004 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
00005 #include "DetectorDescription/Core/interface/DDSplit.h"
00006 #include "DetectorDescription/Core/interface/DDValue.h"
00007 #include "FWCore/Utilities/interface/Exception.h"
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00011
00012 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00013 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00014
00015 MaterialBudgetCastorHistos::MaterialBudgetCastorHistos(const edm::ParameterSet &p){
00016
00017 binEta = p.getUntrackedParameter<int>("NBinEta", 100);
00018 binPhi = p.getUntrackedParameter<int>("NBinPhi", 180);
00019 etaLow = p.getUntrackedParameter<double>("EtaLow", 5.0);
00020 etaHigh = p.getUntrackedParameter<double>("EtaHigh", 7.0);
00021 fillHistos = p.getUntrackedParameter<bool>("FillHisto", true);
00022 printSum = p.getUntrackedParameter<bool>("PrintSummary", false);
00023 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: FillHisto : "
00024 << fillHistos << " PrintSummary " << printSum
00025 << " == Eta plot: NX " << binEta << " Range "
00026 << etaLow << " (" << -etaHigh << ") : "
00027 << etaHigh << " (" << -etaLow <<") Phi plot: "
00028 << "NX " << binPhi << " Range " << -pi << ":"
00029 << pi << " (Eta limit " << etaLow << ":"
00030 << etaHigh <<")";
00031 if (fillHistos) book();
00032
00033 }
00034
00035 MaterialBudgetCastorHistos::~MaterialBudgetCastorHistos() {
00036 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Save user "
00037 << "histos ===";
00038 }
00039
00040 void MaterialBudgetCastorHistos::fillStartTrack(const G4Track* aTrack) {
00041
00042 id1 = id2 = steps = 0;
00043 radLen = intLen = stepLen = 0;
00044
00045 const G4ThreeVector& dir = aTrack->GetMomentum() ;
00046 if (dir.theta() != 0 ) {
00047 eta = dir.eta();
00048 } else {
00049 eta = -99;
00050 }
00051 phi = dir.phi();
00052 double theEnergy = aTrack->GetTotalEnergy();
00053 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
00054
00055 if (printSum) {
00056 matList.clear();
00057 stepLength.clear();
00058 radLength.clear();
00059 intLength.clear();
00060 }
00061
00062 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Track "
00063 << aTrack->GetTrackID() << " Code " << theID
00064 << " Energy " << theEnergy/GeV << " GeV; Eta "
00065 << eta << " Phi " << phi/deg << " PT "
00066 << dir.perp()/GeV << " GeV *****";
00067 }
00068
00069
00070 void MaterialBudgetCastorHistos::fillPerStep(const G4Step* aStep) {
00071
00072 G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
00073 double step = aStep->GetStepLength();
00074 double radl = material->GetRadlen();
00075 double intl = material->GetNuclearInterLength();
00076 double density = material->GetDensity() / (g/cm3);
00077
00078 int id1Old = id1;
00079 int id2Old = id2;
00080 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00081 std::string name = touch->GetVolume(0)->GetName();
00082 std::string matName = material->GetName();
00083 if (printSum) {
00084 bool found = false;
00085 for (unsigned int ii=0; ii<matList.size(); ii++) {
00086 if (matList[ii] == matName) {
00087 stepLength[ii] += step;
00088 radLength[ii] += (step/radl);
00089 intLength[ii] += (step/intl);
00090 found = true;
00091 break;
00092 }
00093 }
00094 if (!found) {
00095 matList.push_back(matName);
00096 stepLength.push_back(step);
00097 radLength.push_back(step/radl);
00098 intLength.push_back(step/intl);
00099 }
00100 edm::LogInfo("MaterialBudget") << name << " " << step << " " << matName
00101 << " " << stepLen << " " << step/radl << " "
00102 << radLen << " " <<step/intl << " " <<intLen;
00103 } else {
00104 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Step at "
00105 << name << " Length " << step << " in "
00106 << matName << " of density " << density
00107 << " g/cc; Radiation Length " <<radl <<" mm;"
00108 << " Interaction Length " << intl << " mm\n"
00109 << " Position "
00110 << aStep->GetPreStepPoint()->GetPosition()
00111 << " Cylindrical R "
00112 <<aStep->GetPreStepPoint()->GetPosition().perp()
00113 << " Length (so far) " << stepLen << " L/X0 "
00114 << step/radl << "/" << radLen << " L/Lambda "
00115 << step/intl << "/" << intLen;
00116 }
00117
00118 int level = ((touch->GetHistoryDepth())+1);
00119 std::string name1="XXXX", name2="XXXX";
00120 if (level>3) name1 = touch->GetVolume(level-4)->GetName();
00121 if (level>4) name2 = touch->GetVolume(level-5)->GetName();
00122 if (name1 == "CAST") {
00123 id1 = 1;
00124 if (name2 == "CAEC") id2 = 2;
00125 else if (name2 == "CAHC") id2 = 3;
00126 else if (name2 == "CEDC") id2 = 4;
00127 else if (name2 == "CHDC") id2 = 5;
00128 else id2 = 0;
00129 } else {
00130 id1 = id2 = 0;
00131 }
00132 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos: Level " << level
00133 << " Volume " << name1 << " and " << name2
00134 << " ID1 " << id1 << " (" << id1Old << ") ID2 "
00135 << id2 << " (" << id2Old << ")";
00136
00137 if (fillHistos) {
00138 if (id1 != id1Old) {
00139 if (id1 == 0) fillHisto(id1Old, 1);
00140 else fillHisto(id1, 0);
00141 }
00142 if (id2 != id2Old) {
00143 if (id2 == 0) fillHisto(id2Old, 1);
00144 else fillHisto(id2, 0);
00145 }
00146 }
00147
00148 stepLen += step;
00149 radLen += step/radl;
00150 intLen += step/intl;
00151 }
00152
00153
00154 void MaterialBudgetCastorHistos::fillEndTrack() {
00155
00156 if (fillHistos) {
00157 if (id1 != 0) fillHisto(id1, 1);
00158 if (id2 != 0) fillHisto(id2, 1);
00159 }
00160 if (printSum) {
00161 for (unsigned int ii=0; ii<matList.size(); ii++) {
00162 edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
00163 << "\t" << radLength[ii] << "\t"
00164 << intLength[ii];
00165 }
00166 }
00167 }
00168
00169 void MaterialBudgetCastorHistos::book() {
00170
00171
00172 edm::Service<TFileService> tfile;
00173
00174 if ( !tfile.isAvailable() )
00175 throw cms::Exception("BadConfig") << "TFileService unavailable: "
00176 << "please add it to config file";
00177
00178 double maxPhi=pi;
00179 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
00180 << "histos === with " << binEta << " bins "
00181 << "in eta from " << etaLow << " to "
00182 << etaHigh << " and " << binPhi << " bins "
00183 << "in phi from " << -maxPhi << " to "
00184 << maxPhi;
00185
00186 char name[10], title[80], tag1[10], tag2[15];
00187
00188 for (int i=0; i<maxSet; i++) {
00189 double minEta=etaLow;
00190 double maxEta=etaHigh;
00191 int ireg = i;
00192 if (i > 9) {
00193 minEta = -etaHigh;
00194 maxEta = -etaLow;
00195 ireg -= 10;
00196 sprintf (tag2, " (-ve Eta Side)");
00197 } else {
00198 sprintf (tag2, " (+ve Eta Side)");
00199 }
00200 if ((i%2) == 0) {
00201 ireg /= 2;
00202 sprintf (tag1, " == Start");
00203 } else {
00204 ireg = (ireg-1)/2;
00205 sprintf (tag1, " == End");
00206 }
00207 sprintf(name, "%d", i+100);
00208 sprintf(title, "MB(X0) prof Eta in region %d%s%s", ireg, tag1, tag2);
00209 me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00210 sprintf(name, "%d", i+200);
00211 sprintf(title, "MB(L0) prof Eta in region %d%s%s", ireg, tag1, tag2);
00212 me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00213 sprintf(name, "%d", i+300);
00214 sprintf(title, "MB(Step) prof Eta in region %d%s%s", ireg, tag1, tag2);
00215 me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00216 sprintf(name, "%d", i+400);
00217 sprintf(title, "Eta in region %d%s%s", ireg, tag1, tag2);
00218 me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
00219 sprintf(name, "%d", i+500);
00220 sprintf(title, "MB(X0) prof Ph in region %d%s%s", ireg, tag1, tag2);
00221 me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00222 sprintf(name, "%d", i+600);
00223 sprintf(title, "MB(L0) prof Ph in region %d%s%s", ireg, tag1, tag2);
00224 me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00225 sprintf(name, "%d", i+700);
00226 sprintf(title, "MB(Step) prof Ph in region %d%s%s", ireg, tag1, tag2);
00227 me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00228 sprintf(name, "%d", i+800);
00229 sprintf(title, "Phi in region %d%s%s", ireg, tag1, tag2);
00230 me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
00231 sprintf(name, "%d", i+900);
00232 sprintf(title, "MB(X0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00233 me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00234 binPhi/2, -maxPhi, maxPhi);
00235 sprintf(name, "%d", i+1000);
00236 sprintf(title, "MB(L0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00237 me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00238 binPhi/2, -maxPhi, maxPhi);
00239 sprintf(name, "%d", i+1100);
00240 sprintf(title, "MB(Step) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00241 me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00242 binPhi/2, -maxPhi, maxPhi);
00243 sprintf(name, "%d", i+1200);
00244 sprintf(title, "Eta vs Phi in region %d%s%s", ireg, tag1, tag2);
00245 me1200[i]= tfile->make<TH2F>(name, title, binEta/2, minEta, maxEta,
00246 binPhi/2, -maxPhi, maxPhi);
00247 }
00248
00249 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
00250 << "histos done ===";
00251
00252 }
00253
00254 void MaterialBudgetCastorHistos::fillHisto(int id, int ix) {
00255
00256 int ii = 2*(id-1) + ix;
00257 double etaAbs = eta;
00258 if (eta < 0) {
00259 etaAbs = -eta;
00260 ii += 10;
00261 }
00262 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos:FillHisto "
00263 << "called with index " << id << ":" << ix
00264 << ":" << ii << " eta " << etaAbs << " ("
00265 << eta << ") integrated step " << stepLen
00266 << " X0 " << radLen << " Lamda " << intLen;
00267
00268 me100[ii]->Fill(eta, radLen);
00269 me200[ii]->Fill(eta, intLen);
00270 me300[ii]->Fill(eta, stepLen);
00271 me400[ii]->Fill(eta);
00272
00273 if (etaAbs >= etaLow && etaAbs <= etaHigh) {
00274 me500[ii]->Fill(phi, radLen);
00275 me600[ii]->Fill(phi, intLen);
00276 me700[ii]->Fill(phi, stepLen);
00277 me800[ii]->Fill(phi);
00278 }
00279
00280 me900[ii]->Fill(eta, phi, radLen);
00281 me1000[ii]->Fill(eta, phi, intLen);
00282 me1100[ii]->Fill(eta, phi, stepLen);
00283 me1200[ii]->Fill(eta, phi);
00284
00285 }
1.7.3