#include "hwalleventplanef.h" #include "walldef.h" #include "hades.h" #include "hruntimedb.h" #include "hcategory.h" #include "hcategorymanager.h" #include "hrecevent.h" #include "hwallhit.h" #include "hwalleventplane.h" #include "hwalleventplanepar.h" #include "TRandom.h" #include "TFile.h" #include "TH1F.h" #include "TMath.h" #include "TVector2.h" using namespace std; //_HADES_CLASS_DESCRIPTION ///////////////////////////////////////////////////////////////////// // // HWallEventPlanef resonstructes the eventplane by HWallHits // ///////////////////////////////////////////////////////////////////// ClassImp(HWallEventPlaneF) Bool_t HWallEventPlaneF::fUseCorrection = kTRUE; HWallEventPlaneF::HWallEventPlaneF(const Text_t *name,const Text_t *title) : HReconstructor(name,title) { fWallHitCat = NULL; fWallEventPlaneCat = NULL; fWallEventPlanePar = NULL; } HWallEventPlaneF::~HWallEventPlaneF() { } Bool_t HWallEventPlaneF::init(void) { fWallHitCat = HCategoryManager::getCategory(catWallHit,kTRUE,"catWallHit"); if (!fWallHitCat) { Warning("init()","catWallHit not available!"); } fWallEventPlaneCat = HCategoryManager::addCategory(catWallEventPlane,"HWallEventPlane",10,"Wall",kFALSE); if(fUseCorrection) fWallEventPlanePar = (HWallEventPlanePar*)gHades->getRuntimeDb()->getContainer("WallEventPlanePar"); else fWallEventPlanePar = NULL; return kTRUE; } Int_t HWallEventPlaneF::execute(void) { if(!fWallHitCat) return 0; Float_t T1_cut = 0.0; Float_t T2_cut = 0.0; Float_t X_shift = 0.0; Float_t Y_shift = 0.0; //Float_t R0_cut = 0.0; Float_t Z1_cut_s = 0.0; Float_t Z1_cut_m = 0.0; Float_t Z1_cut_l = 0.0; //------------------------------------------------------------ // retrieve correction parameters if available if(fWallEventPlanePar) { T1_cut = fWallEventPlanePar->getT1Cut(); T2_cut = fWallEventPlanePar->getT2Cut(); X_shift = fWallEventPlanePar->getXShift(); Y_shift = fWallEventPlanePar->getYShift(); //R0_cut = fWallEventPlanePar->getR0Cut(); Z1_cut_s = fWallEventPlanePar->getZ1_cut_s(); Z1_cut_m = fWallEventPlanePar->getZ1_cut_m(); Z1_cut_l = fWallEventPlanePar->getZ1_cut_l(); } //------------------------------------------------------------ //------------------------------------------------------------ // loop over wall hits and select good cells for eventplane Int_t multWall = 0; Int_t cellNum = 0; Float_t cellTime = 0.0; Float_t cellCharge = 0.0; Float_t dEdxCut = 0.0; Float_t wallX,wallY,wallZ; Float_t wallXc, wallYc; //corrected reCentered TVector2 vect(0.,0.); TVector2 vsum(0.,0.); TVector2 eX(1.,0.); fcellsVect.reset(); // HWallHit* pWallHit = NULL; for(Int_t i = 0; i < fWallHitCat->getEntries(); i ++ ){ pWallHit = HCategoryManager::getObject(pWallHit,fWallHitCat,i); cellNum = pWallHit->getCell(); cellTime = pWallHit->getTime(); cellCharge = pWallHit->getCharge(); if(cellNum < 144 ){ dEdxCut = Z1_cut_s; } if(cellNum >= 144 && cellNum <= 207){ dEdxCut = Z1_cut_m; } if(cellNum > 207 ){ dEdxCut = Z1_cut_l; } if(cellTime > T1_cut && cellTime < T2_cut && cellCharge > dEdxCut){ // spectators selected multWall++; pWallHit->getXYZLab(wallX,wallY,wallZ); wallXc = wallX - X_shift; wallYc = wallY - Y_shift; vect.Set(wallXc, wallYc); vect = vect.Unit(); vsum += vect; fcellsVect.setCellVect(vect); } } //------------------------------------------------------------ //------------------------------------------------------------ // now we go over spectators and make calculations for A/B subevents Int_t NA = 0; Int_t NB = 0; Int_t choiceA = 1; // Int_t choiceB = 1; // Int_t multFWcells = fcellsVect.getNumbOfCells(); Float_t choice; TVector2 vectA(0.,0.), vectB(0.,0.); TVector2 vsumA(0.,0.), vsumB(0.,0.); for(Int_t ic = 0; ic < multFWcells; ic++){ choice = gRandom->Rndm(); // I can be proud that my algorithm has more flat distribution :) Float_t levelA = (multFWcells/2.-choiceA+1.)/Float_t(multFWcells); Float_t levelB = (multFWcells/2.-choiceB+1.)/Float_t(multFWcells); if(choice < levelA/(levelA+levelB)){ vectA = fcellsVect.getCellVector(ic); vsumA += vectA; choiceA++; NA++; } else { vectB = fcellsVect.getCellVector(ic); vsumB += vectB; choiceB++; NB++; } } //------------------------------------------------------------ //------------------------------------------------------------ // calculating eventplane angles Float_t phiEP = vsum.DeltaPhi(eX) *TMath::RadToDeg(); Float_t phiA = vsumA.DeltaPhi(eX) *TMath::RadToDeg(); Float_t phiB = vsumB.DeltaPhi(eX) *TMath::RadToDeg(); Float_t phiAB = vsumA.DeltaPhi(vsumB)*TMath::RadToDeg(); if(multFWcells == 0) phiEP = -1000.; if(NA == 0) phiA = -1000.; if(NB == 0) phiB = -1000.; if(NB == 0 || NA == 0) phiAB = -1000.; //------------------------------------------------------------ //------------------------------------------------------------ // fill output object if(fWallEventPlaneCat){ HWallEventPlane* eventplane = 0; Int_t index = 0; eventplane = HCategoryManager::newObject(eventplane,fWallEventPlaneCat,index); if(eventplane){ eventplane->setPhi ( phiEP); eventplane->setPhiA( phiA ); eventplane->setPhiB( phiB ); eventplane->setPhiAB(phiAB); eventplane->setNA(NA); eventplane->setNB(NB); } } //------------------------------------------------------------ return 0; }