00001
00002
00003
00004
00005
00006
00007
00008 #include "MuonRemoval/MCAnalysis.h"
00009 #include "MessageService/MsgService.h"
00010 #include "MinosObjectMap/MomNavigator.h"
00011 #include "JobControl/JobCModuleRegistry.h"
00012 #include "CandDigit/CandDigitHandle.h"
00013 #include "CandDigit/CandDeMuxDigitHandle.h"
00014 #include "CandDigit/CandDigit.h"
00015 #include "CandDigit/CandDigitListHandle.h"
00016 #include "CandDigit/CandDeMuxDigitListHandle.h"
00017
00018 #include "RecoBase/CandTrackListHandle.h"
00019 #include "RecoBase/CandTrackHandle.h"
00020 #include "RecoBase/CandEventListHandle.h"
00021 #include "RecoBase/CandEventHandle.h"
00022 #include "RecoBase/CandStripHandle.h"
00023 #include "CandData/CandRecord.h"
00024 #include "CandData/CandHeader.h"
00025
00026 #include "DataUtil/Truthifier.h"
00027 #include "DataUtil/TruthHelper.h"
00028 #include "Record/SimSnarlRecord.h"
00029 #include "REROOT_Classes/REROOT_NeuKin.h"
00030
00031 #include "TClonesArray.h"
00032
00033 #include <cmath>
00034 #include <set>
00035 using std::set;
00036
00037
00038 JOBMODULE(MCAnalysis, "MCAnalysis",
00039 "");
00040 CVSID("$Id: MCAnalysis.cxx,v 1.4 2007/02/01 22:10:14 rhatcher Exp $");
00041
00042
00043 MCAnalysis::MCAnalysis() :
00044 fNameListIn(""),
00045
00046 fMakePictures(0)
00047 {
00051 }
00052
00053
00054 MCAnalysis::~MCAnalysis()
00055 {
00059 }
00060
00061
00062
00063 void MCAnalysis::BeginJob()
00064 {
00068 fCanvas = new TCanvas("mcana", "MC Analysis", 800, 600);
00069 fUZAxis = new TH2F("uzaxis", "", 500, 0, 500, 192, 0, 192);
00070 fVZAxis = new TH2F("vzaxis", "", 500, 0, 500, 192, 0, 192);
00071 fCanvas2 = new TCanvas("mcana2", "MC Analysis (Q)", 800, 600);
00072 fUZTrkQ = new TH1F("trkquz", "", 500, 0, 500);
00073 fVZTrkQ = new TH1F("trkqvz", "", 500, 0, 500);
00074
00075 TDirectory* cdir = gDirectory;
00076 fFileOut = TFile::Open("mc_performance.root", "recreate");
00077 fFileOut->cd();
00078 fTreeOut = new TTree("mcp", "MC performance");
00079 fTreeOut->SetDirectory(fFileOut);
00080 fTreeOut->Branch("run", &fRun, "run/I");
00081 fTreeOut->Branch("snarl", &fSnarl, "snarl/I");
00082 fTreeOut->Branch("nu", &fNu, "nu/I");
00083 fTreeOut->Branch("action", &fAction, "action/I");
00084 fTreeOut->Branch("enu", &fENu, "enu/F");
00085 fTreeOut->Branch("emu", &fEMu, "mnu/F");
00086
00087 fTreeOut->Branch("nret", &fNRetained, "nret/I");
00088 fTreeOut->Branch("nretmu", &fNRetainedMuon, "nretmu/I");
00089 fTreeOut->Branch("nretshw", &fNRetainedShw, "nretshw/I");
00090 fTreeOut->Branch("nretboth", &fNRetainedBoth, "nretboth/I");
00091
00092 fTreeOut->Branch("peret", &fPERetained, "peret/F");
00093 fTreeOut->Branch("peretmu", &fPERetainedMuon, "peretmu/F");
00094 fTreeOut->Branch("peretshw", &fPERetainedShw, "peretshw/F");
00095 fTreeOut->Branch("peretboth", &fPERetainedBoth, "peretboth/F");
00096
00097 fTreeOut->Branch("nrej", &fNRejected, "nrej/I");
00098 fTreeOut->Branch("nrejmu", &fNRejectedMuon, "nrejmu/I");
00099 fTreeOut->Branch("nrejshw", &fNRejectedShw, "nrejshw/I");
00100 fTreeOut->Branch("nrejboth", &fNRejectedBoth, "nrejboth/I");
00101
00102 fTreeOut->Branch("nmuondig", &fNMuonDig, "nmuondig/I");
00103 fTreeOut->Branch("nmuondigret", &fNMuonDigRetained, "nmuondigret/I");
00104
00105 fTreeOut->Branch("nshwdig", &fNShwDig, "nshwdig/I");
00106 fTreeOut->Branch("nshwdigret", &fNShwDigRetained, "nshwdigret/I");
00107
00108 fTreeOut->Branch("nshwdigatvtx", &fNShwDigAtVtx, "nshwdigatvtx/I");
00109 fTreeOut->Branch("nshwdigretatvtx", &fNShwDigRetainedAtVtx, "nshwdigretatvtx/I");
00110
00111 fTreeOut->Branch("nshwpe", &fNShwPE, "nshwpe/F");
00112 fTreeOut->Branch("nshwperet", &fNShwPERetained, "nshwpepe/F");
00113
00114 fTreeOut->Branch("nshwpeatvtx", &fNShwPEAtVtx, "nshwpeatvtx/F");
00115 fTreeOut->Branch("nshwperetatvtx", &fNShwPERetainedAtVtx, "nshwpepeatvtx/F");
00116
00117 fTreeOut->Branch("rshwn", &fRShwN, "rshwn/I");
00118 fTreeOut->Branch("rshwe", fRShwE, "rshwe[rshwn]/F");
00119 fTreeOut->Branch("rshwdist", fRShwDist, "rshwdist[rshwn]/I");
00120 fTreeOut->Branch("rshwmuon", fRShwMuon, "rshwmuon[rshwn]/I");
00121 fTreeOut->Branch("rshwtrk", fRShwTrk, "rshwtrk[rshwn]/I");
00122 fTreeOut->Branch("rshwxtalk", fRShwXTalk, "rshwxtalk[rshwn]/I");
00123
00124 fTreeOut->Branch("rshwhittype", fRShwHitType, "rshwhittype[rshwn]/I");
00125 fTreeOut->Branch("rshwhitstatus", fRShwHitStatus, "rshwhitstatus[rshwn]/I");
00126 fTreeOut->Branch("rshwhitmother", fRShwHitMother, "rshwhitmoher[rshwn]/I");
00127
00128 fTreeOut->Branch("fTrkMIPCalibOkay", &fTrkMIPCalibOkay, "fTrkMIPCalibOkay/I");
00129
00130 hMIPTrueMuon = new TH1F("hMIPTrueMuon", "", 100, 0, 10);
00131 hMIPMissTrack = new TH1F("hMIPMissTrack", "", 100, 0, 10);
00132 hMIPMixed = new TH1F("hMIPMixed", "", 100, 0, 10);
00133 hMIPTrueMuon2 = new TH1F("hMIPTrueMuon2", "", 100, 0, 10);
00134 hMIPMissTrack2 = new TH1F("hMIPMissTrack2", "", 100, 0, 10);
00135 hMIPMixed2 = new TH1F("hMIPMixed2", "", 100, 0, 10);
00136
00137 cdir->cd();
00138 }
00139
00140
00141
00142 void MCAnalysis::EndJob()
00143 {
00144 TDirectory* cdir = gDirectory;
00145 fFileOut->cd();
00146 fTreeOut->SetDirectory(fFileOut);
00147 fTreeOut->Write();
00148 hMIPTrueMuon ->SetDirectory(fFileOut);
00149 hMIPMissTrack->SetDirectory(fFileOut);
00150 hMIPMixed ->SetDirectory(fFileOut);
00151 hMIPTrueMuon2->SetDirectory(fFileOut);
00152 hMIPMissTrack2->SetDirectory(fFileOut);
00153 hMIPMixed2 ->SetDirectory(fFileOut);
00154
00155 hMIPTrueMuon ->Write();
00156 hMIPMissTrack->Write();
00157 hMIPMixed ->Write();
00158 hMIPTrueMuon2 ->Write();
00159 hMIPMissTrack2->Write();
00160 hMIPMixed2 ->Write();
00161
00162 cdir->cd();
00163
00164 }
00165
00166
00167
00168 JobCResult MCAnalysis::Ana(const MomNavigator* mom)
00169 {
00170
00171
00172
00173 CandRecord* record = dynamic_cast<CandRecord*>(mom->GetFragment("CandRecord"));
00174 if(!record){
00175 MSG("MCAna",Msg::kError) << " Unable to find CandRecord in mom !!! " << endl;
00176 return JobCResult::kFailed;
00177 }
00178 const VldContext &vldc = *(record->GetVldContext());
00179 if(vldc.GetSimFlag()!=SimFlag::kMC){
00180 MSG("MCAna",Msg::kWarning) << " Event not a MC event"<<endl;
00181 return JobCResult::kPassed;
00182 }
00183 Truthifier* truth = new Truthifier(mom);
00184 TruthHelper* beer = new TruthHelper(mom);
00185 if(!(truth->IsValid())){
00186 MSG("MCAna",Msg::kError) << " The THRUTH is not valid! "<<endl;
00187 return JobCResult::kFailed;
00188 }
00189
00190
00191
00192 const CandEventListHandle * eventlist = dynamic_cast<const CandEventListHandle*>(record->FindCandHandle("CandEventListHandle"));
00193 if(!eventlist || eventlist->GetNDaughters()!=1){
00194 MSG("MCAna",Msg::kWarning) << " Rejecting event as it has no events " << endl;
00195 return JobCResult::kFailed;
00196 }
00197
00198
00199 TIter event_iter(eventlist->GetDaughterIterator());
00200 const CandEventHandle* event = dynamic_cast<const CandEventHandle*>(event_iter());
00201 assert(event);
00202
00203
00204 const CandTrackHandle* track = dynamic_cast<const CandTrackHandle*>(event->GetPrimaryTrack());
00205 assert(track);
00206
00207 const CandDigitListHandle* rawdigitlist = dynamic_cast<const CandDigitListHandle*>(record->FindCandHandle("CandDigitListHandle", "altdemux"));
00208 const CandDigitListHandle* strippeddigitlist = dynamic_cast<const CandDigitListHandle*>(record->FindCandHandle("CandDigitListHandle", "stripdigitlist"));
00209
00210 if(rawdigitlist==NULL || strippeddigitlist==NULL){
00211 MSG("MCAna",Msg::kError) << " Unable to find Merged digit list (" << strippeddigitlist<<") or RawDigitList (" << rawdigitlist <<")"<< endl;
00212 return JobCResult::kFailed;
00213 }
00214
00215 const CandHeader* header = dynamic_cast<const CandHeader*>(record->GetHeader());
00216
00217
00218
00219 fUHits.clear();
00220 fVHits.clear();
00221 fUZTrkQ->Reset();
00222 fVZTrkQ->Reset();
00223
00224
00225
00226
00227 SimSnarlRecord *simsnarl = dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
00228 if (!simsnarl){
00229 MSG("MCAna", Msg::kError)<< " Faied to get the truth (SimSnarl) " <<endl;
00230 return JobCResult::kFailed;
00231 }
00232 const TClonesArray* ctca = dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","StdHep"));
00233 const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","NeuKinList"));
00234 const REROOT_NeuKin* neukin = dynamic_cast<const REROOT_NeuKin*>(neukinarray->At(0));
00235
00236 assert(neukin);
00237 const bool numucc = (abs(neukin->INu())==14 && abs(neukin->IAction()!=0) );
00238
00239
00240
00241
00242 set<Int_t> retainedDigitIds;
00243 TIter strippeddigititer(strippeddigitlist->GetDaughterIterator());
00244 while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(strippeddigititer())) {
00245 if(retainedDigitIds.count(digit->GetRawDigitIndex())!=0){
00246 MSG("MCAna", Msg::kError)<< " Mulitple RAWDIGITID entries " <<endl;
00247 }else{
00248 if(digit->GetQieErrorBits()==0){
00249 retainedDigitIds.insert(digit->GetRawDigitIndex());
00250 }else{
00251 retainedDigitIds.insert(-1*digit->GetRawDigitIndex());
00252 }
00253 }
00254 }
00255
00256
00257
00258
00259 fRun = header->GetRun();
00260 fSnarl = header->GetSnarl();
00261 fNu = neukin->INu();
00262 fAction = neukin->IAction();
00263 fENu = fabs((neukin->P4Neu())[3]);
00264 fEMu = fabs((neukin->P4Mu1())[3]);
00265
00266 fNMuonDig = 0;
00267 fNMuonDigRetained = 0;
00268
00269 fNShwDig = 0;
00270 fNShwDigRetained = 0;
00271
00272 fNShwDigAtVtx = 0;
00273 fNShwDigRetainedAtVtx = 0;
00274
00275 fNShwPE = 0 ;
00276 fNShwPERetained =0;
00277
00278 fNShwPEAtVtx = 0 ;
00279 fNShwPERetainedAtVtx =0;
00280
00281 fNRetained = 0 ;
00282 fNRetainedMuon = 0;
00283 fNRetainedShw = 0;
00284 fNRetainedBoth = 0;
00285
00286 fPERetained = 0 ;
00287 fPERetainedMuon = 0;
00288 fPERetainedShw = 0;
00289 fPERetainedBoth = 0;
00290
00291 fNRejected = 0;
00292 fNRejectedMuon = 0;
00293 fNRejectedShw = 0;
00294 fNRejectedBoth = 0;
00295
00296 fRShwN = 0 ;
00297 for(int i=0; i<500; i++){
00298 fRShwE[i] = 0;
00299 fRShwDist[i] = 0;
00300 fRShwMuon[i] = 0;
00301 fRShwTrk[i] = 0;
00302 fRShwHitType[i] = 0;
00303 fRShwHitStatus[i] = 0;
00304 fRShwHitMother[i] = 0;
00305 }
00306
00307
00308
00309
00310 int lNRejShw = 0;
00311 int lNRejShwMaxTrk = 0 ;
00312 int lNRejShwFakeTrk = 0;
00313 int lNRejShwMix = 0 ;
00314
00315
00316
00317
00318 FillPlnInfo((CandDeMuxDigitListHandle*)rawdigitlist);
00319 FillTrkInfo(track);
00320
00321 const int maxtrkplane = GetMaxTrkPlane();
00322
00323 int minplane = 500;
00324 int maxplane = 0;
00325
00326
00327 bool plane_track[500];
00328 bool plane_muon[500];
00329 bool plane_shw[500];
00330 for(int i =0; i<500; i++){
00331 plane_track[i] = plane_muon[i] = plane_shw[i] = false;
00332 }
00333
00334 TIter rawdigititer(rawdigitlist->GetDaughterIterator());
00335
00336 while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(rawdigititer())) {
00337
00338 const int planeno = digit->GetPlexSEIdAltL().GetBestSEId().GetPlane();
00339 const int stripno = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
00340 const CandDigitHandle tcdh = *digit;
00341 const CandDeMuxDigitHandle* demuxdigit = dynamic_cast<const CandDeMuxDigitHandle*>(digit);
00342 assert(demuxdigit);
00343 bool is_shower = 0;
00344 bool is_muon = 0;
00345 bool is_physics = 0;
00346
00347 const bool is_retained = (retainedDigitIds.count(digit->GetRawDigitIndex())>0);
00348 bool is_track = (fTrkStrips.count(planeno*192+stripno)>0);
00349 const bool is_recoed_xtalk = (demuxdigit->GetDeMuxDigitFlagWord()&CandDeMuxDigit::kXTalk);
00350 const bool is_retained_scaled = (retainedDigitIds.count(-1*digit->GetRawDigitIndex())>0);
00351
00352 const DigiSignal * signal = truth->GetSignal(tcdh);
00353
00354 const DigiScintHit* biggest_hit = truth->GetBiggestHit(tcdh);
00355 TParticle* biggest_part = NULL;
00356 if(biggest_hit!=NULL){
00357 Int_t biggest_track = abs(biggest_hit->TrackId());
00358 biggest_part = dynamic_cast<TParticle*>((*ctca)[biggest_track]);
00359 }
00360
00361 if(signal!=NULL){
00362 is_physics = ((signal->GetTruth() & DigiSignal::kGenuine)!=0);
00363 for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
00364 const DigiScintHit * hit = signal->GetHit(i);
00365 if(hit){
00366 Int_t track = abs(hit->TrackId());
00367 TParticle* part = dynamic_cast<TParticle*>((*ctca)[track]);
00368 if(!beer->IsNeutrino(part) && !beer->IsMuon(part)){
00369 is_shower=true;
00370 }
00371 if(!beer->IsNeutrino(part) && beer->IsMuon(part)){
00372 is_muon=true;
00373 }
00374 }
00375 }
00376 }
00377
00378
00379 if(is_physics){
00380 if(is_shower){
00381 fNShwDig++;
00382 fNShwPE+=digit->GetCharge(CalDigitType::kPE);
00383 if(planeno<maxtrkplane){
00384 fNShwDigAtVtx++;
00385 fNShwPEAtVtx+=digit->GetCharge(CalDigitType::kPE);
00386 }
00387 if(is_retained || is_retained_scaled){
00388 fNShwDigRetained++;
00389 fNShwPERetained+=digit->GetCharge(CalDigitType::kPE);
00390 if(planeno<maxtrkplane){
00391 fNShwDigRetainedAtVtx++;
00392 fNShwPERetainedAtVtx+=digit->GetCharge(CalDigitType::kPE);
00393 }
00394 }else{
00395 if(fRShwN<1000){
00396 fRShwE[fRShwN]=digit->GetCharge(CalDigitType::kPE);
00397 fRShwDist[fRShwN]=maxtrkplane - planeno;
00398 fRShwMuon[fRShwN]=(int)is_muon;
00399 fRShwTrk[fRShwN]=(int)(is_track && !is_muon);
00400 fRShwXTalk[fRShwN]= (int)(is_recoed_xtalk);
00401 if(biggest_part!=NULL){
00402 fRShwHitType[fRShwN] = biggest_part->GetPdgCode();
00403 fRShwHitStatus[fRShwN] = biggest_part->GetStatusCode();
00404 fRShwHitMother[fRShwN] = biggest_part->GetMother(0);
00405 }
00406 fRShwN++;
00407 }
00408 lNRejShw++;
00409 if(planeno>maxtrkplane) lNRejShwMaxTrk++;
00410 if(is_muon) lNRejShwMix++;
00411 if(!is_muon && is_track) lNRejShwFakeTrk++;
00412 }
00413 }
00414 if(is_muon){
00415 fNMuonDig++;
00416 if(is_retained || is_retained_scaled) fNMuonDigRetained++;
00417 }
00418 if(is_retained || is_retained_scaled){
00419 fNRetained++;
00420 fPERetained+=digit->GetCharge(CalDigitType::kPE);
00421 if(is_muon){
00422 fNRetainedMuon++;
00423 fPERetainedMuon+=digit->GetCharge(CalDigitType::kPE);
00424 }
00425 if(is_shower){
00426 fNRetainedShw++;
00427 fPERetainedShw+=digit->GetCharge(CalDigitType::kPE);
00428 }
00429 if(is_muon && is_shower){
00430 fNRetainedBoth++;
00431 fPERetainedBoth+=digit->GetCharge(CalDigitType::kPE);
00432 }
00433 }else{
00434 fNRejected++;
00435 if(is_muon && !is_shower)fNRejectedMuon++;
00436 if(!is_muon && is_shower) fNRejectedShw++;
00437 if(is_muon && is_shower) fNRejectedBoth++;
00438 }
00439
00440
00441
00442 if(is_track){
00443 plane_track[planeno] = true;
00444 if(is_muon) plane_muon[planeno] = true;
00445 if(is_shower) plane_shw[planeno] = true;
00446 }
00447
00448
00449
00450
00451 TMarker marker(planeno, stripno,25);
00452 if( is_muon && !is_shower) marker.SetMarkerStyle(28);
00453 if(!is_muon && is_shower) marker.SetMarkerStyle(24);
00454 if( is_muon && is_shower) marker.SetMarkerStyle(26);
00455 if(is_retained)marker.SetMarkerColor(1);
00456 else if(is_retained_scaled)marker.SetMarkerColor(4);
00457 else marker.SetMarkerColor(2);
00458 if(digit->GetPlexSEIdAltL().GetBestSEId().GetPlaneView()==PlaneView::kU){
00459 fUHits.push_back(marker);
00460 }else{
00461 fVHits.push_back(marker);
00462 }
00463 if(is_physics && is_shower){
00464 if(planeno>maxplane) maxplane = planeno;
00465 if(planeno<minplane) minplane = planeno;
00466 }
00467 }
00468
00469 }
00470
00471
00472
00473 int ntrkcalibokay = 0;
00474 for(int i=0; i<500; i++){
00475 if(plane_track[i]){
00476 if(fTrkMIPDcosz[i]!=0){
00477 ntrkcalibokay++;
00478 if(plane_muon[i] && !plane_shw[i]) hMIPTrueMuon ->Fill(fTrkMIPDcosz[i]);
00479 if(plane_muon[i] && plane_shw[i]) hMIPMixed ->Fill(fTrkMIPDcosz[i]);
00480 if(!plane_muon[i]) hMIPMissTrack->Fill(fTrkMIPDcosz[i]);
00481 if(i<maxtrkplane){
00482 if(plane_muon[i] && !plane_shw[i]) hMIPTrueMuon2 ->Fill(fTrkMIPDcosz[i]);
00483 if(plane_muon[i] && plane_shw[i]) hMIPMixed2 ->Fill(fTrkMIPDcosz[i]);
00484 if(!plane_muon[i]) hMIPMissTrack2->Fill(fTrkMIPDcosz[i]);
00485 }
00486 }
00487 }
00488 }
00489 fTrkMIPCalibOkay = (ntrkcalibokay>4);
00490 fTreeOut->Fill();
00491
00492 if(fMakePictures){
00493 cout << " Event MuCC: " << ((numucc)?"YES":"NO") <<endl;
00494 cout << " Max Trk Pln: "<< maxtrkplane<<endl;
00495 cout << " Muon Digits : " << fNMuonDig<<endl;
00496 cout << " Retained : " << fNMuonDigRetained / (float)fNMuonDig << " (" << fNMuonDigRetained << ")"<<endl;
00497 cout << " Rejected : " << (fNMuonDig - fNMuonDigRetained) / (float)fNMuonDig << " (" << fNMuonDig - fNMuonDigRetained << ")"<<endl;
00498 cout << " Shower Digits : " << fNShwDig <<endl;
00499 cout << " Retained : " << fNShwDigRetained / (float)fNShwDig << " (" << fNShwDigRetained <<")" <<endl;
00500 cout << " Rejected : " << (fNShwDig - fNShwDigRetained) / (float)fNShwDig << " (" << fNShwDig - fNShwDigRetained <<")" <<endl;
00501
00502
00503 cout << " Rejected : " << fNRejected<<endl;
00504 cout << " Muon only: " << fNRejectedMuon/(float)fNRejected << " " << fNRejectedMuon<<endl;
00505 cout << " Shw only : " << fNRejectedShw/(float)fNRejected << " " << fNRejectedShw<<endl;
00506 cout << " both : " << fNRejectedBoth/(float)fNRejected << " " << fNRejectedBoth<<endl;
00507 cout << " Retained : " << fNRetained<<endl;
00508 cout << " Muon only: " << fNRetainedMuon/(float)fNRetained << " " << fNRetainedMuon<<endl;
00509 cout << " Shw only : " << fNRetainedShw/(float)fNRetained << " " << fNRetainedShw<<endl;
00510 cout << " both : " << fNRetainedBoth/(float)fNRetained << " " << fNRetainedBoth<<endl;
00511
00512 cout << " Rejected Shw: " <<endl;
00513 cout << " Total: " << lNRejShw <<endl;
00514 if(lNRejShw>0){
00515 cout << " MAXTRK: " << lNRejShwMaxTrk << " " << lNRejShwMaxTrk / lNRejShw<<endl;
00516 cout << " MIXED: " << lNRejShwMix << " " << lNRejShwMix / lNRejShw<<endl;
00517 cout << " FAKETRK: " << lNRejShwFakeTrk << " " << lNRejShwFakeTrk / lNRejShw<<endl;
00518 }
00519 cout << " Calibrated: " << ((fTrkMIPCalibOkay!=0)? " yes": "NO***") <<endl;
00520
00521 fCanvas ->Clear();
00522 fCanvas->Divide(2,1);
00523 fCanvas->cd(1);
00524 fUZAxis->GetXaxis()->SetRangeUser(minplane-5, maxplane+10);
00525 fUZAxis->Draw();
00526 for(unsigned int i=0; i<fUHits.size(); i++) fUHits[i].Draw();
00527 fCanvas->cd(2);
00528 fVZAxis->GetXaxis()->SetRangeUser(minplane-5, maxplane+10);
00529 fVZAxis->Draw();
00530 for(unsigned int i=0; i<fVHits.size(); i++) fVHits[i].Draw();
00531 fCanvas->Draw();
00532 fCanvas->Update();
00533
00534 fCanvas2->Clear();
00535 fCanvas2->Divide(2,1);
00536 fCanvas2->cd(1);
00537 fUZTrkQ->GetXaxis()->SetRangeUser(minplane-5, maxplane+10);
00538 fUZTrkQ->Draw("hist");
00539 fCanvas2->cd(2);
00540 fVZTrkQ->GetXaxis()->SetRangeUser(minplane-5, maxplane+10);
00541 fVZTrkQ->Draw("hist");
00542 fCanvas2->Draw();
00543 fCanvas2->Update();
00544
00545 char histname[1024];
00546 sprintf(histname, "mc_event_%d_%d.eps", header->GetRun(), header->GetSnarl());
00547 fCanvas->Print(histname);
00548 sprintf(histname, "mc_q_%d_%d.eps", header->GetRun(), header->GetSnarl());
00549 fCanvas2->Print(histname);
00550 }
00551 delete truth;
00552 delete beer;
00553 return JobCResult::kPassed;
00554 }
00555
00556
00557
00558 const Registry& MCAnalysis::DefaultConfig() const
00559 {
00563 static Registry r;
00564
00565
00566 std::string name = this->GetName();
00567 name += ".config.default";
00568 r.SetName(name.c_str());
00569
00570
00571 r.UnLockValues();
00572 r.Set("NameListIn","mergedigitlist");
00573 r.Set("Draw",0);
00574 r.LockValues();
00575
00576 return r;
00577 }
00578
00579
00580
00581 void MCAnalysis::Config(const Registry& r)
00582 {
00586 Int_t tmpi = 0;
00587 fNameListIn = r.GetCharString("NameListIn");
00588 if(r.Get("Draw", tmpi) ){
00589 fMakePictures = tmpi;
00590 }
00591 }
00592
00593
00594 void MCAnalysis::FillTrkInfo(const CandTrackHandle* track){
00595 for(int i=0; i<500; i++){
00596 fTrkQ[i] = 0;
00597 fTrkMIPDcosz[i] = 0;
00598 }
00599 fTrkStrips.clear();
00600
00601 const double nplanes = track->GetEndPlane() - track->GetVtxPlane();
00602 for(int ipln = track->GetVtxPlane(); ipln<=track->GetEndPlane(); ++ipln){
00603 double dcosz =track->GetVtxDirCosZ() + ((ipln - track->GetVtxPlane())/nplanes)*(track->GetEndDirCosZ() - track->GetVtxDirCosZ() );
00604 int ip1 =-1;
00605 int ip2 =-1;
00606 if(track->IsTPosValid(ipln-1) && track->IsTPosValid(ipln+1)){
00607 ip2 = ipln+1; ip1 = ipln-1;
00608 }else if(track->IsTPosValid(ipln) && track->IsTPosValid(ipln+1)){
00609 ip2 = ipln+1; ip1 = ipln;
00610 }else if(track->IsTPosValid(ipln) && track->IsTPosValid(ipln-1)){
00611 ip2 = ipln; ip1 = ipln-1;
00612 }
00613 if(ip1!=-1 && ip2!=-1){
00614 const double l = sqrt((track->GetU(ip2) - track->GetU(ip1))*(track->GetU(ip2) - track->GetU(ip1))
00615 + (track->GetV(ip2) - track->GetV(ip1))*(track->GetV(ip2) - track->GetV(ip1))
00616 + (track->GetZ(ip2) - track->GetZ(ip1))*(track->GetZ(ip2) - track->GetZ(ip1)) );
00617 dcosz = fabs(track->GetZ(ip2) - track->GetZ(ip1) ) / l;
00618 }
00619
00620 fTrkMIPDcosz[ipln] = track->GetPlaneCharge(ipln, CalStripType::kMIP)*dcosz;
00621 }
00622
00623
00624 TIter trkstpIter(track->GetDaughterIterator());
00625 while(const CandStripHandle* strip = dynamic_cast<const CandStripHandle*>(trkstpIter())){
00626 fTrkStrips.insert(strip->GetPlane()*192+strip->GetStrip());
00627 fTrkQ[strip->GetPlane()]+=strip->GetCharge(CalDigitType::kPE);
00628 if(strip->GetPlaneView()==PlaneView::kU){
00629 fUZTrkQ->Fill(strip->GetPlane(), track->GetStripCharge(strip, CalStripType::kMIP) );
00630 }else{
00631 fVZTrkQ->Fill(strip->GetPlane(), track->GetStripCharge(strip, CalStripType::kMIP) );
00632 }
00633 }
00634 }
00635
00636 void MCAnalysis::FillPlnInfo(const CandDeMuxDigitListHandle* digitlist)
00637 {
00638 for(int i=0; i<500; i++) fPlnQ[i] = 0;
00639 TIter digititer(digitlist->GetDaughterIterator());
00640 while (const CandDeMuxDigitHandle *digit = dynamic_cast<const CandDeMuxDigitHandle*>(digititer())) {
00641 if(!(digit->GetDeMuxDigitFlagWord()&CandDeMuxDigit::kXTalk)) {
00642 const int plane = digit->GetPlexSEIdAltL().GetPlane();
00643 if(plane>0 && plane<500) fPlnQ[plane]+=digit->GetCharge(CalDigitType::kPE);
00644
00645 }
00646 }
00647 }
00648
00649 int MCAnalysis::GetMaxTrkPlane(){
00650 int ntracklikeplanes=0;
00651 int maxtracklikeplane = 0;
00652 for(int ipln = 0; ipln<500 && ntracklikeplanes<6; ipln++){
00653 if(fTrkQ[ipln]>0){
00654 if(fTrkQ[ipln]/fPlnQ[ipln]>.8 && fTrkQ[ipln]>3.0){
00655 ntracklikeplanes++;
00656 maxtracklikeplane=ipln;
00657 }
00658 }
00659 }
00660 if(ntracklikeplanes!=6) maxtracklikeplane = 500;
00661 return maxtracklikeplane;
00662 }
00663
00664
00665
