//*-- AUTHOR : Jochen Markert //_HADES_CLASS_DESCRIPTION /////////////////////////////////////////////////////////////////////////////// // // HMdcDeDxHists // // Fills momentum vs mdcdedx (spline or kicktrack) for // full tracks (seg1+seg2), seg1, seg2 , system 0/1 and sector wise. // OPTIONS: // Momentum version:----------------------------------------------------------- // The histograms can be filled by HKicktrack or HSplineTrack: // The selection can be done with // setUseSpline(Bool_t) : kTRUE (default) - Fills from HSplineTrack // : kFALSE - Fills from HKickTrack // dEdx calculation:----------------------------------------------------------- // The mode of HMdcDeDx::calcDeDx(......., mode,module) can be switched // via setDeDxMode(Int_t m) : 0 (default) - Fills from HMdcSeg (only cells kept by the fitter) // 1 - Fills from HMdcClusFit (cells with weight 0 AND 1 are used) // via setDeDxModule(Int_t m): 0 - Fills for 1st module in HMdcSeg // 1 - Fills for 2st module in HMdcSeg // 2 (default) - Fills for both modules in HMdcSeg // // via setDeDxWindow(Float_t w) the window of the truncated mena can be set // in units of sigma (default 5.0); // via setDeDxMinimumWires(Int_t n) the minimal number of wires after the truncated // mean procedure can be set. // via setDeDxUseShifts(Bool_t shift) the normalization procedure of the dEdx calculation // cand be switched on or off (default=ON). // CPR mode:------------------------------------------------------------------- // Single wire distribution / ntuple filled for single/doubles of Jaro's CPR condition // Runs on special pepared input file (Stripped single/doubles with HKickTrack containing // 1 for singles and 2 for doubles ) // setCPR(Bool_t cpr) : kTRUE run in the special mode // setCPR(Bool_t cpr) : kFALSE (default) run in the normal mode // DEBUG mode:----------------------------------------------------------------- // If Debug mode is switched on (default OFF), // HMdcTrkCand.getFirstCandInd() =p*charge; // HMdcTrkCand.getNCandForSeg1() =single_double (only if cpr, 1=single,2=double); // The out put dedx of HMdcTrkCand will be refilled with the calculated values. // setDebug(Bool_t d) // /////////////////////////////////////////////////////////////////////////////// using namespace std; #include "hmdcdedxhists.h" #include "hmdcdedx.h" #include "hmdcseg.h" #include "hmdctrkcand.h" #include "hkicktrack.h" #include "hsplinetrack.h" #include "hmetamatch.h" #include "hmdcdef.h" #include "hmdctrackgdef.h" #include "hdebug.h" #include "hades.h" #include "hiterator.h" #include "hruntimedb.h" #include "hevent.h" #include "hcategory.h" #include "kickdef.h" #include "hmessagemgr.h" #include #include #include ClassImp(HMdcDeDxHists) HMdcDeDxHists::HMdcDeDxHists(void) { // Default constructor . initParameters(); } HMdcDeDxHists::HMdcDeDxHists(Text_t* name,Text_t* title) : HReconstructor(name,title) { // Constructor calls the constructor of class HReconstructor with the name // and the title as arguments. initParameters(); } HMdcDeDxHists::~HMdcDeDxHists(void) { // destructor deletes the iterator if (trkcanditer) delete trkcanditer; trkcanditer=0; if (kickiter) delete kickiter; kickiter=0; if (metamatchiter) delete metamatchiter; metamatchiter=0; } void HMdcDeDxHists::initParameters() { trkcandCat =0; segCat =0; trkcanditer=0; hasPrinted =kFALSE; useSpline =kTRUE; lockick.set(2,0,0); locseg.set(3,0,0,0); kickCat =0; kickiter =0; splineCat =0; metamatchCat =0; metamatchiter=0; online=kFALSE; // no online hists ntuple=kFALSE; // no ntuple debug =kFALSE; // no debugging cpr =kFALSE; // no single/doubles mode =0; // fill from segment module=2; // fill for both modules of a segment useshifts=kTRUE; // normalize dedx minwires=5; // minimum wires after trunce=ated mean window=5.0; // sigma window arround mean } void HMdcDeDxHists::setParContainers() { // creates the parameter containers } Bool_t HMdcDeDxHists::init(void) { // gets the pointers to HMdcTrkC and HMdcSeg categories // calls the function setParContainers() setParContainers(); trkcandCat=(HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcTrkCand)); if (!trkcandCat) { Error("HMdcDeDxHists:init()","HMdcTrkCand Category not available!"); exit(1); } segCat=(HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcSeg)); if (!segCat) { Error("HMdcDeDxHists:init()","HMdcSeg Category not available!"); exit(1); } trkcanditer =(HIterator *)((HCategory*)trkcandCat)->MakeIterator("native"); if(!useSpline) { kickCat=(HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catKickTrack)); if(kickCat) { kickiter=(HIterator *)((HCategory*)kickCat)->MakeIterator("native"); }else{ Warning("init()","No catKickTrack in input!"); } } if(useSpline) { metamatchCat=(HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catMetaMatch)); if(metamatchCat) { metamatchiter=(HIterator *)((HCategory*)metamatchCat)->MakeIterator("native"); }else{ Warning("init()","No catMetaMatch in input!"); } splineCat=(HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catSplineTrack)); if(!splineCat) { Warning("init()","No catSplineTrack in input!"); } } makeHists(); if(!hasPrinted)printStatus(); fActive=kTRUE; return kTRUE; } void HMdcDeDxHists::makeHists() { Char_t name [300]; Char_t title[300]; Int_t nbinX=200; Int_t nbinY=600; Float_t xlow=0.; Float_t xup=200.; Float_t ylow=-1500.; Float_t yup =1500.; output=new TFile("mdcdedx.root","RECREATE"); output->cd(); if(online) { result=new TCanvas("single_drift_cells","single drift cells",1000,800); result->Draw(); } if(ntuple) { celltuple=new TNtuple("dedx","dedx","type:p:s:phi:theta:seg:ncell:ncelltrunc:dedx:dedxtrunc:sigma:sigmatrunc:c0:c1:c2:c3:c4:c5:c6:c7:c8:c9:c10:c11:c12:c13:c14:c15:c16:c17:c18:c19:c20:c21:c22:c23:c24:c25:c26:c27:c28:c29"); } sprintf(name ,"hsinglecells"); sprintf(title,"dEdx of single driftcells"); hdriftcells[0]=new TH1F(name,title,nbinX,xlow,xup); hdriftcells[0]->SetXTitle("t2-t1 [ns]"); hdriftcells[0]->SetLineColor(2); sprintf(name ,"hdoublecells"); sprintf(title,"dEdx of single driftcells"); hdriftcells[1]=new TH1F(name,title,nbinX,xlow,xup); hdriftcells[1]->SetXTitle("t2-t1 [ns]"); hdriftcells[1]->SetLineColor(4); if(online) { result->cd(); hdriftcells[0]->Draw(); hdriftcells[1]->Draw("same"); } sprintf(name,"hdedx_p_all_trk"); sprintf(title,"dEdx_vs momentum all sectors track"); hdedx_p_all_trk=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_all_trk->SetXTitle("t2-t1 [ns]"); hdedx_p_all_trk->SetYTitle("momentum * charge [MeV]"); hdedx_p_all_trk->SetOption("COLZ"); for(Int_t sys=0;sys<2;sys++){ sprintf(name,"hdedx_p_all_trk_sys[%i]",sys); sprintf(title,"dEdx_vs momentum all sectors track sys %i",sys); hdedx_p_all_trk_sys[sys]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_all_trk_sys[sys]->SetXTitle("t2-t1 [ns]"); hdedx_p_all_trk_sys[sys]->SetYTitle("momentum * charge [MeV]"); hdedx_p_all_trk_sys[sys]->SetOption("COLZ"); for(Int_t s=0;s<6;s++){ sprintf(name,"hdedx_p_sec_trk_sys[%i][%i]",s,sys); sprintf(title,"dEdx_vs momentum sector %i track system %i",s,sys); hdedx_p_sec_trk_sys[s][sys]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_sec_trk_sys[s][sys]->SetXTitle("t2-t1 [ns]"); hdedx_p_sec_trk_sys[s][sys]->SetYTitle("momentum * charge [MeV]"); hdedx_p_sec_trk_sys[s][sys]->SetOption("COLZ"); } } for(Int_t seg=0;seg<2;seg++){ sprintf(name,"hdedx_p_all_seg[%i]",seg); sprintf(title,"dEdx_vs momentum all sectors seg %i",seg); hdedx_p_all_seg[seg]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_all_seg[seg]->SetXTitle("t2-t1 [ns]"); hdedx_p_all_seg[seg]->SetYTitle("momentum * charge [MeV]"); hdedx_p_all_seg[seg]->SetOption("COLZ"); for(Int_t sys=0;sys<2;sys++){ sprintf(name,"hdedx_p_all_seg_sys[%i][%i]",seg,sys); sprintf(title,"dEdx_vs momentum all sectors seg %i sys %i",seg,sys); hdedx_p_all_seg_sys[seg][sys]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_all_seg_sys[seg][sys]->SetXTitle("t2-t1 [ns]"); hdedx_p_all_seg_sys[seg][sys]->SetYTitle("momentum * charge [MeV]"); hdedx_p_all_seg_sys[seg][sys]->SetOption("COLZ"); } } for(Int_t s=0;s<6;s++){ sprintf(name,"hdedx_p_sec_trk[%i]",s); sprintf(title,"dEdx_vs momentum sector %i track",s); hdedx_p_sec_trk[s]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_sec_trk[s]->SetXTitle("t2-t1 [ns]"); hdedx_p_sec_trk[s]->SetYTitle("momentum * charge [MeV]"); hdedx_p_sec_trk[s]->SetOption("COLZ"); for(Int_t seg=0;seg<2;seg++){ sprintf(name,"hdedx_p_sec_seg[%i][%i]",s,seg); sprintf(title,"dEdx_vs momentum sector %i ioseg %i",s,seg); hdedx_p_sec_seg[s][seg]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_sec_seg[s][seg]->SetXTitle("t2-t1 [ns]"); hdedx_p_sec_seg[s][seg]->SetYTitle("momentum * charge [MeV]"); hdedx_p_sec_seg[s][seg]->SetOption("COLZ"); for(Int_t sys=0;sys<2;sys++){ sprintf(name,"hdedx_p_sec_seg_sys[%i][%i][%i]",s,seg,sys); sprintf(title,"dEdx_vs momentum sector %i ioseg %i system %i",s,seg,sys); hdedx_p_sec_seg_sys[s][seg][sys]=new TH2F(name,title,nbinX,xlow,xup,nbinY,ylow,yup); hdedx_p_sec_seg_sys[s][seg][sys]->SetXTitle("t2-t1 [ns]"); hdedx_p_sec_seg_sys[s][seg][sys]->SetYTitle("momentum * charge [MeV]"); hdedx_p_sec_seg_sys[s][seg][sys]->SetOption("COLZ"); } } } } void HMdcDeDxHists::printStatus() { // prints the parameters to the screen SEPERATOR_msg("-",60); INFO_msg(10,HMessageMgr::DET_MDC,"HMdcDeDxHistsSetup:"); /* gHades->getMsg()->info(10,DET_MDC,GetName() ,"ModeFlagCal1 = %i : 1 = NoStartandCal, 2 = StartandCal, 3 = NoStartandNoCal",ModeFlagCal1); */ SEPERATOR_msg("-",60); hasPrinted=kTRUE; } Int_t HMdcDeDxHists::execute(void) { // calculates dEdx for all Segments contained in HMdcTrkCand. // The result is filled to HMdctrk. HMdcDeDx* mdcdedx=(HMdcDeDx*)gHades->getRuntimeDb()->getContainer("MdcDeDx"); mdcdedx->setMinimumWires(minwires); mdcdedx->setWindow (window); mdcdedx->setUseShifts (useshifts); Float_t meanold1,sigmaold1,sigma1; UChar_t nwire1,nwiretrunc1; Float_t meanold2,sigmaold2,sigma2; UChar_t nwire2,nwiretrunc2; Float_t* array1; Float_t* array2; HMdcTrkCand* trkcand; Float_t dedx1=0; Float_t dedx2=0; Int_t indseg1,indseg2; Int_t nmatched; Float_t* array=0; Int_t size=-1; Float_t myvalues[42]; Int_t sizeval=42; Int_t firstval=12; Int_t sizetuple=30; HMdcSeg* seg1=0; HMdcSeg*seg2=0; Float_t phi1,theta1,phi2,theta2; if(!useSpline) { trkcanditer->Reset(); while ((trkcand=(HMdcTrkCand *)trkcanditer->Next())!=0) { indseg1=trkcand->getSeg1Ind(); indseg2=trkcand->getSeg2Ind(); array1=0; array2=0; dedx1 =dedx2 =-1; meanold1 =meanold2 =-1; sigmaold1 =sigmaold2 =-1; sigma1 =sigma2 =-1; nwiretrunc1=nwiretrunc2=0; nwire1 =nwire2 =0; seg1=(HMdcSeg*)segCat->getObject(indseg1); dedx1=mdcdedx->calcDeDx(seg1,&meanold1,&sigmaold1,&nwire1,&sigma1,&nwiretrunc1,mode,module); array=mdcdedx->getArray(size); if(size==-1) { Error("HMdcDeDxHists::execute()","Size of wire array not set properly!"); exit(1); } array1= new Float_t [size]; for(Int_t i=0;i-0.9) { seg2=(HMdcSeg*)segCat->getObject(indseg2); dedx2=mdcdedx->calcDeDx(seg2,&meanold2,&sigmaold2,&nwire2,&sigma2,&nwiretrunc2,mode,module); array=mdcdedx->getArray(size); array2= new Float_t [size]; for(Int_t i=0;igetNCandForSeg1(); if(nmatched>=0) { HKickTrack* kick; Int_t segindex; //------------------------looping over kicktrack to get segment kickiter->Reset(); while ((kick=(HKickTrack *)kickiter->Next())!=0) { locseg[0]=kick->getSector(); locseg[1]=0; locseg[2]=kick->getSegment(); Int_t s =kick->getSector(); Int_t sys =kick->getSystem(); //Int_t pid =kick->getPID(); //Int_t quality =kick->getQuality(); Float_t p =kick->getP(); Int_t charge =kick->getCharge(); Int_t single_double=kick->getRingId(); segindex=segCat->getIndex(locseg); if(indseg1==segindex) { phi1 =seg1->getPhi(); theta1=seg1->getTheta(); for(Int_t i=0;i0) { if(cpr)hdriftcells[single_double-1]->Fill(val); else hdriftcells[0] ->Fill(val); if(iFill(myvalues); } if(!cpr&&indseg2!=-1) { // only if no CPR investigation phi2 =seg2->getPhi(); theta2=seg2->getTheta(); // reset working array for ntuple for(Int_t i=0;i0) { hdriftcells[1]->Fill(val); if(iFill(myvalues); } } // draw online pictures if(online&&(gHades->getEventCounter()%100==0)) { result->Modified(); result->Update(); } // write p*charge and single_double flag // to HMdcTrkCand for easier investigation if(debug) { // not standard info will be over written!!!!!! trkcand->setdedxInnerSeg(dedx1); trkcand->setdedxOuterSeg(dedx2); trkcand->setSigmadedxInnerSeg(sigma1); trkcand->setSigmadedxOuterSeg(sigma2); trkcand->setNWirededxInnerSeg(nwire1); trkcand->setNWirededxOuterSeg(nwire2); trkcand->setNWireCutdedxInnerSeg(nwiretrunc1); trkcand->setNWireCutdedxOuterSeg(nwiretrunc2); trkcand->setNCandForSeg1((Short_t)single_double); trkcand->setFirstCandInd((Short_t)(p*charge)); } // filling hists hdedx_p_sec_seg_sys[s][0][sys]->Fill(dedx1,p*charge); hdedx_p_sec_seg [s][0] ->Fill(dedx1,p*charge);; hdedx_p_all_seg_sys[0][sys] ->Fill(dedx1,p*charge); hdedx_p_all_seg [0] ->Fill(dedx1,p*charge); if(indseg2>-0.9) { hdedx_p_sec_seg_sys[s][1][sys]->Fill(dedx2,p*charge); hdedx_p_sec_seg [s][1] ->Fill(dedx2,p*charge);; hdedx_p_sec_trk_sys[s][sys] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_sec_trk [s] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_all_seg_sys[1][sys] ->Fill(dedx2,p*charge); hdedx_p_all_seg [1] ->Fill(dedx2,p*charge); hdedx_p_all_trk_sys[sys] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_all_trk ->Fill((dedx1+dedx2)*0.5,p*charge); } } } //------------------------------------------------------------------- } if(array1)delete [] array1; if(array2)delete [] array2; } } if(useSpline) { HMetaMatch* metamatch; HSplineTrack* splinetrack; Int_t trkindex; Int_t splineindex; //------------------------looping over metamatch to get trkcand and spline track metamatchiter->Reset(); while ((metamatch=(HMetaMatch *)metamatchiter->Next())!=0) { trkindex =metamatch->getTrkCandInd(); trkcand =(HMdcTrkCand*)trkcandCat->getObject(trkindex); nmatched =trkcand->getNCandForSeg1(); if(nmatched==1) { // if only one candidate for inner seg splineindex=metamatch->getSplineInd(); if(splineindex!=-1) { // if spline track exists splinetrack=(HSplineTrack*)splineCat->getObject(splineindex); Int_t s =splinetrack->getSector(); Int_t sys =splinetrack->getSystem(); Float_t p =splinetrack->getP(); Int_t charge =splinetrack->getPolarity(); //Float_t pull =splinetrack->getPull(); if(sys!=-1) { // if meta hit exists indseg1=trkcand->getSeg1Ind(); indseg2=trkcand->getSeg2Ind(); array1=0; array2=0; dedx1 =dedx2 =-1; meanold1 =meanold2 =-1; sigmaold1 =sigmaold2 =-1; sigma1 =sigma2 =-1; nwiretrunc1=nwiretrunc2=0; nwire1 =nwire2 =0; seg1=(HMdcSeg*)segCat->getObject(indseg1); dedx1=mdcdedx->calcDeDx(seg1,&meanold1,&sigmaold1,&nwire1,&sigma1,&nwiretrunc1,mode,module); array=mdcdedx->getArray(size); if(size==-1) { Error("HMdcDeDxHists::execute()","Size of wire array not set properly!"); exit(1); } array1= new Float_t [size]; for(Int_t i=0;i-0.9) { seg2=(HMdcSeg*)segCat->getObject(indseg2); dedx2=mdcdedx->calcDeDx(seg2,&meanold2,&sigmaold2,&nwire2,&sigma2,&nwiretrunc2,mode,module); array=mdcdedx->getArray(size); array2= new Float_t [size]; for(Int_t i=0;igetPhi(); theta1=seg1->getTheta(); for(Int_t i=0;i0) { hdriftcells[0] ->Fill(val); if(iFill(myvalues); } if(indseg2!=-1) { phi2 =seg2->getPhi(); theta2=seg2->getTheta(); // reset working array for ntuple for(Int_t i=0;i0) { hdriftcells[1]->Fill(val); if(iFill(myvalues); } } // draw online pictures if(online&&(gHades->getEventCounter()%100==0)) { result->Modified(); result->Update(); } // write p*charge and single_double flag // to HMdcTrkCand for easier investigation if(debug) { // not standard info will be over written!!!!!! trkcand->setdedxInnerSeg(dedx1); trkcand->setdedxOuterSeg(dedx2); trkcand->setSigmadedxInnerSeg(sigma1); trkcand->setSigmadedxOuterSeg(sigma2); trkcand->setNWirededxInnerSeg(nwire1); trkcand->setNWirededxOuterSeg(nwire2); trkcand->setNWireCutdedxInnerSeg(nwiretrunc1); trkcand->setNWireCutdedxOuterSeg(nwiretrunc2); trkcand->setFirstCandInd((Short_t)(p*charge)); } hdedx_p_sec_seg_sys[s][0][sys]->Fill(dedx1,p*charge); hdedx_p_sec_seg [s][0] ->Fill(dedx1,p*charge);; hdedx_p_all_seg_sys[0][sys] ->Fill(dedx1,p*charge); hdedx_p_all_seg [0] ->Fill(dedx1,p*charge); if(indseg2>-0.9) { hdedx_p_sec_seg_sys[s][1][sys]->Fill(dedx2,p*charge); hdedx_p_sec_seg [s][1] ->Fill(dedx2,p*charge);; hdedx_p_sec_trk_sys[s][sys] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_sec_trk [s] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_all_seg_sys[1][sys] ->Fill(dedx2,p*charge); hdedx_p_all_seg [1] ->Fill(dedx2,p*charge); hdedx_p_all_trk_sys[sys] ->Fill((dedx1+dedx2)*0.5,p*charge); hdedx_p_all_trk ->Fill((dedx1+dedx2)*0.5,p*charge); } if(array1)delete [] array1; if(array2)delete [] array2; } } } } } return 0; } Bool_t HMdcDeDxHists::finalize() { output->cd(); if(ntuple)celltuple->Write(); hdriftcells[0]->Write(); hdriftcells[1]->Write(); hdedx_p_all_trk->Write(); for(Int_t sys=0;sys<2;sys++){ hdedx_p_all_trk_sys[sys]->Write(); for(Int_t s=0;s<6;s++){ hdedx_p_sec_trk_sys[s][sys]->Write(); } } for(Int_t seg=0;seg<2;seg++){ hdedx_p_all_seg[seg]->Write(); for(Int_t sys=0;sys<2;sys++){ hdedx_p_all_seg_sys[seg][sys]->Write(); } } for(Int_t s=0;s<6;s++){ hdedx_p_sec_trk[s]->Write(); for(Int_t seg=0;seg<2;seg++){ hdedx_p_sec_seg[s][seg]->Write(); for(Int_t sys=0;sys<2;sys++){ hdedx_p_sec_seg_sys[s][seg][sys]->Write(); } } } output->Save(); output->Close(); return kTRUE; }