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 LogDebug("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 LogDebug("MaterialBudget") << name << " " << step << " " << matName
00101 << " " << stepLen << " " << step/radl << " "
00102 << radLen << " " << step/intl << " " << intLen;
00103 } else {
00104 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos: Step at " << name
00105 << " Length " << step << " in " << matName
00106 << " of density " << density << " g/cc;"
00107 << " 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 if (fillHistos) {
00156 if (id1 != 0) fillHisto(id1, 1);
00157 if (id2 != 0) fillHisto(id2, 1);
00158 }
00159 if (printSum) {
00160 for (unsigned int ii=0; ii<matList.size(); ii++) {
00161 edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
00162 << "\t" << radLength[ii] << "\t"
00163 << intLength[ii];
00164 }
00165 }
00166 }
00167
00168 void MaterialBudgetCastorHistos::book() {
00169
00170
00171 edm::Service<TFileService> tfile;
00172
00173 if ( !tfile.isAvailable() )
00174 throw cms::Exception("BadConfig") << "TFileService unavailable: "
00175 << "please add it to config file";
00176
00177 double maxPhi=pi;
00178 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
00179 << "histos === with " << binEta << " bins "
00180 << "in eta from " << etaLow << " to "
00181 << etaHigh << " and " << binPhi << " bins "
00182 << "in phi from " << -maxPhi << " to "
00183 << maxPhi;
00184
00185 char name[10], title[80], tag1[10], tag2[15];
00186
00187 for (int i=0; i<maxSet; i++) {
00188 double minEta=etaLow;
00189 double maxEta=etaHigh;
00190 int ireg = i;
00191 if (i > 9) {
00192 minEta = -etaHigh;
00193 maxEta = -etaLow;
00194 ireg -= 10;
00195 sprintf (tag2, " (-ve Eta Side)");
00196 } else {
00197 sprintf (tag2, " (+ve Eta Side)");
00198 }
00199 if ((i%2) == 0) {
00200 ireg /= 2;
00201 sprintf (tag1, " == Start");
00202 } else {
00203 ireg = (ireg-1)/2;
00204 sprintf (tag1, " == End");
00205 }
00206 sprintf(name, "%d", i+100);
00207 sprintf(title, "MB(X0) prof Eta in region %d%s%s", ireg, tag1, tag2);
00208 me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00209 sprintf(name, "%d", i+200);
00210 sprintf(title, "MB(L0) prof Eta in region %d%s%s", ireg, tag1, tag2);
00211 me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00212 sprintf(name, "%d", i+300);
00213 sprintf(title, "MB(Step) prof Eta in region %d%s%s", ireg, tag1, tag2);
00214 me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
00215 sprintf(name, "%d", i+400);
00216 sprintf(title, "Eta in region %d%s%s", ireg, tag1, tag2);
00217 me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
00218 sprintf(name, "%d", i+500);
00219 sprintf(title, "MB(X0) prof Ph in region %d%s%s", ireg, tag1, tag2);
00220 me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00221 sprintf(name, "%d", i+600);
00222 sprintf(title, "MB(L0) prof Ph in region %d%s%s", ireg, tag1, tag2);
00223 me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00224 sprintf(name, "%d", i+700);
00225 sprintf(title, "MB(Step) prof Ph in region %d%s%s", ireg, tag1, tag2);
00226 me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00227 sprintf(name, "%d", i+800);
00228 sprintf(title, "Phi in region %d%s%s", ireg, tag1, tag2);
00229 me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
00230 sprintf(name, "%d", i+900);
00231 sprintf(title, "MB(X0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00232 me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00233 binPhi/2, -maxPhi, maxPhi);
00234 sprintf(name, "%d", i+1000);
00235 sprintf(title, "MB(L0) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00236 me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00237 binPhi/2, -maxPhi, maxPhi);
00238 sprintf(name, "%d", i+1100);
00239 sprintf(title, "MB(Step) prof Eta Phi in region %d%s%s", ireg, tag1, tag2);
00240 me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, minEta, maxEta,
00241 binPhi/2, -maxPhi, maxPhi);
00242 sprintf(name, "%d", i+1200);
00243 sprintf(title, "Eta vs Phi in region %d%s%s", ireg, tag1, tag2);
00244 me1200[i]= tfile->make<TH2F>(name, title, binEta/2, minEta, maxEta,
00245 binPhi/2, -maxPhi, maxPhi);
00246 }
00247
00248 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
00249 << "histos done ===";
00250
00251 }
00252
00253 void MaterialBudgetCastorHistos::fillHisto(int id, int ix) {
00254
00255 int ii = 2*(id-1) + ix;
00256 double etaAbs = eta;
00257 if (eta < 0) {
00258 etaAbs = -eta;
00259 ii += 10;
00260 }
00261 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos:FillHisto called "
00262 << "with index " << id << ":" << ix << ":" << ii
00263 << " eta " << etaAbs << " (" << eta << ")"
00264 << " integrated step " << stepLen << " X0 "
00265 << radLen << " Lamda " << intLen;
00266
00267 me100[ii]->Fill(eta, radLen);
00268 me200[ii]->Fill(eta, intLen);
00269 me300[ii]->Fill(eta, stepLen);
00270 me400[ii]->Fill(eta);
00271
00272 if (etaAbs >= etaLow && etaAbs <= etaHigh) {
00273 me500[ii]->Fill(phi, radLen);
00274 me600[ii]->Fill(phi, intLen);
00275 me700[ii]->Fill(phi, stepLen);
00276 me800[ii]->Fill(phi);
00277 }
00278
00279 me900[ii]->Fill(eta, phi, radLen);
00280 me1000[ii]->Fill(eta, phi, intLen);
00281 me1100[ii]->Fill(eta, phi, stepLen);
00282 me1200[ii]->Fill(eta, phi);
00283
00284 }
1.7.3