// BLCMDntuple.cc /* This source file is part of G4beamline, http://g4beamline.muonsinc.com Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. http://www.gnu.org/copyleft/gpl.html */ #include #include #include #include #include "BLAssert.hh" #include "BLNTuple.hh" #include "BLCommand.hh" #include "BLManager.hh" #include "BLRunManager.hh" typedef unsigned long long Bitmask; const Bitmask Bit = 1; const int MAX_DETECTORS = sizeof(Bitmask)*8; /** class BLCMDntuple - implements an ntuple. * * An ntuple holds the data from multiple virtualdetector-s in a single * NTuple, with one row per event. If multiple tracks hit a given * virtual detector in the ntuple, only the first one is kept. **/ class BLCMDntuple : public BLCommand, public BLManager::EventAction, BLManager::TrackingAction, BLNTupleCallback { G4String name; G4String category; G4String detectors; // comma-separated list: detector patterns G4String required; // comma-separated list: required det patterns G4String veto; // comma-separated list: veto det patterns G4String format; G4String filename; G4int minHit; G4int perTrack; G4int unionFlag; // internal data: G4String fields; // colon-separated list: fields BLNTuple *ntuple; double *data; int ndata; std::map place; // map from caller to index std::vector start; // index in data[] std::set vetoSet; // set of veto NTuples Bitmask need; Bitmask got; int nHit; int vetoed; public: /// Default constructor. BLCMDntuple(); /// Destructor. virtual ~BLCMDntuple() { } /// Copy constructor. BLCMDntuple(const BLCMDntuple& r); /// commandName() returns "ntuple". G4String commandName() { return "ntuple"; } /// command() implements the ntuple command. int command(BLArgumentVector& argv, BLArgumentMap& namedArgs); /// defineNamedArgs() defines the named arguments for the command. void defineNamedArgs(); /// setupFirstEvent() will setup the NTuple (called at first event) /// Done here to be sure all other NTUples have been created (e.g. in /// callbacks). void setupFirstEvent(); /// from BLManager::EventAction void BeginOfEventAction(const G4Event* event); void EndOfEventAction(const G4Event* event); /// from BLManager::TrackingAction virtual void PreUserTrackingAction(const G4Track *track); virtual void PostUserTrackingAction(const G4Track *track); /// from BLNTupleCallback virtual void ntupleCallback(BLNTuple *caller, double data[], int ndata); /// returns true if the string s matches one of the patterns in pattern. static bool matches(G4String s, G4String pattern); }; BLCMDntuple defaultNtuple; BLCMDntuple::BLCMDntuple() : BLCommand(), BLManager::EventAction(), place(), start(), vetoSet() { registerCommand(BLCMDTYPE_DATA); setSynopsis("Define an NTuple containing multiple detectors."); setDescription("An ntuple holds the data from multiple " "detectors in a single NTuple, with one row (entry) per " "event or per track. This permits the generation of plots " "that compare different detectors. " "Up to 64 detectors can be used. " "While 'detector' is used in this description, any existing " "NTuple can be used, generated by any command, such as: " "virtualdetector, zntuple, beamlossntuple, timentuple, and " "newparticlentuple.\n\n" "There are two ways the detectors can be combined: the default " "method is to construct a new NTuple that combines the fields " "of all the detectors, each prepended with the detector name. " "If union=1 is given, then all detectors must have the same " "list of fields, and that list is used for this NTuple; any " "hit in any detector is simply copied to this NTuple -- this " "permits multiple detectors to be combined into a single " "NTuple (you may want noSingles=1 for the detectors). With " "union=1, the 'required', 'veto' and 'minHit' arguments cannot " "be used.\n\n" "By default, an entry in this ntuple is made for each event " "satisfying the require, veto, and minHit conditions; if " "perTrack=1 then an entry is made for each track that " "satisfies the require, veto, and minHit conditions. " "Any hit in any detector matching the patterns in veto will " "prevent the event/track from being entered into the ntuple." "\n\n" "If multiple hits occur in a given detector during the " "event or track, only the first one is kept in this ntuple. " "Detectors are specified by patterns identical to " "UNIX file-matching, so '*Det*' matches any detector " "with 'Det' anywhere in its name, etc. The 'required' argument " "permits events to be omitted unless all of the matching " "detectors were hit at least once in the event or track. The " "patterns in 'required' are applied only to detectors " "in the ntuple, so a simple '*' only matches detectors named " "in the 'detectors' argument.\n\n" "If union is nonzero, the hits in each detector are entered " "into this NTuple as they occur; each hit in any detector " "is included as a row in this NTuple. All detectors " "must have the same fields, which are used for this NTuple.\n\n" "NOTE: the name of a detector is by default the " "concatenation of its ancestors' names before its own (except " "World), unless " "rename=NAME was used in its place command. The patterns are " "applied to the names of the virtualdetectors as they were " "placed (including rename), not the bare name of the " "virtualdetector command. If 'rename=det#' was used when " "placing " "the virtualdetector-s, you probably want a * to match the #, " "or list them individually (det1,det2,det3...).\n\n" "NOTE: This command does not work correctly in collective " "tracking mode, unless union=1.\n\n" "Valid Formats (ignore case): "); description += BLNTuple::getFormatList(); // initial default values: name = ""; category = "NTuple"; detectors = ""; required = "*"; veto = ""; format = ""; filename = ""; minHit = 0; perTrack = 0; unionFlag = 0; fields = ""; ntuple = 0; data = 0; ndata = 0; need = 0; got = 0; nHit = 0; vetoed = 0; } BLCMDntuple::BLCMDntuple(const BLCMDntuple &r) : BLCommand(r), BLManager::EventAction(), BLManager::TrackingAction(), BLNTupleCallback(), place(), start(), vetoSet() { name = ""; category = r.category; detectors = r.detectors; required = r.required; veto = r.veto; format = r.format; filename = r.filename; minHit = r.minHit; perTrack = r.perTrack; unionFlag = r.unionFlag; fields = r.fields; ntuple = 0; data = 0; ndata = 0; need = r.need; got = 0; nHit = 0; vetoed = 0; } int BLCMDntuple::command(BLArgumentVector& argv, BLArgumentMap& namedArgs) { if(argv.size() != 1) { printError("Invalid ntuple command -- need name"); return -1; } if(argv[0] == "default") { handleNamedArgs(namedArgs); return -1; } BLCMDntuple *t = new BLCMDntuple(defaultNtuple); int retval = t->handleNamedArgs(namedArgs); t->name = argv[0]; if(t->unionFlag) { if(t->required != "*" || t->veto.size() > 0 || t->minHit > 1) printError("ntuple: In union mode, 'required', 'veto'," " and 'minHit' cannot be used."); } t->print(argv[0]); if(t->perTrack == 0) BLManager::getObject()->registerEventAction(t,false); else // at front, so newparticlentuple will work BLManager::getObject()->registerTrackingAction(t,true); return retval; } void BLCMDntuple::defineNamedArgs() { argString(category,"category","The category of the NTuple."); argString(detectors,"detectors","A comma-separated list of detector patterns."); argString(required,"required","A comma-separated list of required detector patterns(default=*)."); argString(veto,"veto","A comma-separated list of detector patterns, any hit cancels entry into the NTuple (default='')."); argString(format,"format","NTuple format (see above for list)."); argString(filename,"filename","Name of file."); argInt(minHit,"minHit","Minimum number of detectors hit (default 0)."); argInt(perTrack,"perTrack","Nonzero for an entry per track; " "0 for an entry per event (default 0)."); argInt(unionFlag,"union","Set nonzero to perform a union of the " "detectors, rather than combining them."); argString(required,"require","Synonym for required."); argInt(minHit,"minHits","Synonym for minHit."); argString(filename,"file","Synonym for filename."); } void BLCMDntuple::setupFirstEvent() { std::vector list = BLNTuple::getList(); for(unsigned i=0; igetName(); if(matches(detName,detectors)) { G4String detFields = detNtuple->getFields(); detNtuple->registerCallback(this); int j = start.size(); if(unionFlag) { if(fields.size() == 0) { fields = detFields; ndata = detNtuple->getNData(); } if(fields != detFields) { G4Exception("ntuple command", "Field mismatch", FatalException, "Detectors have different fields"); } continue; } place[detNtuple] = j; start.push_back(ndata); if(matches(detName,required)) { need |= (Bit< 0 && matches(detName,veto)) { detNtuple->registerCallback(this); vetoSet.insert(detNtuple); } } if(ndata <= 0) { printError("ntuple: no detectors in ntuple '%s'! " "Possibilities are:",name.c_str()); for(unsigned i=0; igetName(); printf("%s ",detName.c_str()); } printf("\n"); data = new double[1]; need = ~(Bitmask)0; return; } data = new double[ndata]; ntuple = BLNTuple::create(format,category,name,fields,filename); } void BLCMDntuple::BeginOfEventAction(const G4Event* event) { if(!data) { setupFirstEvent(); if(unionFlag) return; if(BLRunManager::getObject()->getCollectiveMode()) G4Exception("ntuple command","Collective mode error", FatalException, "Does not work in collective tracking mode"); } got = (Bitmask)0; nHit = 0; vetoed = 0; for(int i=0; i= minHit && !vetoed) ntuple->appendRow(data,ndata); } void BLCMDntuple::PreUserTrackingAction(const G4Track *track) { BeginOfEventAction(0); } void BLCMDntuple::PostUserTrackingAction(const G4Track *track) { EndOfEventAction(0); } void BLCMDntuple::ntupleCallback(BLNTuple *caller,double *detData,int ndetData) { if(unionFlag) { ntuple->appendRow(detData,ndetData); return; } if(place.count(caller) > 0) { unsigned j = place[caller]; BLAssert(j < start.size()); Bitmask me = (Bit<= 0 && k < ndata); data[k++] = detData[l]; } } } if(vetoSet.count(caller) > 0) ++vetoed; } bool BLCMDntuple::matches(G4String s, G4String pattern) { G4String::size_type place = 0; while(place != pattern.npos) { G4String::size_type next = pattern.find(",",place); G4String word = pattern.substr(place, (next==pattern.npos ? pattern.npos : next-place)); if(fnmatch(word.c_str(),s.c_str(),0) == 0) return true; place = next; if(next != pattern.npos) place = next + 1; } return false; }