// Based on alice_esd.C from ROOT tutorials // Modified for hyperk by Alex Finch #include "fitqun_def.h" #ifdef FITQUNEXISTS fitQunDisplay fiTQun; #endif #include "TColor.h" #include "HKEVInfo.h" /* Function definitions */ void make_gui(); void load_event(); void createGeometry(TGeoManager* geom,bool flatTube=kFALSE); void wcsim_load_event(int iTrigger,bool &firstTrackIsNeutrino,bool &secondTrackIsTarget ); void wcsim_load_cherenkov(int iTrigger); void wcsim_load_truth_tracks(int iTrigger,bool &firstTrackIsNeutrino,bool &secondTrackIsTarget); void fitqun_load_event(int entry); void UnrollView(double* pmtX ,double* pmtY,double* pmtZ,int location,float maxY,float maxZ,bool debug=false); TGeoVolume* createCylinder(TString name, float radius,float height); TEveViewer* New2dView(TString name,TGLViewer::ECameraType type, TEveScene* scene); void update_html_summary(int iTrigger,WCSimRootTrigger * wcsimrootTrigger,bool firstTrackIsNeutrino,bool secondTrackIsTarget); void selectTruthTrack(WCSimRootTrack * wcsimroottrack,bool& keep,bool& draw); TString ParticleName(int pdgCode); TString paletteText(); /* Global variables (!) */ float minTruthTrackEnergyToKeep=0.0; float minTruthTrackEnergyToDraw=100.0; int fPredefinedPalette{kRainBow}; WCSimRootGeom *wcsimrootgeom; WCSimRootEvent *wcsimrootEvent; WCSimRootTrigger *wcsimrootTrigger; TFile* WCSimFile, *fiTQunFile; TTree * wcsimGeoT; TTree * wcsimT; TTree* fiTQunTree; float maxX,maxY,maxZ,minZ; TGeoVolume *WorldVolume,*SimpleVolume; TEveElementList *FlatGeometry; TEveScene *UnrolledScene; TEveScene *flatGeometryScene; TEveViewer *UnrolledView; bool FITQUN,fAccumulateEvents{kTRUE},fDigitIsTime,fFirstTime{kTRUE}; bool fSimpleGeometry; TGTextButton *fGeometrySwitchButton; TGCheckButton *fAccumulateEventsBox; TGCheckButton *fColourIsTimeBox; TGNumberEntry *fTruthTrackMinimumEnergyWidget; TEveRGBAPalette *fHitPalette=0; TEveRGBAPaletteOverlay *fPalletteOverlay=0; TGLAnnotation *fPaletteTitle=0; TGLAnnotation *fPaletteTitle2=0; HKEVInfo *fHKEVInfo=nullptr; HKEVInfo *fHKEVInfo2D=nullptr; Int_t event_id= 0; // Current event id. TEveTrackList *gTrackList = 0; TEveGeoTopNode *geomRoot; #include "hyperk_esd_html_summary.h" HtmlSummary *fgHtmlSummary = 0; TGHtml *fgHtml = 0; #include "Picker.h" Picker* fPicker; Picker* fPicker2D; /******************************************************************************/ // Initialization and steering functions /******************************************************************************/ /* Entry point. Set everything up and then load the first event */ void hyperk_esd() { // Main function, initializes the application. /* Open files */ const char *filetypes[] = { "ROOT files", "*.root", "All files", "*", 0, 0 }; TString CurrentDirectory=gSystem->pwd(); TString originalDirectory=CurrentDirectory; TGFileInfo fi; fi.fFileTypes = filetypes; fi.fIniDir = StrDup(CurrentDirectory); cout<<" Please choose your WCSim file "<GetRoot(), 0, kFDOpen, &fi); if (!fi.fFilename) { cout<<" No WCSim file chosen "<cd(originalDirectory); fi.fIniDir = StrDup(CurrentDirectory); #ifdef FITQUNEXISTS FITQUN=kTRUE; #else FITQUN=kFALSE; #endif if(FITQUN) { cout<<" Please choose your fiTQun file, cancel if you don't have one "<GetRoot(), 0, kFDOpen, &fi); if (!fi.fFilename) { cout<<" No fiTQun file chosen "<Get("fiTQun"); fiTQun.Init(fiTQunTree); } } if(FITQUN) { fiTQunTree=(TTree *) fiTQunFile->Get("fiTQun"); fiTQun.Init(fiTQunTree); } gSystem->cd(originalDirectory); /* Initialise WCSim */ wcsimGeoT=(TTree *) WCSimFile->Get("wcsimGeoT"); wcsimrootgeom = new WCSimRootGeom (); wcsimGeoT -> SetBranchAddress ("wcsimrootgeom" ,&wcsimrootgeom ); wcsimGeoT->GetEntry(0) ; /* Initialise Eve */ gStyle->SetPalette(fPredefinedPalette ,0); TColor::InvertPalette(); // Default options: "FIV" - file-browser, command-line, GL-viewer. TEveManager::Create(kTRUE,"V"); gEve->GetBrowser()->SetTabTitle("3D View",TRootBrowser::kRight,0); gEve->GetDefaultViewer()->SetName("3D View"); TGLViewer* GLViewer =gEve->GetDefaultViewer()->GetGLViewer(); GLViewer->SetCurrentCamera(TGLViewer::kCameraPerspXOY); GLViewer->CurrentCamera().SetExternalCenter(kTRUE); GLViewer->CurrentCamera().SetCenterVec(0,0,0); GLViewer->UpdateScene(); // Create a clip box centree at 0, length 6000 TGLClipSet* thisClipSet=GLViewer->GetClipSet(); thisClipSet->SetClipType(TGLClip::kClipBox); TGLClipBox* thisClip = (TGLClipBox*) thisClipSet->GetCurrentClip(); thisClip->SetMode(TGLClip::kOutside); TGLVector3 min_point,max_point; const float maxDimension{4500.0}; min_point.Set(-maxDimension,-maxDimension, -maxDimension); max_point.Set(maxDimension,maxDimension,maxDimension); thisClip->Setup( min_point,max_point); GLViewer->SetClipAutoUpdate(kFALSE); GLViewer->RefreshPadEditor(); GLViewer->RequestDraw(); GLViewer->SetResetCamerasOnUpdate(kFALSE); // An overlay for extra information // /* Add information object */ fHKEVInfo = new HKEVInfo(GLViewer,0.0,0.99); fHKEVInfo->SetPermanent(" HyperK Event Display"); /* set up picking */ fPicker= new Picker(fHKEVInfo); GLViewer->Connect("Clicked(TObject*)", "Picker", fPicker, "Picked(TObject*)"); /* Initialise the html summary tab */ fgHtmlSummary = new HtmlSummary("HyperK Event Display Summary Table"); auto slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight()); fgHtml = new TGHtml(0, 100, 100); TEveWindowFrame *wf = slot->MakeFrame(fgHtml); fgHtml->MapSubwindows(); wf->SetElementName("Summary"); /* Initialise the geometry objects */ gSystem->Load("libGeom"); auto geom = new TGeoManager("HyperK", "HyperK Detector"); createGeometry(geom,kTRUE); TGeoNode* node = gGeoManager->GetTopNode(); geomRoot = new TEveGeoTopNode(gGeoManager, node); geomRoot->SetVisLevel(4); geomRoot->GetNode()->GetVolume()->SetVisibility(kFALSE); gEve->AddGlobalElement(geomRoot); gEve->Redraw3D(kTRUE); /* Set up the event control GUI */ gEve->GetBrowser()->GetTabRight()->SetTab(1); make_gui(); /* Initialise the unrolled geometry scene and view */ gEve->GetBrowser()->GetTabRight()->SetTab(0); UnrolledScene = gEve->SpawnNewScene("Unrolled Event"); flatGeometryScene = gEve->SpawnNewScene("Unrolled Geometry"); UnrolledView = New2dView("Unrolled View",TGLViewer::kCameraOrthoXnOY,UnrolledScene); flatGeometryScene->AddElement(FlatGeometry); UnrolledView->AddScene(flatGeometryScene); TEveSceneInfo* gSI= (TEveSceneInfo*) (UnrolledView->FindChild("SI - Geometry scene")); if(gSI!=NULL)gSI->Delete(); auto UGLV=UnrolledView->GetGLViewer(); fHKEVInfo2D = new HKEVInfo(UGLV,0.0,0.99); fHKEVInfo2D->SetPermanent(" HyperK Event Display"); fPicker2D= new Picker(fHKEVInfo2D); UGLV->Connect("Clicked(TObject*)", "Picker", fPicker2D, "Picked(TObject*)"); UGLV->AddOverlayElement(fPalletteOverlay); fPaletteTitle2 = new TGLAnnotation(UGLV,paletteText(), 0.7+(0.3/2.0)-0.05, 0.05+0.02+0.04); fHitPalette->SetInterpolate(kTRUE); fHitPalette->SetupColorArray(); fHitPalette->IncRefCount(); fHitPalette->IncRefCount(); /* Initialise FitQun */ if(FITQUN){ fiTQun.setLimits(maxX,maxY,maxZ,minZ); fiTQun.SetWCSimGeom(wcsimrootgeom); // fiTQun.CD.maxY=maxY; } /* Get the WCSim information */ wcsimT = (TTree*)WCSimFile ->Get("wcsimT"); wcsimrootEvent = new WCSimRootEvent (); wcsimT -> SetBranchAddress ("wcsimrootevent" ,&wcsimrootEvent); wcsimT->GetBranch("wcsimrootevent")->SetAutoDelete(kTRUE); wcsimT->GetEvent(0); cout<<" there are "<GetEntries()<<" events in this file."<GetDefaultGLViewer()->UpdateScene(); // Reset camera after the first event has been shown. gEve->Redraw3D(kTRUE); /* Load the first event */ load_event(); } /* Cartesian to Polar conversion */ void CartesianToPolar(double &mom,double &theta,double &phi,double p[3]){ mom=sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]); theta=acos(p[2]/mom); phi=atan2(p[1],p[0]); } /* Load an event. */ void load_event() { fHKEVInfo->ClearTemporary(); fHKEVInfo2D->ClearTemporary(); // Load event specified in global event_id. // The contents of previous event are removed. if(fFirstTime) { gStyle->SetPalette(fPredefinedPalette ,0); TColor::InvertPalette(); fFirstTime=kFALSE; } TEveEventManager* CurrentEvent =gEve->GetCurrentEvent(); if(! fAccumulateEvents ){ if( CurrentEvent != 0)CurrentEvent->DestroyElements(); if( UnrolledScene !=0 )UnrolledScene->DestroyElements(); } /* Get the next event from WCSim, loop over its triggers and put it into Eve */ wcsimT->GetEvent(event_id); fgHtmlSummary->Clear("D"); for(int iTrigger=0;iTriggerGetNumberOfEvents();iTrigger++) { wcsimrootTrigger = wcsimrootEvent->GetTrigger(iTrigger); bool firstTrackIsNeutrino,secondTrackIsTarget; wcsim_load_event(iTrigger,firstTrackIsNeutrino,secondTrackIsTarget); update_html_summary(iTrigger,wcsimrootTrigger,firstTrackIsNeutrino,secondTrackIsTarget); } if(FITQUN){ fitqun_load_event(event_id); fiTQun.CreateTable( fgHtmlSummary); } fgHtmlSummary->SetTitle(Form("HyperK Event Display Summary, for entry %i",event_id)); fgHtmlSummary->Build(); fgHtml->Clear(); fgHtml->ParseText((char*)fgHtmlSummary->Html().Data()); fgHtml->Layout(); TEveElement* top = gEve->GetCurrentEvent(); gEve->Redraw3D(kFALSE, kTRUE); } /******************************************************************************/ // GUI /******************************************************************************/ /* EvNavHandler class is needed to connect GUI signals to matching actions. */ class EvNavHandler { public: void Fwd() { if (event_id < wcsimT->GetEntries()) { if(UnrolledScene != 0){ TEveElement* CherenkovHits=UnrolledScene->FindChild(Form("CherenkovHits (Unrolled version) %i",event_id)); if(CherenkovHits!=NULL)UnrolledScene->RemoveElement(CherenkovHits); } ++event_id; load_event(); } else { printf("Already at last event.\n"); } } void Bck() { if (event_id > 0) { --event_id; load_event(); } else { printf("Already at first event.\n"); } } void AccumulateEventsBoxChanged(Bool_t s) { fAccumulateEvents=s; } void TruthTrackMinimumEnergyChanged() { minTruthTrackEnergyToDraw = fTruthTrackMinimumEnergyWidget->GetNumber(); cout<<" New minimum energy to draw truth tracks is "<GetNumberOfEvents();iTrigger++) { wcsimrootTrigger = wcsimrootEvent->GetTrigger(iTrigger); reload_wcsim_hits(iTrigger); } fPaletteTitle ->SetText(paletteText()); fPaletteTitle2->SetText(paletteText()); gEve->GetDefaultGLViewer()->UpdateScene(); UnrolledView->GetGLViewer()->UpdateScene(); } void reload_wcsim_hits(int iTrigger) { TEveEventManager* CurrentEvent =gEve->GetCurrentEvent(); if(CurrentEvent != 0){ TEveElement* CherenkovHits=CurrentEvent->FindChild(Form("Cherenkov Hits subevent %i",iTrigger)); if(CherenkovHits !=NULL)CurrentEvent->RemoveElement(CherenkovHits); } if(UnrolledScene != 0){ TEveElement* CherenkovHits=UnrolledScene->FindChild(Form("CherenkovHits (Unrolled version) %i",iTrigger)); if(CherenkovHits!=NULL)UnrolledScene->RemoveElement(CherenkovHits); } wcsim_load_cherenkov(iTrigger); } void SwitchGeometry() { if(fSimpleGeometry) { // simple to full fSimpleGeometry=kFALSE; gGeoManager->SetTopVolume(WorldVolume); TGeoNode* CurrentNode = gGeoManager->GetCurrentNode(); if(geomRoot!=NULL)gEve->GetGlobalScene()->RemoveElement(geomRoot); geomRoot = new TEveGeoTopNode(gGeoManager, CurrentNode); gEve->AddGlobalElement(geomRoot); fGeometrySwitchButton->SetText("Simple Geometry"); fGeometrySwitchButton->SetToolTipText("Switch to Simple Geometry"); } else { // full to simple fSimpleGeometry=kTRUE; gGeoManager->SetTopVolume(SimpleVolume); TGeoNode* CurrentNode = gGeoManager->GetCurrentNode(); if(geomRoot!=NULL)gEve->GetGlobalScene()->RemoveElement(geomRoot); geomRoot = new TEveGeoTopNode(gGeoManager, CurrentNode); gEve->AddGlobalElement(geomRoot); fGeometrySwitchButton->SetText("Full Geometry"); fGeometrySwitchButton->SetToolTipText("Switch to Full Geometry"); } gEve->GetDefaultGLViewer()->UpdateScene(); } void ClearOverlays() { fHKEVInfo->ClearTemporary(); fHKEVInfo2D->ClearTemporary(); } }; //______________________________________________________________________________ void make_gui() { // Create minimal GUI for event navigation, etc. TEveBrowser* browser = gEve->GetBrowser(); browser->StartEmbedding(TRootBrowser::kLeft); TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600); frmMain->SetWindowName("HyperK Event Display"); frmMain->SetCleanup(kDeepCleanup); TGCanvas* fCanvasWindow = new TGCanvas(frmMain, 10, 240); TGCompositeFrame* fFrame = new TGCompositeFrame(fCanvasWindow->GetViewPort(), 10, 100, kVerticalFrame); fFrame->SetLayoutManager(new TGVerticalLayout(fFrame)); fCanvasWindow->SetContainer(fFrame); // use hierarchical cleaning for container fFrame->SetCleanup(kDeepCleanup); EvNavHandler *fh = new EvNavHandler; TGLayoutHints* hfHints = new TGLayoutHints(kLHintsCenterX | kLHintsExpandX ,5,5,3,4); TGLayoutHints* itemHints = new TGLayoutHints(kLHintsCenterX | kLHintsCenterY ,5,5,3,4); TGGroupFrame* Group; { Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Event Navigation"); TGHorizontalFrame* hf = new TGHorizontalFrame(Group); { TString icondir( Form("%s/icons/", gSystem->Getenv("ROOTSYS")) ); TGPictureButton* b = 0; b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif")); b->SetToolTipText("Go back one event"); hf->AddFrame(b,itemHints); b->Connect("Clicked()", "EvNavHandler", fh, "Bck()"); b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoForward.gif")); b->SetToolTipText("Go forward one event"); hf->AddFrame(b,itemHints); b->Connect("Clicked()", "EvNavHandler", fh, "Fwd()"); } Group->AddFrame(hf); fCanvasWindow->AddFrame(Group,hfHints); /* Control to toggle from full to simple geometry */ Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Geometry Choice"); fSimpleGeometry=kTRUE; hf = new TGHorizontalFrame(Group); { if(!fSimpleGeometry) { fGeometrySwitchButton = new TGTextButton(Group, " Simple Geometry "); fGeometrySwitchButton->SetToolTipText("Switch to Simple Geometry"); } else { fGeometrySwitchButton = new TGTextButton(Group, " Full Geometry "); fGeometrySwitchButton->SetToolTipText("Switch to Full Geometry"); } fGeometrySwitchButton->Connect("Clicked()","EvNavHandler", fh,"SwitchGeometry()"); } Group->AddFrame( fGeometrySwitchButton); fCanvasWindow->AddFrame(Group,hfHints); /* Control to toggle event accumulation on or off */ fAccumulateEvents=kFALSE; Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Event accumulation"); hf = new TGHorizontalFrame(Group); { fAccumulateEventsBox = new TGCheckButton(hf, "Accumulate events.",1); fAccumulateEventsBox->SetState(kButtonUp); fAccumulateEventsBox->SetToolTipText("If this is checked, don't unload events from Eve."); fAccumulateEventsBox->Connect("Toggled(Bool_t)", "EvNavHandler", fh, "AccumulateEventsBoxChanged(Bool_t)"); hf->AddFrame(fAccumulateEventsBox); } Group->AddFrame(hf); fCanvasWindow->AddFrame(Group,hfHints); /* Control to toggle event colour is time */ fDigitIsTime=kTRUE; Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Colour for Digits"); hf = new TGHorizontalFrame(Group); { fColourIsTimeBox = new TGCheckButton(hf, "Digit Colour is Time.",1); fColourIsTimeBox->SetState(kButtonDown); fColourIsTimeBox->SetToolTipText("If this is checked, digit colour represents time, otherwise it is charge."); fColourIsTimeBox->Connect("Toggled(Bool_t)", "EvNavHandler", fh, "ColourIsTimeBoxChanged(Bool_t)"); hf->AddFrame(fColourIsTimeBox); } Group->AddFrame(hf); fCanvasWindow->AddFrame(Group,hfHints); /* Control to clear temporary information from the overlay */ Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Display"); hf = new TGHorizontalFrame(Group,500); { //TString icondir(gSystem->Getenv("$HKEVENTDISPLAY_ROOT/") ); TGPictureButton* b = 0; b = new TGPictureButton(hf, gClient->GetPicture("$HKEVENTDISPLAY_ROOT/eraser-icon-48.png")); b->SetToolTipText("Clear information about picked events from all event displays."); hf->AddFrame(b); b->Connect("Clicked()", "EvNavHandler", fh, "ClearOverlays()"); } Group->AddFrame(hf); fCanvasWindow->AddFrame(Group,hfHints); /* Control to set energy range for display of truth tracks */ Group = new TGGroupFrame(fCanvasWindow->GetContainer(),"Truth Track Display "); hf = new TGHorizontalFrame(Group); { //"Minimum Energy" (MeV)" TGLayoutHints* LH = new TGLayoutHints(kLHintsCenterX | kLHintsCenterY,5,5,3,4); TGLabel *label = new TGLabel(hf, "Minimum"); hf->AddFrame(label,LH); fTruthTrackMinimumEnergyWidget = new TGNumberEntry(hf,minTruthTrackEnergyToDraw, 9, 1, TGNumberFormat::kNESReal, //style TGNumberFormat::kNEAPositive, //input value filter TGNumberFormat::kNELNoLimits, //specify limits 0.,10000.); fTruthTrackMinimumEnergyWidget->Connect("ValueSet(Long_t)", "EvNavHandler", fh, "TruthTrackMinimumEnergyChanged()"); hf->AddFrame(fTruthTrackMinimumEnergyWidget,LH); label = new TGLabel(hf, "MeV"); hf->AddFrame(label,LH); // "(MeV" } Group->AddFrame(hf); fCanvasWindow->AddFrame(Group,hfHints); } frmMain->AddFrame(fCanvasWindow,new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 0, 0, 20, 20)); frmMain->MapSubwindows(); frmMain->Resize(); frmMain->MapWindow(); browser->StopEmbedding(); browser->SetTabTitle("Event Control", 0); /* Set up color palette object */ gStyle->SetPalette(fPredefinedPalette ,0); fHitPalette = new TEveRGBAPalette(0, 10000); fHitPalette->SetFixColorRange(kFALSE); fHitPalette->SetUnderflowAction( TEveRGBAPalette::kLA_Clip); fHitPalette->SetOverflowAction( TEveRGBAPalette::kLA_Clip); fHitPalette->SetupColorArray(); /* Create palette overlay once */ fPalletteOverlay =new TEveRGBAPaletteOverlay(fHitPalette, 0.7, 0.05, 0.3, 0.02); gEve->GetDefaultGLViewer()->AddOverlayElement(fPalletteOverlay); fPaletteTitle = new TGLAnnotation(gEve->GetDefaultGLViewer(),paletteText(), 0.7+(0.3/2.0)-0.05, 0.05+0.02+0.04); fPaletteTitle->SetState(TGLOverlayElement::kDisabled); fHitPalette->IncRefCount(); fHitPalette->IncRefCount(); } TString paletteText() { if(fDigitIsTime) return "PMT Hit Times"; else return "PMT Charge"; } /******************************************************************************/ // Event Loading /******************************************************************************/ //______________________________________________________________________________ void fitqun_load_event(int entry) { fiTQun.Process(entry); } void wcsim_load_event(int iTrigger,bool &firstTrackIsNeutrino,bool &secondTrackIsTarget ) { wcsim_load_cherenkov(iTrigger); wcsim_load_truth_tracks(iTrigger,firstTrackIsNeutrino,secondTrackIsTarget); } /* Call back function for displaying information about Cherenkov hits */ void PrintCherenkovHitInformation(Picker* picker,TEveDigitSet* ds, Int_t idx, TObject* obj) { if(obj) { WCSimRootCherenkovDigiHit *cDigiHit=(WCSimRootCherenkovDigiHit*) obj; if(cDigiHit) { int tubeId=cDigiHit->GetTubeId(); float Q=cDigiHit->GetQ(); float Time=cDigiHit->GetT(); WCSimRootPMT pmt = wcsimrootgeom -> GetPMT ( tubeId ); int location=pmt.GetCylLoc(); picker->Output(TString::Format("WCSim Cherenkov Hit. Id: %d \n Charge: %f, Time : %f\n Position : (%8.5G,%8.5G,%8.5G)",tubeId,cDigiHit->GetQ(), cDigiHit->GetT(),pmt.GetPosition (0),pmt.GetPosition (1),pmt.GetPosition (2))); } else cout<<" PrintCherenkovHitInformation called without a WCSimRootCherenkovDigiHit object."<ClearTemporary(); PrintCherenkovHitInformation(fPicker,ds, idx, obj); gEve->GetDefaultViewer()->GetGLViewer()->UpdateScene(); } void PrintCherenkovHitInformation2D(TEveDigitSet* ds, Int_t idx, TObject* obj) { fHKEVInfo2D->ClearTemporary(); PrintCherenkovHitInformation(fPicker2D, ds, idx, obj); UnrolledView->GetGLViewer()->UpdateScene(); } /* Loop over cherenkov digits and add them to Eve */ void wcsim_load_cherenkov(int iTrigger){ /* Set up box set for 3D view */ TEveBoxSet* CherenkovHits= new TEveBoxSet(Form("Cherenkov Hits subevent %i",iTrigger)); CherenkovHits->SetPalette(fHitPalette); CherenkovHits->SetEmitSignals(kTRUE); CherenkovHits->SetCallbackFoo(&PrintCherenkovHitInformation3D); CherenkovHits->Reset(TEveBoxSet::kBT_Cone, kFALSE, 64); /* set up box set for unrolled view */ TEveBoxSet* CherenkovHits2= new TEveBoxSet(Form("CherenkovHits (Unrolled version) %i",iTrigger)); CherenkovHits2->SetPalette(fHitPalette); CherenkovHits2->SetEmitSignals(kTRUE); CherenkovHits2->SetCallbackFoo(&PrintCherenkovHitInformation2D); CherenkovHits2->Reset(TEveBoxSet::kBT_Cone, kFALSE, 64); float PMTRadius = wcsimrootgeom->GetWCPMTRadius() ; int max=wcsimrootTrigger->GetNcherenkovdigihits(); float maxCharge=0; float minT=1e10; float maxT=-1e10; /* Loop over hits and add them to the digit set objects We need the mean and standard deviation of the data to set a sensible range for the colours used in the display/ */ float count=0; float mean=0; float M2=0; for (int i = 0; iGetCherenkovDigiHits()->At(i); int tubeId=cDigiHit->GetTubeId(); float Q=cDigiHit->GetQ(); float Time=cDigiHit->GetT(); maxCharge=TMath::Max(Q,maxCharge); minT = TMath::Min(minT,Time); maxT = TMath::Max(maxT,Time); WCSimRootPMT pmt = wcsimrootgeom -> GetPMT ( tubeId ); double pmtX = pmt.GetPosition (0); double pmtY = pmt.GetPosition (1); double pmtZ = pmt.GetPosition (2); float R=sqrt(pmtX*pmtX+pmtY*pmtY+pmtZ*pmtZ); int location=pmt.GetCylLoc(); double pmtX2=pmtX; double pmtY2=pmtY; double pmtZ2=pmtZ; if(R<100 || location <0 || location >2) cout<<"INVALID photo tube location . Tube id:" <GetNtrack(); if(ntrack>0) { int npar = wcsimrootTrigger->GetNpar(); int i; if(iTrigger==0) { firstTrackIsNeutrino=kFALSE; secondTrackIsTarget=kFALSE; // Take a first look and decide if we have a neutrino interaction. // If so the neutrino and target will need their start position // adjusting. TObject *element = (wcsimrootTrigger->GetTracks())->At(0); WCSimRootTrack *wcsimroottrack = dynamic_cast(element); int pdgCode=wcsimroottrack->GetIpnu(); if(abs(pdgCode)==12 || abs(pdgCode)==14) firstTrackIsNeutrino=kTRUE; if(ntrack>1) { TObject *element = (wcsimrootTrigger->GetTracks())->At(1); WCSimRootTrack *wcsimroottrack = dynamic_cast(element); int pdgCode=wcsimroottrack->GetIpnu(); if(abs(pdgCode)==2212 && wcsimroottrack->GetStart(0)==0.0 && wcsimroottrack->GetStart(1)==0.0 && wcsimroottrack->GetStart(2)==0.0 ) secondTrackIsTarget=kTRUE; } } // Loop through elements in the TClonesArray of WCSimTracks TEveElementList* TrueTracks = new TEveElementList(Form("Truth Tracks subevent %i",iTrigger)); for (i=0; iGetTracks())->At(i); WCSimRootTrack *wcsimroottrack = dynamic_cast(element); bool keep; bool draw; selectTruthTrack(wcsimroottrack,keep,draw); if(!keep)continue; int pdgCode=wcsimroottrack->GetIpnu(); if(pdgCode==0)continue; TString Name("Truth "+ParticleName(pdgCode)); if(iTrigger==0) { if(i==1 && secondTrackIsTarget) Name+=" (target) "; if(i==0 && firstTrackIsNeutrino) Name+=" (beam) "; } Name=Name+Form(" (%d) ",i); THKMCTrack* track=new THKMCTrack(Name); float Start[3]; float Stop[3]; float Mass=0; float Momentum=0; double Theta=0; double Phi=0; float Energy=0; track->SetElementTitle(Name); track->SetMainColor(kWhite); if(abs(pdgCode)==11)track->SetMainColor(kYellow); if(abs(pdgCode)==12)track->SetMainColor(kBlue); if(abs(pdgCode)==13)track->SetMainColor(kMagenta); if(abs(pdgCode)==14)track->SetMainColor(kGreen); if(abs(pdgCode)==2212)track->SetMainColor(kRed); Stop[0]=wcsimroottrack->GetStop(0); Stop[1]=wcsimroottrack->GetStop(1); Stop[2]=wcsimroottrack->GetStop(2); Start[0]=wcsimroottrack->GetStart(0); Start[1]=wcsimroottrack->GetStart(1); Start[2]=wcsimroottrack->GetStart(2); /* Implement fix for neutrino and target start positions */ if(i==0 && firstTrackIsNeutrino && iTrigger == 0) { track->SetNextPoint(wcsimroottrack->GetStop(0),wcsimroottrack->GetStop(1),maxZ*-1.1); Start[0]=wcsimroottrack->GetStop(0); Start[1]=wcsimroottrack->GetStop(1); Start[2]=maxZ*-1.1; } else { if(i==1 && secondTrackIsTarget && iTrigger == 0 ) { track->SetNextPoint(wcsimroottrack->GetStop(0),wcsimroottrack->GetStop(1),wcsimroottrack->GetStop(2)); Start[0]=wcsimroottrack->GetStop(0); Start[1]=wcsimroottrack->GetStop(1); Start[2]=wcsimroottrack->GetStop(2); } else { track->SetNextPoint(wcsimroottrack->GetStart(0),wcsimroottrack->GetStart(1),wcsimroottrack->GetStart(2)); } } if(abs(pdgCode)==12 || abs(pdgCode)==14 || abs(pdgCode)==16 ||abs(pdgCode)==22 ||abs(pdgCode)==111 ||abs(pdgCode)==310) track->SetLineStyle(2); track->SetNextPoint(wcsimroottrack->GetStop(0),wcsimroottrack->GetStop(1),wcsimroottrack->GetStop(2)); Mass=wcsimroottrack->GetM(); Momentum=wcsimroottrack->GetP(); Energy=wcsimroottrack->GetE(); double PT[3]; double PTtot=0; if(wcsimroottrack->GetP()>0){ PT[0]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(0); PT[1]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(1); PT[2]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(2); CartesianToPolar(PTtot,Theta,Phi,PT); } track->SetValues(Start,Stop,Mass,Momentum,Energy,Theta,Phi,draw); cout<<" truth track : momentum: "<AddElement(cone); //se3d->AddElement(track); //se2d->AddElement(cone2D); if(!draw) { track->SetRnrSelf(false); } TrueTracks->AddElement(track); } // End of loop over tracks if(ntrack>0)gEve->AddElement(TrueTracks); //Hists->Update(); } } /******************************************************************************/ // Geometry /******************************************************************************/ /// Create the geometry obects. /// There are Full , and Simple (3D) /// and unrolled (2D) void createGeometry(TGeoManager* geom, bool UseSimpleGeometry ) { /* Find the maximum values of X,Y,Z */ maxX=0; maxY=0; maxZ=0; minZ=0; // I assume z goes negative! float maxR=0; for(int tubeId=0;tubeIdGetWCNumPMT();tubeId++) { WCSimRootPMT pmt = wcsimrootgeom -> GetPMT ( tubeId ); int location=pmt.GetCylLoc(); if(location !=0 && location !=1 && location !=2)continue; double pmtX = pmt.GetPosition (0); double pmtY = pmt.GetPosition (1); double pmtZ = pmt.GetPosition (2); if(pmtX>maxX)maxX=pmtX; if(pmtY>maxY)maxY=pmtY; if(pmtZ>maxZ)maxZ=pmtZ; if(pmtZmaxR)maxR=R; } /// Now we know the maximum extent // MATERIALS, MIXTURES AND TRACKING MEDIA // Material: world TGeoMaterial* worldMaterial = new TGeoMaterial("world", 0.0,0.0,0.0); worldMaterial->SetIndex(0); worldMaterial->SetTransparency(60); // Medium: medium0 Int_t numed = 0; // medium number auto worldMedium = new TGeoMedium("worldMedium", numed,worldMaterial); //, par); float dx = 2*maxX; float dy = 2*maxY; float dz = 2*maxZ; TGeoShape *pworldbox_1 = new TGeoBBox("worldbox", dx,dy,dz); WorldVolume = new TGeoVolume("FullGeometry",pworldbox_1, worldMedium); WorldVolume->SetVisLeaves(kTRUE); SimpleVolume = new TGeoVolume("SimpleGeometry",pworldbox_1, worldMedium); SimpleVolume->SetVisLeaves(kTRUE); //UnrolledVolume = new TGeoVolume("UnrolledGeometry",pworldbox_1, pMed1); //UnrolledVolume->SetVisLeaves(kTRUE); //UnrolledVolume = new TGeoVolume("UnrolledGeometry",pworldbox_1, pMed1); //UnrolledVolume->SetVisLeaves(kTRUE); FlatGeometry = new TEveElementList("Flat Geometry"); TEveGeoShape *Shape; // SET TOP VOLUME OF GEOMETRY geom->SetMaxVisNodes(wcsimrootgeom->GetWCNumPMT()+1000); if(UseSimpleGeometry ) geom->SetTopVolume(SimpleVolume); else geom->SetTopVolume(WorldVolume); //FlatGeometry->SetTopVolume(UnrolledVolume); // SHAPES, VOLUMES AND GEOMETRICAL HIERARCHY // Shape: phototube type: TGeoSphere float PMTRadius = wcsimrootgeom->GetWCPMTRadius() ; Double_t rmin = 0.900000*PMTRadius; Double_t rmax = PMTRadius; Double_t theta1 = 0.000000; Double_t theta2 = 90.000000; Double_t phi1 = 0.000000; Double_t phi2 = 360.000000; /* Full geometry uses hemisphere for tube, UseSimpleGeometry=>just a disk */ Double_t zSize=0.1*PMTRadius; cout<<" create a phototube shape of dimenstion 0.0,"<SetVisLeaves(kTRUE); PhotoTubeVolume->SetLineColor(kYellow-5); auto PhotoSphereVolume = new TGeoVolume("photosphere",PhotoSphereShape, worldMedium); PhotoSphereVolume->SetVisLeaves(kTRUE); PhotoSphereVolume->SetLineColor(kYellow-5); TGeoRotation rota("rot",10,20,30); TGeoTranslation trans; TGeoCombiTrans *c1 = new TGeoCombiTrans(trans,rota); /* Read in location of ALL phototubes and add to the different scenes */ Double_t theta, phi; double pmtX2=0; double pmtY2=0; double pmtZ2=0; for(int tubeId=0;tubeIdGetWCNumPMT();tubeId++) { WCSimRootPMT pmt = wcsimrootgeom -> GetPMT ( tubeId ); double pmtX = pmt.GetPosition (0); double pmtY = pmt.GetPosition (1); double pmtZ = pmt.GetPosition (2); int location= pmt.GetCylLoc(); pmtX2=pmtX; pmtY2=pmtY; pmtZ2=pmtZ; // Work out x,y in unrolled view UnrollView(&pmtX2,&pmtY2,&pmtZ2,location,maxY,maxZ); theta=0.0; phi=0.0; float rad2deg=57.3; if(location==1) { float lengthXY=sqrt(pmtX*pmtX+pmtY*pmtY); float xd=pmtX/lengthXY; float yd=pmtY/lengthXY; TVector3 D(xd,yd,0.0); theta=D.Theta()*rad2deg; phi=D.Phi()*rad2deg; } else { theta=0.0; phi=0.0; if(location==2)theta=180.0; } if(location==0 || location ==1 || location ==2) { /* create a fake geometry object out of tevegeoshapes for the rolled out view*/ /* TGeoTranslation PhototubeUnrolledPositionMatrix("ExlodedShift",pmtX2,pmtY2,pmtZ2); auto shape = new TEveGeoShape(Form("Phototube %i",tubeId)); shape->SetShape(PhotoTubeShape); shape->SetTransMatrix(PhototubeUnrolledPositionMatrix); shape->SetMainColor(kYellow-5); shape->SetMainTransparency(70); if(tubeId<100) { cout<<" adding object to FlatGeometry at "<AddElement(shape); } */ /* now the 'normal' root geometry objects */ TGeoRotation TubeRotation("rotation",phi-90.0,theta,0.0); //D.Phi(),D.Theta,0.0); TGeoTranslation PhototubePositionMatrix("shift",pmtX,pmtY,pmtZ); TGeoCombiTrans *ShiftAndTwist = new TGeoCombiTrans(PhototubePositionMatrix,TubeRotation); //SimpleVolume->AddNode(PhotoTubeVolume, tubeId, ShiftAndTwist); WorldVolume-> AddNode(PhotoSphereVolume, tubeId, ShiftAndTwist); } } TGeoVolume* OutlineVolume = createCylinder("Inner Detector",maxR,maxZ); OutlineVolume->SetVisLeaves(kTRUE); OutlineVolume->SetLineColor(kGray); WorldVolume->AddNode(OutlineVolume,0); SimpleVolume->AddNode(OutlineVolume,0); geom->CloseGeometry(); /* Create a simple background box object to represent the unrolled cylinder in the unrolled view */ TEveBox* cylinder=new TEveBox("UnrolledBarell"); pmtX2=0.0; pmtY2=maxY; pmtZ2=maxZ; // Work out x,y in unrolled view UnrollView(&pmtX2,&pmtY2,&pmtZ2,1,maxY,maxZ); pmtZ2=10000; cylinder->SetVertex(0,-pmtX2,-pmtY2,pmtZ2-100); cylinder->SetVertex(1,-pmtX2,+pmtY2,pmtZ2-100); cylinder->SetVertex(2, pmtX2,+pmtY2,pmtZ2-100); cylinder->SetVertex(3, pmtX2,-pmtY2,pmtZ2-100); cylinder->SetVertex(4,-pmtX2,-pmtY2,pmtZ2+100); cylinder->SetVertex(5,-pmtX2,+pmtY2,pmtZ2+100); cylinder->SetVertex(6, pmtX2,+pmtY2,pmtZ2+100); cylinder->SetVertex(7, pmtX2,-pmtY2,pmtZ2+100); cylinder->SetMainColor(kGray); cylinder->SetMainTransparency(80); //cylinder->SetLineColor(kGray); FlatGeometry->AddElement(cylinder); for(int location=0;location<3;location+=2) { auto arr=new TEveArrow(); arr->SetTubeR(maxR); arr->SetConeR(0.1); arr->SetConeL(0.0); arr->SetDrawQuality(100); arr->SetMainColor(kGray); //arr->SetLineColor(kGray); arr->SetMainTransparency(80); pmtX2=0.0; pmtY2=0.0; pmtZ2=maxZ; UnrollView(&pmtX2,&pmtY2,&pmtZ2,location,maxY,maxZ); arr->SetOrigin(pmtX2,pmtY2,10000); FlatGeometry->AddElement(arr); } } // Convert x,y,z of phototubes to position in Unrolled view void UnrollView(double* pmtX ,double* pmtY,double* pmtZ,int location,float maxY,float maxZ,bool debug) { if(location==0) { // cout<<" add 2* "<xshift)x=x-(xshift*2); *pmtX=x; //cout<<" x,y,z of Unrolled phototube "<<*pmtX<<" "<<*pmtY<<" "<<*pmtZ<SetIndex(0); pMaterial->SetTransparency(86); TGeoMedium* pMedium = new TGeoMedium(name, 1,pMaterial); TGeoXtru *xtru = new TGeoXtru(2); int nvertices=1000; double x[nvertices]; double y[nvertices]; #define PI 3.14159 float stepSize=2*PI/(float)nvertices; float angle=0.0; for(int step=0;stepDefinePolygon(nvertices,x,y); xtru->DefineSection(0,-height); xtru->DefineSection(1,height); return new TGeoVolume(name,xtru,pMedium); } TEveViewer* New2dView(TString name,TGLViewer::ECameraType type, TEveScene* scene) { TEveViewer* View =gEve->SpawnNewViewer(name,name); View->AddScene(scene); // add the special scene that only applies to this view View->AddScene(gEve->GetGlobalScene()); // add the geometry information View->GetGLViewer()->SetCurrentCamera(type); View->GetGLViewer()->GetLightSet()->SetLight( TGLLightSet::kLightFront,kFALSE); View->GetGLViewer()->SetResetCamerasOnUpdate(kFALSE); return View; } //______________________________________________________________________________ void update_html_summary(int iTrigger,WCSimRootTrigger * wcsimrootTrigger,bool firstTrackIsNeutrino,bool secondTrackIsTarget) { // Update summary of current event. int ntrack = wcsimrootTrigger->GetNtrack(); if(ntrack==0)return; int nused=0; HtmlObjTable *table; if(iTrigger==0) table = fgHtmlSummary->AddTable(Form("Truth Tracks subevent: %i",iTrigger), 13,kTRUE,"first"); else table = fgHtmlSummary->AddTable(Form("Truth Tracks subevent: %i",iTrigger), 13,kTRUE,"other"); table->SetLabel(0, "Type"); table->SetLabel(1, "Start X"); table->SetLabel(2, "Start Y"); table->SetLabel(3, "Start Z"); table->SetLabel(4, "Stop X"); table->SetLabel(5, "Stop Y"); table->SetLabel(6, "Stop Z"); table->SetLabel(7, "Mass "); table->SetLabel(8, "Momentum "); table->SetLabel(9, "Energy "); table->SetLabel(10, "Theta"); table->SetLabel(11, "Phi"); table->SetLabel(12, "Drawn"); int i; int used=0; // Loop through elements in the TClonesArray of WCSimTrack for (i=0; iGetTracks())->At(i); WCSimRootTrack *wcsimroottrack = dynamic_cast(element); int pdgCode=wcsimroottrack->GetIpnu(); if(pdgCode==0)continue; TString Name(Form(" #%d (%d) ",i,pdgCode)); if(pdgCode==11)Name+="electron"; if(pdgCode==-11)Name+="positron"; if(pdgCode==12)Name+="electron neutrino"; if(pdgCode==-12)Name+="electron anti-neutrino"; if(pdgCode==13)Name+="muon"; if(pdgCode==-13)Name+="anti-muon"; if(pdgCode==14)Name+="muon neutrino"; if(pdgCode==-14)Name+="muon anti-neutrino"; if(pdgCode==2212)Name+="proton"; if(pdgCode==-2212)Name+="anti-proton"; if(pdgCode==22)Name+="photon"; //Name=Name+Form(" (%d) ",i); if(i==0 && firstTrackIsNeutrino) Name+=" (beam) "; if(i==1 && secondTrackIsTarget) Name+=" (target) "; table->SetRowName(used, Name); for(int l=0;l<3;l++){ table->SetValue(l+1,used,wcsimroottrack->GetStart(l)); table->SetValue(l+4,used,wcsimroottrack->GetStop(l)); } if(i==0 && firstTrackIsNeutrino && iTrigger==0) { table->SetValue(1,used,wcsimroottrack->GetStop(0)); table->SetValue(2,used,wcsimroottrack->GetStop(1)); table->SetValue(3,used,-1.1*maxZ); } if(i==1 && secondTrackIsTarget && iTrigger==0) { table->SetValue(1,used,wcsimroottrack->GetStop(0)); table->SetValue(2,used,wcsimroottrack->GetStop(1)); table->SetValue(3,used,wcsimroottrack->GetStop(2)); } // Int_t GetIpnu() const { return fIpnu;} table->SetValue(7,used,wcsimroottrack->GetM()); table->SetValue(8,used,wcsimroottrack->GetP()); table->SetValue(9,used,wcsimroottrack->GetE()); double PT[3]; double PTtot=0; double Theta=0; double Phi=0; if(wcsimroottrack->GetP()>0){ PT[0]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(0); PT[1]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(1); PT[2]=wcsimroottrack->GetP()*wcsimroottrack->GetPdir(2); CartesianToPolar(PTtot,Theta,Phi,PT); } table->SetValue(10,used,Theta); table->SetValue(11,used,Phi); bool keep; bool draw; selectTruthTrack(wcsimroottrack,keep,draw); if(draw) table->SetStringValue(12,used,"drawn"); else table->SetStringValue(12,used,"not drawn"); //for(int l=0;l<3;l++) // table->SetValue(10+l,used,wcsimroottrack->GetPdir(l)); used++; } // End of loop over tracks // } void selectTruthTrack(WCSimRootTrack * wcsimroottrack,bool& keep,bool& draw) { float Energy=wcsimroottrack->GetE(); keep=true; draw=true; if(Energy