/scratch0/jsnuveri/BDSIM/BDSIMgit/bdsim/src/BDSOutputROOT.cc

00001 #ifdef USE_ROOT
00002 #include "BDSOutputROOT.hh"
00003 
00004 #include "BDSDebug.hh"
00005 #include "BDSExecOptions.hh"
00006 #include "BDSSampler.hh"
00007 #include "BDSSamplerCylinder.hh"
00008 #include "BDSTrajectory.hh"
00009 #include "BDSUtilities.hh"
00010 #include "BDSHistogram.hh"
00011 
00012 
00013 BDSOutputROOT::BDSOutputROOT():BDSOutputBase()
00014 {
00015 #ifdef BDSDEBUG
00016   G4cout<<"output format ROOT"<<G4endl;
00017 #endif
00018   theRootOutputFile = NULL;
00019   PrecisionRegionEnergyLossTree = NULL;
00020   EnergyLossTree    = NULL;
00021   PrimaryHitsTree   = NULL;
00022   PrimaryLossTree   = NULL;
00023   Init(); // activate the output
00024 }
00025 
00026 BDSOutputROOT::~BDSOutputROOT()
00027 {
00028   if (theRootOutputFile && theRootOutputFile->IsOpen()) {
00029     theRootOutputFile->Write();
00030   }
00031 }
00032 
00033 void BDSOutputROOT::BuildSamplerTree(G4String name){
00034   TTree* SamplerTree = new TTree(name, "Sampler output");
00035   
00036   SamplerTree->Branch("E0",&E0,"E0/F"); // (GeV)
00037   SamplerTree->Branch("x0",&x0,"x0/F"); // (mum)
00038   SamplerTree->Branch("y0",&y0,"y0/F"); // (mum)
00039   SamplerTree->Branch("z0",&z0,"z0/F"); // (m)
00040   SamplerTree->Branch("xp0",&xp0,"xp0/F"); // (rad)
00041   SamplerTree->Branch("yp0",&yp0,"yp0/F"); // (rad)
00042   SamplerTree->Branch("zp0",&zp0,"zp0/F"); // (rad)
00043   SamplerTree->Branch("t0",&t0,"t0/F"); // (ns)
00044 
00045   SamplerTree->Branch("E_prod",&E_prod,"E_prod/F"); // (GeV)
00046   SamplerTree->Branch("x_prod",&x_prod,"x_prod/F"); // (mum)
00047   SamplerTree->Branch("y_prod",&y_prod,"y_prod/F"); // (mum)
00048   SamplerTree->Branch("z_prod",&z_prod,"z_prod/F"); // (m)
00049   SamplerTree->Branch("xp_prod",&xp_prod,"xp_prod/F"); // (rad)
00050   SamplerTree->Branch("yp_prod",&yp_prod,"yp_prod/F"); // (rad)
00051   SamplerTree->Branch("zp_prod",&zp_prod,"zp_prod/F"); // (rad)
00052   SamplerTree->Branch("t_prod",&t_prod,"t_prod/F"); // (ns)
00053 
00054   SamplerTree->Branch("E_lastScat",&E_lastScat,"E_lastScat/F"); // (GeV)
00055   SamplerTree->Branch("x_lastScat",&x_lastScat,"x_lastScat/F"); // (mum)
00056   SamplerTree->Branch("y_lastScat",&y_lastScat,"y_lastScat/F"); // (mum)
00057   SamplerTree->Branch("z_lastScat",&z_lastScat,"z_lastScat/F"); // (m)
00058   SamplerTree->Branch("xp_lastScat",&xp_lastScat,"xp_lastScat/F"); // (rad)
00059   SamplerTree->Branch("yp_lastScat",&yp_lastScat,"yp_lastScat/F"); // (rad)
00060   SamplerTree->Branch("zp_lastScat",&zp_lastScat,"zp_lastScat/F"); // (rad)
00061   SamplerTree->Branch("t_lastScat",&t_lastScat,"t_lastScat/F"); // (ns)
00062   
00063   SamplerTree->Branch("E",&E,"E/F"); // (GeV)
00064   SamplerTree->Branch("x",&x,"x/F"); // (mum)
00065   SamplerTree->Branch("y",&y,"y/F"); // (mum)
00066   SamplerTree->Branch("z",&z,"z/F"); // (m)
00067   SamplerTree->Branch("xp",&xp,"xp/F"); // (rad)
00068   SamplerTree->Branch("yp",&yp,"yp/F"); // (rad)
00069   SamplerTree->Branch("zp",&zp,"zp/F"); // (rad)
00070   SamplerTree->Branch("t",&t,"t/F"); // (ns)
00071   
00072   SamplerTree->Branch("X",&X,"X/F"); // (mum)
00073   SamplerTree->Branch("Y",&Y,"Y/F"); // (mum)
00074   SamplerTree->Branch("Z",&Z,"Z/F"); // (m)
00075   SamplerTree->Branch("Xp",&Xp,"Xp/F"); // (rad)
00076   SamplerTree->Branch("Yp",&Yp,"Yp/F"); // (rad)
00077   SamplerTree->Branch("Zp",&Zp,"Zp/F"); // (rad)
00078   
00079   SamplerTree->Branch("s",&s,"s/F"); // (m)
00080   
00081   SamplerTree->Branch("weight",&weight,"weight/F");
00082   SamplerTree->Branch("partID",&part,"partID/I");
00083   SamplerTree->Branch("nEvent",&nev,"nEvent/I");
00084   SamplerTree->Branch("parentID",&pID,"parentID/I");
00085   SamplerTree->Branch("trackID",&track_id,"trackID/I");
00086   SamplerTree->Branch("turnnumber",&turnnumber,"turnnumber/I");
00087 }
00088 
00089 void BDSOutputROOT::Init()
00090 {
00091   // set up the root file
00092   filename = BDSExecOptions::Instance()->GetOutputFilename() + "_" 
00093     + BDSGlobalConstants::Instance()->StringFromInt(outputFileNumber) + ".root";
00094   
00095   G4cout<<"Setting up new file: "<<filename<<G4endl;
00096   theRootOutputFile=new TFile(filename,"RECREATE", "BDS output file");
00097 
00098   //build sampler tree
00099   G4String primariesSamplerName="primaries";
00100 #ifdef BDSDEBUG
00101   G4cout << __METHOD_NAME__ << " building sampler tree named: " << primariesSamplerName << G4endl;
00102 #endif
00103   BuildSamplerTree(primariesSamplerName);
00104   for(G4int i=0;i<BDSSampler::GetNSamplers();i++)
00105     {
00106 #ifdef BDSDEBUG
00107       G4cout << __METHOD_NAME__ << " building sampler tree number: " << i << G4endl;
00108 #endif
00109       G4String name=BDSSampler::outputNames[i];
00110 #ifdef BDSDEBUG
00111       G4cout << __METHOD_NAME__ << " named: " << name << G4endl;
00112 #endif
00113       BuildSamplerTree(name);
00114     }
00115   for(G4int i=0;i<BDSSamplerCylinder::GetNSamplers();i++)
00116     {
00117       //G4String name="samp"+BDSGlobalConstants::Instance()->StringFromInt(i+1);
00118       G4String name=BDSSamplerCylinder::outputNames[i];
00119       BuildSamplerTree(name);
00120     }
00121 
00122   if(BDSGlobalConstants::Instance()->GetStoreTrajectory() || BDSGlobalConstants::Instance()->GetStoreMuonTrajectories() || BDSGlobalConstants::Instance()->GetStoreNeutronTrajectories()) 
00123     // create a tree with trajectories
00124     {
00125       TTree* TrajTree = new TTree("Trajectories", "Trajectories");
00126       TrajTree->Branch("x",&x,"x/F"); // (mum)
00127       TrajTree->Branch("y",&y,"y/F"); // (mum)
00128       TrajTree->Branch("z",&z,"z/F"); // (m)
00129       TrajTree->Branch("part",&part,"part/I");
00130     }
00131 
00132   // build energy loss histogram  
00133   EnergyLossTree= new TTree("ElossTree", "Energy Loss");
00134   EnergyLossTree->Branch("s",&S_el,"s/F"); // (m)
00135   EnergyLossTree->Branch("E",&E_el,"E/F"); // (GeV)
00136 
00137   //Primary loss tree setup
00138   PrimaryLossTree  = new TTree("PlossTree", "Primary Losses");
00139   PrimaryLossTree->Branch("X",          &X_pl,          "X/F"); // (m)
00140   PrimaryLossTree->Branch("Y",          &Y_pl,          "Y/F"); // (m)
00141   PrimaryLossTree->Branch("Z",          &Z_pl,          "Z/F"); // (m)
00142   PrimaryLossTree->Branch("S",          &S_pl,          "S/F"); // (m)
00143   PrimaryLossTree->Branch("x",          &x_pl,          "x/F"); // (m)
00144   PrimaryLossTree->Branch("y",          &y_pl,          "y/F"); // (m)
00145   PrimaryLossTree->Branch("z",          &z_pl,          "z/F"); // (m)
00146   PrimaryLossTree->Branch("E",          &E_pl,          "s/F"); // (GeV)
00147   PrimaryLossTree->Branch("weight",     &weight_pl,     "weight/F");
00148   PrimaryLossTree->Branch("partID",     &part_pl,       "partID/I");
00149   PrimaryLossTree->Branch("turnnumber", &turnnumber_pl, "turnnumber/I");
00150   PrimaryLossTree->Branch("eventNo",    &eventno_pl,    "eventNo/I");
00151 
00152   //Primary hits tree setup
00153   PrimaryHitsTree  = new TTree("PhitsTree", "Primary Hits");
00154   PrimaryHitsTree->Branch("X",          &X_ph,          "X/F"); // (m)
00155   PrimaryHitsTree->Branch("Y",          &Y_ph,          "Y/F"); // (m)
00156   PrimaryHitsTree->Branch("Z",          &Z_ph,          "Z/F"); // (m)
00157   PrimaryHitsTree->Branch("S",          &S_ph,          "S/F"); // (m)
00158   PrimaryHitsTree->Branch("x",          &x_ph,          "x/F"); // (m)
00159   PrimaryHitsTree->Branch("y",          &y_ph,          "y/F"); // (m)
00160   PrimaryHitsTree->Branch("z",          &z_ph,          "z/F"); // (m)
00161   PrimaryHitsTree->Branch("E",          &E_ph,          "E/F"); // (GeV)
00162   PrimaryHitsTree->Branch("weight",     &weight_ph,     "weight/F");
00163   PrimaryHitsTree->Branch("partID",     &part_ph,       "partID/I");
00164   PrimaryHitsTree->Branch("turnnumber", &turnnumber_ph, "turnnumber/I");
00165   PrimaryHitsTree->Branch("eventNo",    &eventno_ph,    "eventNo/I");
00166   
00167   //Precision region energy loss tree setup
00168   PrecisionRegionEnergyLossTree= new TTree("PrecisionRegionElossTree", "Energy Loss");
00169   PrecisionRegionEnergyLossTree->Branch("X",          &X_el_p,          "X/F"); // (m)
00170   PrecisionRegionEnergyLossTree->Branch("Y",          &Y_el_p,          "Y/F"); // (m)
00171   PrecisionRegionEnergyLossTree->Branch("Z",          &Z_el_p,          "Z/F"); // (m)
00172   PrecisionRegionEnergyLossTree->Branch("S",          &S_el_p,          "S/F"); // (m)
00173   PrecisionRegionEnergyLossTree->Branch("x",          &x_el_p,          "x/F"); // (m)
00174   PrecisionRegionEnergyLossTree->Branch("y",          &y_el_p,          "y/F"); // (m)
00175   PrecisionRegionEnergyLossTree->Branch("z",          &z_el_p,          "z/F"); // (m)
00176   PrecisionRegionEnergyLossTree->Branch("E",          &E_el_p,          "E/F"); // (GeV)
00177   PrecisionRegionEnergyLossTree->Branch("weight",     &weight_el_p,     "weight/F");
00178   PrecisionRegionEnergyLossTree->Branch("partID",     &part_el_p,       "partID/I");
00179   PrecisionRegionEnergyLossTree->Branch("volumeName", &volumeName_el_p, "volumeName/C");
00180   PrecisionRegionEnergyLossTree->Branch("turnnumber", &turnnumber_el_p, "turnnumber/I");
00181   PrecisionRegionEnergyLossTree->Branch("eventNo",    &eventno_el_p,    "eventNo/I");
00182 }
00183 
00184 void BDSOutputROOT::WriteRootHit(G4String Name, 
00185                                  G4double InitMom, 
00186                                  G4double InitX, 
00187                                  G4double InitY, 
00188                                  G4double InitZ, 
00189                                  G4double InitXPrime, 
00190                                  G4double InitYPrime, 
00191                                  G4double InitZPrime, 
00192                                  G4double InitT, 
00193                                  G4double ProdMom, 
00194                                  G4double ProdX, 
00195                                  G4double ProdY, 
00196                                  G4double ProdZ, 
00197                                  G4double ProdXPrime, 
00198                                  G4double ProdYPrime, 
00199                                  G4double ProdZPrime, 
00200                                  G4double ProdT, 
00201                                  G4double LastScatMom, 
00202                                  G4double LastScatX, 
00203                                  G4double LastScatY, 
00204                                  G4double LastScatZ, 
00205                                  G4double LastScatXPrime, 
00206                                  G4double LastScatYPrime, 
00207                                  G4double LastScatZPrime, 
00208                                  G4double LastScatT, 
00209                                  G4double Mom, 
00210                                  G4double X, 
00211                                  G4double Y, 
00212                                  G4double Z, 
00213                                  G4double XPrime, 
00214                                  G4double YPrime, 
00215                                  G4double ZPrime, 
00216                                  G4double T, 
00217                                  G4double GlobalX, 
00218                                  G4double GlobalY, 
00219                                  G4double GlobalZ, 
00220                                  G4double GlobalXPrime, 
00221                                  G4double GlobalYPrime, 
00222                                  G4double GlobalZPrime, 
00223                                  G4double S, 
00224                                  G4double Weight, 
00225                                  G4int    PDGtype, 
00226                                  G4int    EventNo, 
00227                                  G4int    ParentID,
00228                                  G4int    TrackID, 
00229                                  G4int    TurnsTaken){
00230 
00231   TTree* sTree=(TTree*)gDirectory->Get(Name);
00232   if(!sTree) G4Exception("BDSOutputROOT: ROOT Sampler not found!", "-1", FatalException, "");
00233   E0=InitMom / CLHEP::GeV;
00234   x0=InitX / CLHEP::micrometer;
00235   y0=InitY / CLHEP::micrometer;
00236   z0=InitZ / CLHEP::m;
00237   xp0=InitXPrime / CLHEP::radian;
00238   yp0=InitYPrime / CLHEP::radian;
00239   zp0=InitZPrime / CLHEP::radian;
00240   t0=InitT / CLHEP::ns;
00241   E_prod=ProdMom / CLHEP::GeV;
00242   x_prod=ProdX / CLHEP::micrometer;
00243   y_prod=ProdY / CLHEP::micrometer;
00244   z_prod=ProdZ / CLHEP::m;
00245   xp_prod=ProdXPrime / CLHEP::radian;
00246   yp_prod=ProdYPrime / CLHEP::radian;
00247   zp_prod=ProdZPrime / CLHEP::radian;
00248   t_prod=ProdT / CLHEP::ns;
00249   E_lastScat=LastScatMom / CLHEP::GeV;
00250   x_lastScat=LastScatX / CLHEP::micrometer;
00251   y_lastScat=LastScatY / CLHEP::micrometer;
00252   z_lastScat=LastScatZ / CLHEP::m;
00253   xp_lastScat=LastScatXPrime / CLHEP::radian;
00254   yp_lastScat=LastScatYPrime / CLHEP::radian;
00255   zp_lastScat=LastScatZPrime / CLHEP::radian;
00256   t_lastScat=LastScatT / CLHEP::ns;
00257 
00258   E=Mom / CLHEP::GeV;
00259   //Edep=Edep / CLHEP::GeV;
00260   x=X / CLHEP::micrometer;
00261   y=Y / CLHEP::micrometer;
00262   z=Z / CLHEP::m;
00263   xp=XPrime / CLHEP::radian;
00264   yp=YPrime / CLHEP::radian;
00265   zp=ZPrime / CLHEP::radian;
00266   t=T / CLHEP::ns;
00267   X=GlobalX / CLHEP::m;
00268   Y=GlobalY / CLHEP::m;
00269   Z=GlobalZ / CLHEP::m;
00270   Xp=GlobalXPrime / CLHEP::radian;
00271   Yp=GlobalYPrime / CLHEP::radian;
00272   Zp=GlobalZPrime / CLHEP::radian;
00273   s=S / CLHEP::m;
00274   weight=Weight;
00275   part=PDGtype; 
00276   nev=EventNo; 
00277   pID=ParentID; 
00278   track_id=TrackID;
00279   turnnumber=TurnsTaken;
00280   sTree->Fill();
00281 }
00282 
00283 void BDSOutputROOT::WritePrimary(G4String samplerName, 
00284                                  G4double E,
00285                                  G4double x0,
00286                                  G4double y0,
00287                                  G4double z0,
00288                                  G4double xp,
00289                                  G4double yp,
00290                                  G4double zp,
00291                                  G4double t,
00292                                  G4double weight,
00293                                  G4int    PDGType, 
00294                                  G4int    nEvent, 
00295                                  G4int    TurnsTaken){
00296   WriteRootHit(samplerName, 
00297                E, 
00298                x0, y0, z0, 
00299                xp, yp, zp, 
00300                t, 
00301                E, 
00302                x0, y0, z0, 
00303                xp, yp, zp, 
00304                t, 
00305                E, 
00306                x0, y0, z0, 
00307                xp, yp, zp, 
00308                t, 
00309                E, 
00310                x0, y0, z0, 
00311                xp, yp, zp, 
00312                t, 
00313                0,0,0,0,0,0,0,
00314                weight, 
00315                PDGType, 
00316                nEvent, 
00317                0, 
00318                1, 
00319                TurnsTaken);
00320 }
00321 
00322 void BDSOutputROOT::WriteHits(BDSSamplerHitsCollection *hc)
00323 {
00324   G4String name;
00325 #ifdef BDSDEBUG
00326   G4cout << __METHOD_NAME__ << " hc->entries() = " << hc->entries() << G4endl;
00327 #endif
00328   for (G4int i=0; i<hc->entries(); i++)
00329     {
00330       G4String name = (*hc)[i]->GetName();
00331 #ifdef BDSDEBUG
00332       G4cout << "Writing hit to sampler " << name << G4endl;
00333 #endif
00334       WriteRootHit(name,
00335                    (*hc)[i]->GetInitMom(),
00336                    (*hc)[i]->GetInitX(),
00337                    (*hc)[i]->GetInitY(),
00338                    (*hc)[i]->GetInitZ(),
00339                    (*hc)[i]->GetInitXPrime(),
00340                    (*hc)[i]->GetInitYPrime(),
00341                    (*hc)[i]->GetInitZPrime(),
00342                    (*hc)[i]->GetInitT(),
00343                    (*hc)[i]->GetProdMom(),
00344                    (*hc)[i]->GetProdX(),
00345                    (*hc)[i]->GetProdY(),
00346                    (*hc)[i]->GetProdZ(),
00347                    (*hc)[i]->GetProdXPrime(),
00348                    (*hc)[i]->GetProdYPrime(),
00349                    (*hc)[i]->GetProdZPrime(),
00350                    (*hc)[i]->GetProdT(),
00351                    (*hc)[i]->GetLastScatMom(),
00352                    (*hc)[i]->GetLastScatX(),
00353                    (*hc)[i]->GetLastScatY(),
00354                    (*hc)[i]->GetLastScatZ(),
00355                    (*hc)[i]->GetLastScatXPrime(),
00356                    (*hc)[i]->GetLastScatYPrime(),
00357                    (*hc)[i]->GetLastScatZPrime(),
00358                    (*hc)[i]->GetLastScatT(),
00359                    (*hc)[i]->GetMom(),
00360                    (*hc)[i]->GetX(),
00361                    (*hc)[i]->GetY(),
00362                    (*hc)[i]->GetZ(),
00363                    (*hc)[i]->GetXPrime(),
00364                    (*hc)[i]->GetYPrime(),
00365                    (*hc)[i]->GetZPrime(),
00366                    (*hc)[i]->GetT(),
00367                    (*hc)[i]->GetGlobalX(),
00368                    (*hc)[i]->GetGlobalY(),
00369                    (*hc)[i]->GetGlobalZ(),
00370                    (*hc)[i]->GetGlobalXPrime(),
00371                    (*hc)[i]->GetGlobalYPrime(),
00372                    (*hc)[i]->GetGlobalZPrime(),
00373                    (*hc)[i]->GetS(),
00374                    (*hc)[i]->GetWeight(),
00375                    (*hc)[i]->GetPDGtype(), 
00376                    (*hc)[i]->GetEventNo(), 
00377                    (*hc)[i]->GetParentID(), 
00378                    (*hc)[i]->GetTrackID(),
00379                    (*hc)[i]->GetTurnsTaken()
00380                    );
00381     }
00382 }
00383 
00385 void BDSOutputROOT::WriteTrajectory(std::vector<BDSTrajectory*> &TrajVec){
00386   //  G4int tID;
00387   G4TrajectoryPoint* TrajPoint;
00388   G4ThreeVector TrajPos;
00389   
00390   TTree* TrajTree;
00391   
00392   G4String name = "Trajectories";
00393   
00394   TrajTree=(TTree*)gDirectory->Get(name);
00395   
00396   if(TrajTree == NULL) { G4cerr<<"TrajTree=NULL"<<G4endl; return;}
00397   
00398   if(TrajVec.size())
00399     {
00400       std::vector<BDSTrajectory*>::iterator iT;
00401       
00402       for(iT=TrajVec.begin();iT<TrajVec.end();iT++)
00403         {
00404           G4Trajectory* Traj=(G4Trajectory*)(*iT);
00405           
00406           //      tID=Traj->GetTrackID();             
00407           
00408           G4int parentID=Traj->GetParentID();
00409           part = Traj->GetPDGEncoding();
00410           
00411           G4bool saveTrajectory=true;
00412           
00413           if(!((parentID==0)&&(BDSGlobalConstants::Instance()->GetStoreTrajectory()))){ 
00414             if(!((std::abs(part)==13)&&(BDSGlobalConstants::Instance()->GetStoreMuonTrajectories()))){ 
00415               if(!((part==2112)&&(BDSGlobalConstants::Instance()->GetStoreNeutronTrajectories()))){ 
00416                 saveTrajectory=false;
00417               }
00418             }
00419           }
00420           
00421           if(saveTrajectory){
00422             for(G4int j=0; j<Traj->GetPointEntries(); j++)
00423               {
00424                 TrajPoint=(G4TrajectoryPoint*)Traj->GetPoint(j);
00425                 TrajPos=TrajPoint->GetPosition();
00426                 
00427                 x = TrajPos.x() / CLHEP::m;
00428                 y = TrajPos.y() / CLHEP::m;
00429                 z = TrajPos.z() / CLHEP::m;
00430                 
00431                 TrajTree->Fill();
00432               }
00433           }
00434         }
00435     }
00436 }
00437 
00438 void BDSOutputROOT::WriteEnergyLoss(BDSEnergyCounterHitsCollection* hc)
00439 {
00440   G4int n_hit = hc->entries();
00441   for (G4int i=0;i<n_hit;i++)
00442     {
00443       //all regions fill the energy loss tree....
00444       E_el = (*hc)[i]->GetEnergy()/CLHEP::GeV;
00445       S_el = (*hc)[i]->GetS()/CLHEP::m;
00446       EnergyLossTree->Fill();
00447       
00448       if((*hc)[i]->GetPrecisionRegion()){ //Only the precision region fills this tree, preserving every hit, its position and weight, instead of summing weighted energy in each beam line component.
00449         weight_el_p  = (*hc)[i]->GetWeight();
00450         E_el_p       = (*hc)[i]->GetEnergy()/CLHEP::GeV;
00451         X_el_p       = (*hc)[i]->GetX()/CLHEP::m;
00452         Y_el_p       = (*hc)[i]->GetY()/CLHEP::m;
00453         Z_el_p       = (*hc)[i]->GetZ()/CLHEP::m;
00454         S_el_p       = (*hc)[i]->GetS()/CLHEP::m;
00455         x_el_p       = (*hc)[i]->Getx()/CLHEP::m;
00456         y_el_p       = (*hc)[i]->Gety()/CLHEP::m;
00457         z_el_p       = (*hc)[i]->Getz()/CLHEP::m;
00458         part_el_p    = (*hc)[i]->GetPartID();
00459         turnnumber   = (*hc)[i]->GetTurnsTaken();
00460         eventno_el_p = (*hc)[i]->GetEventNo();
00461         //name - convert to char array for root
00462         G4String temp = (*hc)[i]->GetName()+'\0';
00463         strncpy(volumeName_el_p,temp.c_str(),sizeof(volumeName_el_p)-1);
00464         PrecisionRegionEnergyLossTree->Fill();
00465       }
00466     }
00467 }
00468 
00469 void BDSOutputROOT::WritePrimaryLoss(BDSEnergyCounterHit* hit)
00470 {
00471   //copy variables from hit to root variables
00472   X_pl          = hit->GetX()/CLHEP::m;
00473   Y_pl          = hit->GetY()/CLHEP::m;
00474   Z_pl          = hit->GetZ()/CLHEP::m;
00475   S_pl          = hit->GetS()/CLHEP::m;
00476   x_pl          = hit->Getx()/CLHEP::m;
00477   y_pl          = hit->Gety()/CLHEP::m;
00478   z_pl          = hit->Getz()/CLHEP::m;
00479   E_pl          = hit->GetEnergy()/CLHEP::GeV;
00480   weight_pl     = hit->GetWeight();
00481   part_pl       = hit->GetPartID();
00482   turnnumber_pl = hit->GetTurnsTaken();
00483   eventno_pl    = hit->GetEventNo();
00484   
00485   //write to file
00486   PrimaryLossTree->Fill();
00487 }
00488 
00489 void BDSOutputROOT::WritePrimaryHit(BDSEnergyCounterHit* hit)
00490 {
00491   //copy variables from hit to root variables
00492   X_ph          = hit->GetX()/CLHEP::m;
00493   Y_ph          = hit->GetY()/CLHEP::m;
00494   Z_ph          = hit->GetZ()/CLHEP::m;
00495   S_ph          = hit->GetS()/CLHEP::m;
00496   x_ph          = hit->Getx()/CLHEP::m;
00497   y_ph          = hit->Gety()/CLHEP::m;
00498   z_ph          = hit->Getz()/CLHEP::m;
00499   E_ph          = hit->GetEnergy()/CLHEP::GeV;
00500   weight_ph     = hit->GetWeight();
00501   part_ph       = hit->GetPartID();
00502   turnnumber_ph = hit->GetTurnsTaken();
00503   eventno_ph    = hit->GetEventNo();
00504   
00505   //write to file
00506   PrimaryHitsTree->Fill();
00507 }
00508 
00509 void BDSOutputROOT::WriteHistogram(BDSHistogram1D* hIn)
00510 {
00511   G4String hname = hIn->GetName();
00512   hname = BDS::PrepareSafeName(hname);
00513 
00514   std::vector<G4double> binLowerEdges = hIn->GetBinLowerEdges();
00515   if (binLowerEdges.size() < 1)
00516     {return;} //no bins - don't write the histogram
00517   binLowerEdges.push_back(hIn->GetLastBin()->GetUpperEdge()); //root requires nbins+1
00518 
00519   //always construct the histogram by bin edges - works both with constant
00520   //and variable bin width histograms
00521   // &vector[0] gives an array to the contents of the vector - ensured as
00522   // standard is that the vector's contents are contiguous
00523   TH1D* h = new TH1D(hname, hIn->GetTitle(), hIn->GetNBins(), &binLowerEdges[0]);
00524 
00525   G4int i;
00526   for(hIn->first(),i = 1;!hIn->isDone();hIn->next(), i++)
00527     {
00528       BDSBin* currentBin = hIn->currentBin();
00529       h->SetBinContent(i,currentBin->GetValue());
00530       h->SetBinError(i,  currentBin->GetError());
00531     }
00532   //over / underflow manually set
00533   h->SetBinContent(0,hIn->GetUnderflowBin()->GetValue()); //underflow
00534   h->SetBinContent(0,hIn->GetUnderflowBin()->GetError());
00535   h->SetBinContent(h->GetNbinsX()+1,hIn->GetOverflowBin()->GetValue()); //overflow
00536   h->SetBinContent(h->GetNbinsX()+1,hIn->GetOverflowBin()->GetError());
00537 
00538   h->SetEntries(hIn->GetNEntries());
00539   
00540   h->Write(); // as commit actually closes a file as does write..
00541   delete h;
00542 }
00543 
00544 void BDSOutputROOT::Commit()
00545 {
00546   Write();
00547   outputFileNumber++;
00548   Init();
00549 }
00550 
00551 void BDSOutputROOT::Write()
00552 {
00553 #ifdef BDSDEBUG
00554   G4cout << __METHOD_NAME__ << G4endl;
00555 #endif
00556 
00557   if(theRootOutputFile && theRootOutputFile->IsOpen())
00558     {
00559 #ifdef BDSDEBUG
00560       G4cout << __METHOD_NAME__ << " - ROOT file found and open, writing." << G4endl;
00561 #endif
00562       //Dump all other quantities to file...
00563       theRootOutputFile->Write();
00564       theRootOutputFile->Close();
00565       delete theRootOutputFile;
00566       theRootOutputFile=NULL;
00567     }
00568   G4cout << __METHOD_NAME__ << " ...finished." << G4endl;
00569 }
00570 #endif
Generated on 28 Jun 2015 for BDSIM by  doxygen 1.4.7