# include # include # include # include # include # include # include # include "hmdcmetaaligner.h" # include "hmdcaligner.h" # include "hades.h" # include "hdetector.h" # include "hmdcdetector.h" # include "hmdcmodulegeometry.h" # include "hmdcgeomstruct.h" # include "hevent.h" # include "hmatrixcategory.h" # include "hmatrixcatiter.h" # include "hmdchit.h" # include "hmdchitaux.h" # include "hruntimedb.h" # include "hspectrometer.h" # include "mdcsdef.h" # include "../tof/tofdef.h" # include "../shower/showerdef.h" ClassImp(HMdcMetaAligner) //*-- AUTHOR : Hector Alvarez-Pol //*-- Date: 08/2001 //*-- Last Update: 25/08/01 //*-- Copyright: GENP (Univ. Santiago de Compostela) //////////////////////////////////////////////////////////////////////// // // HMdcMetaAligner (ongoing work) // // Performs the TOF/SHOWER software alignment (first part). // // First, a set of Hits from MDCs is used to follow the tracks // up to the SHOWER/TOF. The Hit selection repeates the scheme // from hmdcalignernew.cc. Then the differences in X, Y // gives information about the position // //////////////////////////////////////////////////////////////////////// Int_t HMdcMetaAligner::fRecCount=0; TFile* HMdcMetaAligner::fOutRoot=0; Int_t MA_DEBUG=0; HMdcMetaAligner::HMdcMetaAligner(void) : HReconstructor() { // // Default constructor. // fLoc.setNIndex(5); fHits = new HMdcMetaHit(4); fLoc.set(5,0,0,1,2,3); //dummy sector 0 fNumMods = 4; initDefaults(); } HMdcMetaAligner::HMdcMetaAligner(Int_t sector, Int_t modA, Int_t modB) : HReconstructor() { // // Constructor including module election. Alignment procedure // proyects hits of modA and modB on META coordinate system and compares // fLoc.setNIndex(3); fHits = new HMdcMetaHit(2); fLoc.set(3,sector,modA,modB); fNumMods = 2; initDefaults(); } HMdcMetaAligner::HMdcMetaAligner(Text_t* name, Text_t* title, Int_t sector, Int_t modA, Int_t modB, Int_t modC=-1, Int_t modD=-1) : HReconstructor(name, title) { // // Constructor including module election (and name and title, which // seems to be very important). Alignment procedure // proyects hits of modA and modB on META coordinate system and compares // if(modC == -1) { fHits = new HMdcMetaHit(2); fLoc.setNIndex(3); fLoc.set(3,sector,modA,modB); fNumMods = 2; } else if(modD == -1){ fHits = new HMdcMetaHit(3); fLoc.setNIndex(4); fLoc.set(4,sector,modA,modB,modC); fNumMods = 3; } else { fHits = new HMdcMetaHit(4); fLoc.setNIndex(5); fLoc.set(5,sector,modA,modB,modC,modD); fNumMods = 4; } initDefaults(); } HMdcMetaAligner::HMdcMetaAligner(Text_t* name, Text_t* title, Int_t sector, Int_t modA) : HReconstructor(name, title) { // // Constructor for only one MDC // fHits=new HMdcMetaHit(1); fLoc.setNIndex(3); fLoc.set(3,sector,modA,0); //the last one is a dummy fNumMods = 1; initDefaults(); } void HMdcMetaAligner::initDefaults(void) { // // Inits common defaults // if(fNumMods>1){ fRotMat = new HGeomRotation[fNumMods-1]; fTranslation = new HGeomVector[fNumMods-1]; fEuler = new HGeomVector[fNumMods-1]; } fMRotMat = new HGeomRotation[9]; fMTranslation = new HGeomVector[9]; fMEuler = new HGeomVector[9]; if(fNumMods>1) fDiscart = new Int_t[fNumMods-1]; fError = new Double_t[6]; fSDiscart = 0; fTDiscart = new Int_t[8]; if(fNumMods>1){ fHitsMdc = new Int_t[fNumMods]; fHitsFoundInWindow = new Int_t[fNumMods-1]; fHitsFoundAndUnique = new Int_t[fNumMods-1]; } else fHitsMdc = new Int_t[1]; fHitsFoundInTofWindow = new Int_t[8]; fHitsFoundInTofAndUnique = new Int_t[8]; if(fNumMods>1){ for(Int_t i=0;i 90% of gauss distribution fYSigmas = YSigmas; //1.96 -> 95% of gauss distribution fS0Sigmas = S0Sigmas; //2.58 -> 99% of gauss distribution fS1Sigmas = S1Sigmas; //3.29 -> 99.9% of gauss distribution } void HMdcMetaAligner::setHistograms(void) { // // Inits histograms // Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; fRecCount++; Char_t title[50], tmp[50]; if(!fOutRoot) fOutRoot = new TFile(fNameRoot,"UPDATE"); //Histograms for studying residuals in projections if(fNumMods>3) { CvsDinDCS_X = new TH1F*[fHistNum]; CvsDinDCS_Y = new TH1F*[fHistNum]; CvsDinDCS_XSlope = new TH1F*[fHistNum]; CvsDinDCS_YSlope = new TH1F*[fHistNum]; BvsDinDCS_X = new TH1F*[fHistNum]; BvsDinDCS_Y = new TH1F*[fHistNum]; BvsDinDCS_XSlope = new TH1F*[fHistNum]; BvsDinDCS_YSlope = new TH1F*[fHistNum]; AvsDinDCS_X = new TH1F*[fHistNum]; AvsDinDCS_Y = new TH1F*[fHistNum]; AvsDinDCS_XSlope = new TH1F*[fHistNum]; AvsDinDCS_YSlope = new TH1F*[fHistNum]; DvsCinCCS_X = new TH1F*[fHistNum]; DvsCinCCS_Y = new TH1F*[fHistNum]; DvsCinCCS_XSlope = new TH1F*[fHistNum]; DvsCinCCS_YSlope = new TH1F*[fHistNum]; DvsAinACS_X = new TH1F*[fHistNum]; DvsAinACS_Y = new TH1F*[fHistNum]; DvsAinACS_XSlope = new TH1F*[fHistNum]; DvsAinACS_YSlope = new TH1F*[fHistNum]; DvsBinBCS_X = new TH1F*[fHistNum]; DvsBinBCS_Y = new TH1F*[fHistNum]; DvsBinBCS_XSlope = new TH1F*[fHistNum]; DvsBinBCS_YSlope = new TH1F*[fHistNum]; SqrDistToD = new TH1F*[fHistNum]; } if(fNumMods>2) { BvsCinCCS_X = new TH1F*[fHistNum]; BvsCinCCS_Y = new TH1F*[fHistNum]; BvsCinCCS_XSlope = new TH1F*[fHistNum]; BvsCinCCS_YSlope = new TH1F*[fHistNum]; AvsCinCCS_X = new TH1F*[fHistNum]; AvsCinCCS_Y = new TH1F*[fHistNum]; AvsCinCCS_XSlope = new TH1F*[fHistNum]; AvsCinCCS_YSlope = new TH1F*[fHistNum]; CvsBinBCS_X = new TH1F*[fHistNum]; CvsBinBCS_Y = new TH1F*[fHistNum]; CvsBinBCS_XSlope = new TH1F*[fHistNum]; CvsBinBCS_YSlope = new TH1F*[fHistNum]; CvsAinACS_X = new TH1F*[fHistNum]; CvsAinACS_Y = new TH1F*[fHistNum]; CvsAinACS_XSlope = new TH1F*[fHistNum]; CvsAinACS_YSlope = new TH1F*[fHistNum]; XChi2Hist = new TH1F*[fHistNum]; YChi2Hist = new TH1F*[fHistNum]; TotalChi2 = new TH1F*[fHistNum]; SqrDistToA = new TH1F*[fHistNum]; SqrDistToB = new TH1F*[fHistNum]; SqrDistToC = new TH1F*[fHistNum]; SqrDist = new TH1F*[fHistNum]; } if(fNumMods>1) { AvsBinBCS_X = new TH1F*[fHistNum]; AvsBinBCS_Y = new TH1F*[fHistNum]; AvsBinBCS_XSlope = new TH1F*[fHistNum]; AvsBinBCS_YSlope = new TH1F*[fHistNum]; BvsAinACS_X = new TH1F*[fHistNum]; BvsAinACS_Y = new TH1F*[fHistNum]; BvsAinACS_XSlope = new TH1F*[fHistNum]; BvsAinACS_YSlope = new TH1F*[fHistNum]; } fResShowerX = new TH1F*[fMetaHistNum]; fResShowerY = new TH1F*[fMetaHistNum]; fResTof0X = new TH1F*[fMetaHistNum]; fResTof0Y = new TH1F*[fMetaHistNum]; fResTof1X = new TH1F*[fMetaHistNum]; fResTof1Y = new TH1F*[fMetaHistNum]; fResTof2X = new TH1F*[fMetaHistNum]; fResTof2Y = new TH1F*[fMetaHistNum]; fResTof3X = new TH1F*[fMetaHistNum]; fResTof3Y = new TH1F*[fMetaHistNum]; fResTof4X = new TH1F*[fMetaHistNum]; fResTof4Y = new TH1F*[fMetaHistNum]; fResTof5X = new TH1F*[fMetaHistNum]; fResTof5Y = new TH1F*[fMetaHistNum]; fResTof6X = new TH1F*[fMetaHistNum]; fResTof6Y = new TH1F*[fMetaHistNum]; fResTof7X = new TH1F*[fMetaHistNum]; fResTof7Y = new TH1F*[fMetaHistNum]; if (fNumMods>3){ Int_t modC = fLoc[3]; Int_t modD = fLoc[4]; sprintf(title,"%s%i%s%i%s%i%s%i%s%i","Sector_",sector,"_ModA_",modA, "_ModB_",modB,"_ModC_",modC,"_ModD_",modD); sprintf(tmp,"%s%s%i%s%i%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC,"_ModD_",modD); fAlignMeta = new TTree(tmp,title); fAlignMeta->Branch("hits","HMdcMetaHit",&fHits,64000); sprintf(tmp,"%s%s%i%s%i%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC,"_ModD_",modD); fAlignMetaCut = new TTree(tmp,title); fAlignMetaCut->Branch("hits","HMdcMetaHit",&fHits,64000); static Char_t newDirName[255]; sprintf(newDirName,"%s%s%i%s%i%s%i%s%i%s%i","MetaAligner_","S_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC,"_ModD_",modD); fOutRoot->mkdir(newDirName,newDirName); fOutRoot->cd(newDirName); //binning Int_t bin=200, binS=200, binChi=200, binDist=200; Int_t min=-100,max=100,minS=-1,maxS=1; Int_t minChi=0, maxChi=10, minDist=0,maxDist=100; for(Int_t index=0;index2){ Int_t modC = fLoc[3]; sprintf(title,"%s%i%s%i%s%i%s%i","Sector_",sector,"_ModA_",modA, "_ModB_",modB,"_ModC_",modC); sprintf(tmp,"%s%s%i%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC); fAlignMeta = new TTree(tmp,title); fAlignMeta->Branch("hits","HMdcMetaHit",&fHits,64000); sprintf(tmp,"%s%s%i%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC); fAlignMetaCut = new TTree(tmp,title); fAlignMetaCut->Branch("hits","HMdcMetaHit",&fHits,64000); static Char_t newDirName[255]; sprintf(newDirName,"%s%s%i%s%i%s%i%s%i","MetaAligner_","S_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC); fOutRoot->mkdir(newDirName,newDirName); fOutRoot->cd(newDirName); //binning Int_t bin=200, binS=200, binChi=200, binDist=200; Int_t min=-100,max=100,minS=-1,maxS=1; Int_t minChi=0, maxChi=10, minDist=0,maxDist=100; for(Int_t index=0;index1){ sprintf(title,"%s%i%s%i%s%i","Sector_",sector, "_ModA_",modA,"_ModB_",modB); sprintf(tmp,"%s%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fAlignMeta = new TTree(tmp,title); fAlignMeta->Branch("hits","HMdcMetaHit",&fHits,64000); sprintf(tmp,"%s%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fAlignMetaCut = new TTree(tmp,title); fAlignMetaCut->Branch("hits","HMdcMetaHit",&fHits,64000); static Char_t newDirName[255]; sprintf(newDirName,"%s%s%i%s%i%s%i","MetaAligner_","S_",sector, "_ModA_",modA,"_ModB_",modB); fOutRoot->mkdir(newDirName,newDirName); fOutRoot->cd(newDirName); //binning Int_t bin=200; Int_t binS=200; Int_t min=-100,max=100,minS=-1,maxS=1; for(Int_t index=0;indexBranch("hits","HMdcMetaHit",&fHits,64000); sprintf(tmp,"%s%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA); fAlignMetaCut = new TTree(tmp,title); fAlignMetaCut->Branch("hits","HMdcMetaHit",&fHits,64000); static Char_t newDirName[255]; sprintf(newDirName,"%s%s%i%s%i","MetaAligner_","S_",sector, "_ModA_",modA); fOutRoot->mkdir(newDirName,newDirName); fOutRoot->cd(newDirName); } //Histogramas comunes if(fNumMods>1){ sprintf(tmp,"%s%s%i%s%i%s%i","fResX","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fResX = new TH1F(tmp,title,2000,-1000,1000); sprintf(tmp,"%s%s%i%s%i%s%i","fResY","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fResY = new TH1F(tmp,title,2000,-1000,1000); sprintf(tmp,"%s%s%i%s%i%s%i","fResS0","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fResS0 = new TH1F(tmp,title,2000,-1,1); sprintf(tmp,"%s%s%i%s%i%s%i","fResS1","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); fResS1 = new TH1F(tmp,title,2000,-1,1); sprintf(tmp,"%s%s%i%s%i%s%i","fAccResX","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); } //binning Int_t Sbin=250,Tbin=250; Int_t Smin=-125,Smax=125,Tmin=-125,Tmax=125; sprintf(tmp,"%s%s%i","ShowerPlaneProj","_Sector_",sector); ShowerPlaneProj = new TH2F(tmp,title,500,-1000,1000,500,-800,1200); for(Int_t index=0;indexcd(); } void HMdcMetaAligner::fitHistograms(Int_t index) { // //Fits to a gaussian the four relevant histograms //and obtains the fit parameters for further data selection // //Till now, only working for two MDCs //Antes de hacerlo para mas, pensar que es necesario cortar! if(fNumMods==2){ Float_t XNewAreaA, XNewAreaB, YNewAreaA, YNewAreaB; Float_t S0NewAreaA, S0NewAreaB, S1NewAreaA, S1NewAreaB; Float_t XNewMeanA, XNewMeanB, YNewMeanA, YNewMeanB; Float_t S0NewMeanA, S0NewMeanB, S1NewMeanA, S1NewMeanB; TF1 *f1X = new TF1("f1X","gaus",-fXArea,fXArea); TF1 *f1Y = new TF1("f1Y","gaus",-fYArea,fYArea); TF1 *f1S = new TF1("f1S","gaus",-1,1); TF1 *totalX = new TF1("totalX","gaus(0)+pol2(3)",-fXArea,fXArea); TF1 *totalY = new TF1("totalY","gaus(0)+pol2(3)",-fYArea,fYArea); TF1 *totalS = new TF1("totalS","gaus(0)+pol2(3)",-1,1); Double_t par[6]; if(MA_DEBUG>1) cout << endl <<"**** fitHistograms() results ****" << endl; if(MA_DEBUG>1) cout << endl <<"**** Gauss fit: mean, sigma ****" << endl <<"**** Gauss+pol: mean, sigma ****" << endl; AvsBinBCS_X[index]->Fit("f1X","RQNW"); Float_t fitPar0 = f1X->GetParameter(0); // constant Float_t fitPar1 = f1X->GetParameter(1); // mean Float_t fitPar2 = f1X->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " AvsBinBCS_X[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1X->GetParameters(&par[0]); par[3] = par[4] = par[5] = 0.; totalX->SetParameters(par); AvsBinBCS_X[index]->Fit("totalX","RQNW"); fitPar0 = totalX->GetParameter(0); fitPar1 = totalX->GetParameter(1); fitPar2 = totalX->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; XNewAreaA = fXSigmas * fitPar2; XNewMeanA = fitPar1; AvsBinBCS_Y[index]->Fit("f1Y","RQNW"); fitPar0 = f1Y->GetParameter(0); // constant fitPar1 = f1Y->GetParameter(1); // mean fitPar2 = f1Y->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " AvsBinBCS_Y[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1Y->GetParameters(&par[0]); totalY->SetParameters(par); AvsBinBCS_Y[index]->Fit("totalY","RQNW"); fitPar1 = totalY->GetParameter(1); fitPar2 = totalY->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; YNewAreaA = fYSigmas * fitPar2; YNewMeanA = fitPar1; AvsBinBCS_XSlope[index]->Fit("f1S","RQNW"); fitPar0 = f1S->GetParameter(0); // constant fitPar1 = f1S->GetParameter(1); // mean fitPar2 = f1S->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " AvsBinBCS_XSlope[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1S->GetParameters(&par[0]); totalS->SetParameters(par); AvsBinBCS_XSlope[index]->Fit("totalS","RQNW"); fitPar1 = totalS->GetParameter(1); fitPar2 = totalS->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; S0NewAreaA = fS0Sigmas * fitPar2; S0NewMeanA = fitPar1; AvsBinBCS_YSlope[index]->Fit("f1S","RQNW"); fitPar0 = f1S->GetParameter(0); // constant fitPar1 = f1S->GetParameter(1); // mean fitPar2 = f1S->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " AvsBinBCS_YSlope[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1S->GetParameters(&par[0]); totalS->SetParameters(par); AvsBinBCS_YSlope[index]->Fit("totalS","RQNW"); fitPar1 = totalS->GetParameter(1); fitPar2 = totalS->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; S1NewAreaA = fS1Sigmas * fitPar2; S1NewMeanA = fitPar1; BvsAinACS_X[index]->Fit("f1X","RQNW"); fitPar0 = f1X->GetParameter(0); // constant fitPar1 = f1X->GetParameter(1); // mean fitPar2 = f1X->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " BvsAinACS_X[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1X->GetParameters(&par[0]); totalX->SetParameters(par); BvsAinACS_X[index]->Fit("totalX","RQNW"); fitPar1 = totalX->GetParameter(1); fitPar2 = totalX->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; XNewAreaB = fXSigmas * fitPar2; XNewMeanB = fitPar1; BvsAinACS_Y[index]->Fit("f1Y","RQNW"); fitPar0 = f1Y->GetParameter(0); // constant fitPar1 = f1Y->GetParameter(1); // mean fitPar2 = f1Y->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " BvsAinACS_Y[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1Y->GetParameters(&par[0]); totalY->SetParameters(par); BvsAinACS_Y[index]->Fit("totalY","RQNW"); fitPar1 = totalY->GetParameter(1); fitPar2 = totalY->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; YNewAreaB = fYSigmas * fitPar2; YNewMeanB = fitPar1; BvsAinACS_XSlope[index]->Fit("f1S","RQNW"); fitPar0 = f1S->GetParameter(0); // constant fitPar1 = f1S->GetParameter(1); // mean fitPar2 = f1S->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " BvsAinACS_XSlope[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1S->GetParameters(&par[0]); totalS->SetParameters(par); BvsAinACS_XSlope[index]->Fit("totalS","RQNW"); fitPar1 = totalS->GetParameter(1); fitPar2 = totalS->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; S0NewAreaB = fS0Sigmas * fitPar2; S0NewMeanB = fitPar1; BvsAinACS_YSlope[index]->Fit("f1S","RQNW"); fitPar0 = f1S->GetParameter(0); // constant fitPar1 = f1S->GetParameter(1); // mean fitPar2 = f1S->GetParameter(2); // sigma if(MA_DEBUG>1) cout << " BvsAinACS_YSlope[" << index << "]: " << fitPar1 << ", " << fitPar2<< ", " ; f1S->GetParameters(&par[0]); totalS->SetParameters(par); BvsAinACS_YSlope[index]->Fit("totalS","RQNW"); fitPar1 = totalS->GetParameter(1); fitPar2 = totalS->GetParameter(2); if(MA_DEBUG>1) cout << fitPar1 << ", " << fitPar2 << endl; S1NewAreaB = fS1Sigmas * fitPar2; S1NewMeanB = fitPar1; Stat_t entries = fAlignMeta->GetEntries(); HMdcHit* hitA; HMdcHit* hitB; HGeomVector projPoint; Float_t* projSlopes = new Float_t[2]; Float_t* origSlopes = new Float_t[2]; HGeomRotation rotaux; HGeomVector transaux; for (Int_t i=0;iGetEntry(i); hitA = (HMdcHit*) fHits->getMdcHitA(); hitB = (HMdcHit*) fHits->getMdcHitB(); //if(fNumMods>2) hitC = (HMdcHit*) fHits->At(2); //if(fNumMods>3) hitD = (HMdcHit*) fHits->At(3); Float_t resInAvsBinBCS_X, resInAvsBinBCS_Y; Float_t resInAvsBinBCS_XSlope, resInAvsBinBCS_YSlope; Float_t resInBvsAinACS_X, resInBvsAinACS_Y; Float_t resInBvsAinACS_XSlope, resInBvsAinACS_YSlope; //projecting on MDC B projPoint = findProjPoint(hitA,fRotMat[0],fTranslation[0],projSlopes); resInAvsBinBCS_X = hitB->getX() - projPoint.getX(); resInAvsBinBCS_Y = hitB->getY() - projPoint.getY(); origSlopes = transformToSlopes(hitB); resInAvsBinBCS_XSlope = origSlopes[0] - projSlopes[0]; resInAvsBinBCS_YSlope = origSlopes[1] - projSlopes[1]; //then projecting on MDC A rotaux = fRotMat[0].inverse(); transaux = -(rotaux * fTranslation[0]); projPoint = findProjPoint(hitB,rotaux,transaux,projSlopes); resInBvsAinACS_X = hitA->getX() - projPoint.getX(); resInBvsAinACS_Y = hitA->getY() - projPoint.getY(); origSlopes = transformToSlopes(hitA); resInBvsAinACS_XSlope = origSlopes[0] - projSlopes[0]; resInBvsAinACS_YSlope = origSlopes[1] - projSlopes[1]; //strong condition: cutting in all histograms if(MA_DEBUG>3){ cout << endl <<"Cuts in fitHistograms(): " << endl; cout << fabs(resInAvsBinBCS_X - XNewMeanA) << " vs " << XNewAreaA << endl; cout << fabs(resInAvsBinBCS_Y - YNewMeanA) << " vs " << YNewAreaA << endl; cout << fabs(resInAvsBinBCS_XSlope - S0NewMeanA) << " vs " << S0NewAreaA << endl; cout << fabs(resInAvsBinBCS_YSlope - S1NewMeanA) << " vs " << S1NewAreaA << endl; cout << fabs(resInBvsAinACS_X - XNewMeanB) << " vs " << XNewAreaB << endl; cout << fabs(resInBvsAinACS_Y - YNewMeanB) << " vs " << YNewAreaB << endl; cout << fabs(resInBvsAinACS_XSlope - S0NewMeanB) << " vs " << S0NewAreaB << endl; cout << fabs(resInBvsAinACS_YSlope - S1NewMeanB) << " vs " << S1NewAreaB << endl; } if(fabs(resInAvsBinBCS_X - XNewMeanA) < XNewAreaA && fabs(resInAvsBinBCS_Y - YNewMeanA) < YNewAreaA && fabs(resInAvsBinBCS_XSlope - S0NewMeanA) < S0NewAreaA && fabs(resInAvsBinBCS_YSlope - S1NewMeanA) < S1NewAreaA && fabs(resInBvsAinACS_X - XNewMeanB) < XNewAreaB && fabs(resInBvsAinACS_Y - YNewMeanB) < YNewAreaB && fabs(resInBvsAinACS_XSlope - S0NewMeanB) < S0NewAreaB && fabs(resInBvsAinACS_YSlope - S1NewMeanB) < S1NewAreaB ) { if(fUseSlopeCut){ // This cut makes the sample near indep. of Z translations // and results useful for X and Y minimization. // The cut is effective only in MDCB, because fTranslation // is represented in MDCB coordinates. Then, tracks passing // the cut are almost parallel to Z direction if(( fabs(hitB->getXDir()) < fSlopeCut) && ( fabs(hitB->getYDir()) < fSlopeCut)){ if(MA_DEBUG>3) cout << " -- CUT PASSED (fSlopeCut=" << fSlopeCut << " ) --" << endl; //Aqui no hago nada ??? Comprobar que realmente es lo correcto //Creo que fHits contiene todo lo que necesito tal como //viene de fAlignMeta. Solo hay que llenar el Cut si cumple el if fHits->print(); fAlignMetaCut->Fill(); fCountCut++; } } else{ if(MA_DEBUG>3) cout << " --------- CUT PASSED ------------" << endl; fAlignMetaCut->Fill(); fCountCut++; } } } } } TTree* HMdcMetaAligner::getTree(void) { // // return a tree to the minimization functional. // if(fNumMods == 1) return fAlignMeta; else { if(fUseCut) return fAlignMetaCut; return fAlignMeta; } } void HMdcMetaAligner::setTree(void) { // // Inits TNtuple // Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; fRecCount++; Char_t title[50], tmp[50], tmp2[50]; if(!fOutRoot) fOutRoot = new TFile(fNameRoot,"UPDATE"); if (fNumMods>3){ Int_t modC = fLoc[3]; Int_t modD = fLoc[4]; sprintf(title,"%s%i%s%i%s%i%s%i%s%i","Sector_",sector,"_ModA_",modA, "_ModB_",modB,"_ModC_",modC,"_ModD_",modD); sprintf(tmp,"%s%s%i%s%i%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC,"_ModD_",modD); sprintf(tmp2,"%s%s%i%s%i%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC,"_ModD_",modD); } else if(fNumMods>2){ Int_t modC = fLoc[3]; sprintf(title,"%s%i%s%i%s%i%s%i","Sector_",sector,"_ModA_",modA, "_ModB_",modB,"_ModC_",modC); sprintf(tmp,"%s%s%i%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC); sprintf(tmp2,"%s%s%i%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB,"_ModC_",modC); } else if(fNumMods>1){ sprintf(title,"%s%i%s%i%s%i","Sector_",sector, "_ModA_",modA,"_ModB_",modB); sprintf(tmp,"%s%s%i%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); sprintf(tmp2,"%s%s%i%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA,"_ModB_",modB); } else{ sprintf(title,"%s%i%s%i","Sector_",sector, "_ModA_",modA); sprintf(tmp,"%s%s%i%s%i","Meta","_Sector_",sector, "_ModA_",modA); sprintf(tmp2,"%s%s%i%s%i","MetaCut","_Sector_",sector, "_ModA_",modA); } fAlignMeta = new TTree(tmp,title); fAlignMeta->Branch("hits","HMdcMetaHit",&fHits,64000); fAlignMetaCut = new TTree(tmp2,title); fAlignMetaCut->Branch("hits","HMdcMetaHit",&fHits,64000); } Bool_t HMdcMetaAligner::init(void) { // // Inits hitaux category and calls setParContainers() // //Inicialization of hit category fHitCat=gHades->getCurrentEvent()->getCategory(catMdcHit); if (!fHitCat) { fHitCat=gHades->getSetup()->getDetector("Mdc")->buildCategory(catMdcHit); if (!fHitCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catMdcHit,fHitCat,"Mdc"); } fShowerHitCat=gHades->getCurrentEvent()->getCategory(catShowerHit); if (!fShowerHitCat) { fShowerHitCat=gHades->getSetup()->getDetector("Shower") ->buildCategory(catShowerHit); if (!fShowerHitCat) return kFALSE; else gHades->getCurrentEvent() ->addCategory(catShowerHit,fShowerHitCat,"Shower"); } fTofHitCat=gHades->getCurrentEvent()->getCategory(catTofHit); if (!fTofHitCat) { fTofHitCat=gHades->getSetup()->getDetector("Tof") ->buildCategory(catTofHit); if (!fTofHitCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catTofHit,fTofHitCat,"Tof"); } fIter1 = (HIterator *)fHitCat->MakeIterator(); fIter2 = (HIterator *)fHitCat->MakeIterator(); fIter3 = (HIterator *)fHitCat->MakeIterator(); fIter4 = (HIterator *)fHitCat->MakeIterator(); fIter5 = (HIterator *)fShowerHitCat->MakeIterator(); fIter6 = (HIterator *)fTofHitCat->MakeIterator(); setParContainers(); if(fHistoOff!=kTRUE) setHistograms(); else setTree(); return kTRUE; } Bool_t HMdcMetaAligner::reinit(void) { // // Call the functions returning auxiliar geometrical parameters // if(fAuto == kFALSE) setGeomAuxPar(); //Esta linea llamaba a la funcion cuando no debia hacerlo. //Afortunadamente no hace nada la funcion sin argumentos. //else if(fAuto == kTRUE) setGeomAuxParSim(); else if(fAuto == kTRUE) ; return kTRUE; } void HMdcMetaAligner::setGeomParAutoOn(void) { // //Turn on the automatic geometrical parameter input // fAuto =kTRUE; cout << "WARNING in HMdcMetaAligner::setGeomParManOn(): " << "introducing manually Geometrical" << endl; cout << "Parameters without containers. " << "Parameters should be in the macro" << endl; } void HMdcMetaAligner::setControlHistoOff(void) { // // Disables control histograms output (saving memory) // fHistoOff = kTRUE; } void HMdcMetaAligner::setMinimizationOff(void) { // // Disables minimization // fMinOff = kTRUE; } void HMdcMetaAligner::setUnitErrors(void) { // // introduce unit errors in Hits // fUseUnitErrors = kTRUE; } void HMdcMetaAligner::setFix(Int_t fix) { // // Fix a parameter set in minimization // fFix = fix; } void HMdcMetaAligner::setNoCut(void) { // // Disables the cuts after fitting the histograms // fUseCut = kFALSE; } void HMdcMetaAligner::setSlopeCut(Float_t SlopeCut=0.1) { // // Enables the Slope cuts after fitting the histograms. // This cut makes the sample near indep. of Z translations // and results useful for X and Y minimization. // The cut is effective only in MDCB, because fTranslation // is represented in MDCB coordinates. Then, tracks passing // the cut are almost parallel to Z direction fUseSlopeCut = kTRUE; fSlopeCut = SlopeCut; } void HMdcMetaAligner::setParContainers(void) { // // Loads the parameter containers it uses later // fMdcGeomPar=(HMdcGeomPar*)gHades->getRuntimeDb()->getContainer("MdcGeomPar"); //El SHOWER no es standard. Su geometria es hija de HShowerParSet, //el cual deriva directamente de HParSet, sin ser hijo (como en //los otros casos) de HDetGeomPar. fShowerGeometry=(HShowerGeometry*)gHades->getRuntimeDb() ->getContainer("ShowerGeometry"); if(!fShowerGeometry){ fShowerGeometry = new HShowerGeometry; if (!fShowerGeometry) { Error ("init","Unable to create Shower Geometry container"); //return kFALSE; } gHades->getRuntimeDb()->addContainer(fShowerGeometry); } fTofGeomPar=(HTofGeomPar*)gHades->getRuntimeDb()->getContainer("TofGeomPar"); } void HMdcMetaAligner::setGeomAuxPar(void) { // // Defines the parameter containers it uses later // // The transformations are: X = R X' + T // in the geometrical parameters. // To obtain relative transformations: X(B) = M X(A) + V // we follow: // X = R(A) X(A) + T(A) // X = R(B) X(B) + T(B) // // X(B) = [R(B)]E-1 R(A) X(A) + [R(B)]E-1 [T(A)-T(B)] // // M = [R(B)]E-1 R(A) // V = [R(B)]E-1 [T(A)-T(B)] // // From M it is easy to get the relative Euler angles // // CONVENTION: fRotMat[i] and fTranslation[i] store the // matrix M and vector V of the transformations made like this: // 4 modules: [0] X(D) = fRotMat[0] X(C) + fTranslation[0] // [1] X(D) = fRotMat[1] X(B) + fTranslation[1] // [2] X(D) = fRotMat[2] X(A) + fTranslation[2] // // 3 modules: [0] X(C) = fRotMat[0] X(B) + fTranslation[0] // [1] X(C) = fRotMat[1] X(A) + fTranslation[1] // // 2 modules: [0] X(B) = fRotMat[0] X(A) + fTranslation[0] // Int_t sector = fLoc[0]; Int_t moduleA = fLoc[1]; Int_t moduleB = fLoc[2]; Int_t moduleC=-1,moduleD=-1; HGeomVector euler; HGeomTransform transformA; transformA = fMdcGeomPar->getModule(sector,moduleA)->getLabTransform(); HGeomRotation rotA; rotA = transformA.getRotMatrix(); HGeomVector vectorA; vectorA = transformA.getTransVector(); HGeomRotation rotB,rotC,rotD; HGeomVector vectorB,vectorC, vectorD; if(fNumMods>3){ moduleC = fLoc[3]; moduleD = fLoc[4]; HGeomTransform transformB; transformB = fMdcGeomPar->getModule(sector,moduleB)->getLabTransform(); HGeomTransform transformC; transformC = fMdcGeomPar->getModule(sector,moduleC)->getLabTransform(); HGeomTransform transformD; transformD = fMdcGeomPar->getModule(sector,moduleD)->getLabTransform(); rotB = transformB.getRotMatrix(); vectorB = transformB.getTransVector(); rotC = transformC.getRotMatrix(); vectorC = transformC.getTransVector(); rotD = transformD.getRotMatrix(); vectorD = transformD.getTransVector(); if(MA_DEBUG>0){ cout << endl <<"Debugging output in HMdcMetaAligner::setGeomAuxPar" << endl; cout << "Original transformation from container" << endl; cout << " ------ SECTOR " << sector << " ------ " << endl; cout << "MDC A (Module " << moduleA << ")" << endl; transformA.print(); cout << "MDC B (Module " << moduleB << ")" << endl; transformB.print(); cout << "MDC C (Module " << moduleC << ")" << endl; transformC.print(); cout << "MDC D (Module " << moduleD << ")" << endl; transformD.print(); } //From the previous transformation, get the relative transformation // M = [R(D)]E-1 R(C) // V = [R(D)]E-1 [T(C)-T(D)] HGeomRotation rotDinv = rotD.inverse(); HGeomRotation relrot = rotDinv * rotC; HGeomVector relvector = rotDinv * (vectorC - vectorD); if(MA_DEBUG>0){ //Printing relative transformation (MDCC -> MDCD) cout << endl <<"Relative transformation: (MDCC -> MDCD)" << endl; relrot.print(); relvector.print(); } //Filling the first rotation (MDCC and MDCD) euler = findEulerAngles(relrot); fillRotMatrix(0,euler.getX(),euler.getY(),euler.getZ()); fillTranslation(0,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(0,euler.getX(),euler.getY(),euler.getZ()); //The same for (MDCB and MDCD) relrot = rotDinv * rotB; relvector = rotDinv * (vectorB - vectorD); if(MA_DEBUG>0){ //Printing relative transformation (MDCB -> MDCD) cout << endl <<"Relative transformation: (MDCB -> MDCD)" << endl; relrot.print(); relvector.print(); } //Filling the second rotation (MDCB and MDCD) euler = findEulerAngles(relrot); fillRotMatrix(1,euler(0),euler(1),euler(2)); fillTranslation(1,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(1,euler(0),euler(1),euler(2)); //The same for (MDCA and MDCD) relrot = rotDinv * rotA; relvector = rotDinv * (vectorA - vectorD); if(MA_DEBUG>0){ //Printing relative transformation (MDCA -> MDCD) cout << endl <<"Relative transformation: (MDCA -> MDCD)" << endl; relrot.print(); relvector.print(); } //Filling the second rotation (MDCA and MDCD) euler = findEulerAngles(relrot); fillRotMatrix(2,euler(0),euler(1),euler(2)); fillTranslation(2,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(2,euler(0),euler(1),euler(2)); cout << "**************************************************" << endl; cout << "* HMdcMetaAligner::setGeomAuxPar: from geom params: *" << endl; cout << "* Sector: "<< sector << " ModA: " << moduleA << " ModB: " << moduleB << " ModC: " << moduleC << " ModD: " << moduleD << endl; for(Int_t c=0;c<3;c++){ cout << "* Rotation(" << c << "): " << fEuler[c].getX() << ", " << fEuler[c].getY() << ", " << fEuler[c].getZ() << endl; cout << "* Translation(" << c << "): " << fTranslation[c].getX() << ", " << fTranslation[c].getY() << ", " << fTranslation[c].getZ() << endl; } cout << "**************************************************" << endl; } else if(fNumMods>2){ moduleC = fLoc[3]; HGeomTransform transformB; transformB = fMdcGeomPar->getModule(sector,moduleB)->getLabTransform(); HGeomTransform transformC; transformC = fMdcGeomPar->getModule(sector,moduleC)->getLabTransform(); rotB = transformB.getRotMatrix(); vectorB = transformB.getTransVector(); rotC = transformC.getRotMatrix(); vectorC = transformC.getTransVector(); if(MA_DEBUG>0){ cout << endl <<"Debugging output in HMdcMetaAligner::setGeomAuxPar" << endl; cout << "Original transformation from container" << endl; cout << " ------ SECTOR " << sector << " ------ " << endl; cout << "MDC A (Module " << moduleA << ")" << endl; transformA.print(); cout << "MDC B (Module " << moduleB << ")" << endl; transformB.print(); cout << "MDC C (Module " << moduleC << ")" << endl; transformC.print(); } //From the previous transformation, get the relative transformation // M = [R(B)]E-1 R(A) // V = [R(B)]E-1 [T(A)-T(B)] HGeomRotation rotCinv = rotC.inverse(); HGeomRotation relrot = rotCinv * rotB; HGeomVector relvector = rotCinv * (vectorB - vectorC); if(MA_DEBUG>0){ //Printing relative transformation (MDCB -> MDCC) cout << endl <<"Relative transformation: (MDCB -> MDCC)" << endl; relrot.print(); relvector.print(); } //Filling the second rotation (MDCB and MDCC) euler = findEulerAngles(relrot); fillRotMatrix(0,euler(0),euler(1),euler(2)); fillTranslation(0,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(0,euler(0),euler(1),euler(2)); //The same for (MDCA and MDCC) relrot = rotCinv * rotA; relvector = rotCinv * (vectorA - vectorC); if(MA_DEBUG>0){ //Printing relative transformation (MDCA -> MDCC) cout << endl <<"Relative transformation: (MDCA -> MDCC)" << endl; relrot.print(); relvector.print(); } //Filling the second rotation (MDCA and MDCC) euler = findEulerAngles(relrot); fillRotMatrix(1,euler(0),euler(1),euler(2)); fillTranslation(1,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(1,euler(0),euler(1),euler(2)); cout << "**************************************************" << endl; cout << "* HMdcMetaAligner::setGeomAuxPar: from geom params: *" << endl; cout << "* Sector: "<< sector << " ModA: " << moduleA << " ModB: " << moduleB << " ModC: " << moduleC << endl; for(Int_t c=0;c<2;c++){ cout << "* Rotation(" << c << "): " << fEuler[c].getX() << ", " << fEuler[c].getY() << ", " << fEuler[c].getZ() << endl; cout << "* Translation(" << c << "): " << fTranslation[c].getX() << ", " << fTranslation[c].getY() << ", " << fTranslation[c].getZ() << endl; } cout << "**************************************************" << endl; } else if(fNumMods>1){ HGeomTransform transformB; transformB = fMdcGeomPar->getModule(sector,moduleB)->getLabTransform(); rotB = transformB.getRotMatrix(); vectorB = transformB.getTransVector(); if(MA_DEBUG>0){ cout << endl <<"Debugging output in HMdcMetaAligner::setGeomAuxPar" << endl; cout << "Original transformation from container" << endl; cout << " ------ SECTOR " << sector << " ------ " << endl; cout << "MDC A (Module " << moduleA << ")" << endl; transformA.print(); cout << "MDC B (Module " << moduleB << ")" << endl; transformB.print(); } //From the previous transformation, get the relative transformation // M = [R(B)]E-1 R(A) // V = [R(B)]E-1 [T(A)-T(B)] HGeomRotation rotBinv = rotB.inverse(); HGeomRotation relrot = rotBinv * rotA; HGeomVector relvector = rotBinv * (vectorA - vectorB); if(MA_DEBUG>0){ //Printing relative transformation cout << endl <<"Relative transformation: " << endl; relrot.print(); relvector.print(); } //Filling the first rotation (MDCA and MDCB) euler = findEulerAngles(relrot); fillRotMatrix(0,euler(0),euler(1),euler(2)); fillTranslation(0,relvector.getX(),relvector.getY(),relvector.getZ()); setEulerAngles(0,euler(0),euler(1),euler(2)); cout << "**************************************************" << endl; cout << "* HMdcMetaAligner::setGeomAuxPar: from geom params: *" << endl; cout << "* Sector: "<< sector << " ModA: " << moduleA << " ModB: " << moduleB << endl; cout << "* Rotation(0): " << fEuler[0].getX() << ", " << fEuler[0].getY() << ", " << fEuler[0].getZ() << endl; cout << "* Translation: " << relvector.getX() << ", " << relvector.getY() << " , " << relvector.getZ() << endl; cout << "**************************************************" << endl; } else{//1MDC if(MA_DEBUG>0){ cout << endl <<"Debugging output in HMdcMetaAligner::setGeomAuxPar" << endl; cout << "Original transformation from container" << endl; cout << " ------ SECTOR " << sector << " ------ " << endl; cout << "MDC A (Module " << moduleA << ")" << endl; transformA.print(); cout << " NO RELATIVE TRANSFORMATION (ONLY ONE MDC) " << endl; } } // Relative transformation of Shower vs. MDCs // The transformations are: X(SHOWER) = Rrel X(MDC) + Trel // where // X(SHOWER) = R(SHOWER)E-1 (X(LAB) - T(SHOWER)) // X(LAB) = R(MDC) X(MDC) + T(MDC) // // X(SHOWER)=[R(SHOWER)]E-1 R(MDC) X(MDC)+[R(SHOWER)]E-1 [T(MDC)-T(SHOWER)] // // Rrel = [R(SHOWER)]E-1 R(MDC) // Trel = [R(SHOWER)]E-1 [T(MDC)-T(SHOWER)] // // From M it is easy to get the relative Euler angles HGeomTransform transformS = fShowerGeometry->getTransform(sector,0); HGeomRotation rotS; rotS = transformS.getRotMatrix(); HGeomVector vectorS; vectorS = transformS.getTransVector(); HGeomTransform transformT; HGeomRotation rotT[8]; HGeomVector vectorT[8]; for(Int_t module=0;module<8;module++) { transformT = fTofGeomPar ->getModule(sector,module)->getLabTransform(); rotT[module] = transformT.getRotMatrix(); vectorT[module] = transformT.getTransVector(); } HGeomRotation rotTinv; HGeomRotation rotSinv = rotS.inverse(); HGeomRotation relSrot,relTrot; HGeomVector relSvector,relTvector; if (fNumMods>3){ relSrot = rotSinv * rotD; relSvector = rotSinv * (vectorD - vectorS); } else if(fNumMods>2){ relSrot = rotSinv * rotC; relSvector = rotSinv * (vectorC - vectorS); } else if(fNumMods>1){ relSrot = rotSinv * rotB; relSvector = rotSinv * (vectorB - vectorS); } else{ //1 MDC relSrot = rotSinv * rotA; relSvector = rotSinv * (vectorA - vectorS); } for(Int_t module=0;module<8;module++) { rotTinv = rotT[module].inverse(); if (fNumMods>3){ relTrot = rotTinv * rotD; relTvector = rotTinv * (vectorD - vectorT[module]); } else if(fNumMods>2){ relTrot = rotTinv * rotC; relTvector = rotTinv * (vectorC - vectorT[module]); } else if(fNumMods>2){ relTrot = rotTinv * rotB; relTvector = rotTinv * (vectorB - vectorT[module]); } else{ relTrot = rotTinv * rotA; relTvector = rotTinv * (vectorA - vectorT[module]); } euler = findEulerAngles(relTrot); fillMetaRotMatrix(module,euler(0),euler(1),euler(2)); fillMetaTranslation(module,relTvector.getX(), relTvector.getY(),relTvector.getZ()); setMetaEulerAngles(module,euler(0),euler(1),euler(2)); } //Filling the Shower rotation euler = findEulerAngles(relSrot); fillMetaRotMatrix(8,euler(0),euler(1),euler(2)); fillMetaTranslation(8,relSvector.getX(),relSvector.getY(),relSvector.getZ()); setMetaEulerAngles(8,euler(0),euler(1),euler(2)); cout << "**************************************************" << endl; cout << "* HMdcMetaAligner::setGeomAuxPar: from geom params: *" << endl; cout << "* Sector: "<< sector << " Meta vs. "; if(fNumMods==4) cout << " MDC module: " << moduleD << endl; if(fNumMods==3) cout << " MDC module: " << moduleC << endl; if(fNumMods==2) cout << " MDC module: " << moduleB << endl; else cout << " MDC module: " << moduleA << endl; cout << "From 0 to 7: TOF modules. 8 corresponds to SHOWER" << endl;; for(Int_t module=0;module<9;module++){ cout << module << " Rotation: " << fMEuler[module].getX() << ", " << fMEuler[module].getY() << ", " << fMEuler[module].getZ() << endl; cout << module << " Translation: " << fMTranslation[module].getX() << ", " << fMTranslation[module].getY() << " , " << fMTranslation[module].getZ() << endl; cout << "**************************************************" << endl; } } HGeomVector HMdcMetaAligner::findEulerAngles(HGeomRotation rot){ // // From the relative rotation, get the euler angles (radians) // // From an Euler rotation (see Dahlinger paper for the convention) // the euler angles are: // euler[0] = atan2(rot(5)/sin(euler[1]),rot(2)/sin(euler[1])) // euler[1] = acos (rot(8)) with possible change of sign // euler[2] = atan2(rot(7)/sin(euler[1]),-rot(6)/sin(euler[1])) // Double_t euler[3]; HGeomVector eul; //Checking if rot(8) is in the acos() domain if(rot(8)< 0.99999 && rot(8)>-0.99999) euler[1] = acos(rot(8)); else euler[1]=0; Double_t sinaux; if(euler[1] == 0.){ //euler[0] and euler[2] are equivalent. Write all in euler[0] euler[0]= (TMath::Pi()/2)+acos(rot(0)); euler[2]=-(TMath::Pi()/2); } else{ //IS AN EULER MATRIX sinaux = sin(euler[1]); euler[0] = atan2(rot(5)/sinaux,rot(2)/sinaux); euler[2] = atan2(rot(7)/sinaux,-rot(6)/sinaux); } if(MA_DEBUG>0){ cout << endl <<"Euler angles: " << euler[0] << ", " << euler[1] << ", " << euler[2] << endl; } HGeomRotation test; test.setEulerAngles(euler[0]*180./TMath::Pi(), euler[1]*180./TMath::Pi(), euler[2]*180./TMath::Pi()); if(MA_DEBUG>0){ cout << endl <<"Rotation from Euler angles (first attemp): " << endl; test.print(); } //Now solving the problem when euler[1]<0 Double_t eps = 0.0001; //check precision if( (fabs(test(0)-rot(0))>eps) || (fabs(test(1)-rot(1))>eps) || (fabs(test(3)-rot(3))>eps) || (fabs(test(4)-rot(4))>eps) ) { if(MA_DEBUG>0) cout << endl << "Bad election in first euler angle! " << "Trying again. "<< endl; euler[1] = - euler[1]; sinaux = sin(euler[1]); euler[0] = atan2(rot(5)/sinaux,rot(2)/sinaux); euler[2] = atan2(rot(7)/sinaux,-rot(6)/sinaux); test.setEulerAngles(euler[0]*180./TMath::Pi(), euler[1]*180./TMath::Pi(), euler[2]*180./TMath::Pi()); if(MA_DEBUG>0){ cout << "Rotation from Euler angles (second attemp): " << endl; test.print(); } //testing if euler angles are rigth now if( (fabs(test(0)-rot(0))>eps) || (fabs(test(1)-rot(1))>eps) || (fabs(test(3)-rot(3))>eps) || (fabs(test(4)-rot(4))>eps) ){ cout << endl << "FATAL ERROR: Bad election in euler angles! "<< endl; cout << "Original rot matrix: "<< endl; rot.print(); cout << "From obtained euler angles: " << endl; test.print(); //What to do?? } } eul.setX(euler[0]); eul.setY(euler[1]); eul.setZ(euler[2]); return eul; } void HMdcMetaAligner::setGeomAuxParSim(Int_t ind=0, Float_t eu1=0, Float_t eu2=0, Float_t eu3=0, Float_t tr1=0, Float_t tr2=0, Float_t tr3=0){ // // Entering geometrical parameters. // // To be used introducing in the macro the parameters for // testing, optimization ... // cout << "WARNING: Introducing automatically Geometrical" << endl; cout << "Parameters without containers" << endl; if(ind!=0 || eu1!=0 || eu2!=0 || eu3!=0 || tr1!=0 || tr2!=0 || tr3!=0){ fillRotMatrix(ind,eu1,eu2,eu3); fillTranslation(ind,tr1,tr2,tr3); setEulerAngles(ind,eu1,eu2,eu3); if(MA_DEBUG>0){ cout << "Transformation[" << ind << "]:" << endl; cout <<" Euler angles: " << eu1 << ", " << eu2 << ", " << eu3 << endl; cout << " Translation: " << tr1 << ", " << tr2 << ", " << tr3 << endl; } } else{ // Parameters could enter automatically? cout << " !! Sorry, not implemented !! " << endl; } } void HMdcMetaAligner::setMetaGeomAuxParSim(Int_t ind=0, Float_t eu1=0, Float_t eu2=0, Float_t eu3=0, Float_t tr1=0, Float_t tr2=0, Float_t tr3=0) { // // Entering geometrical parameters. // // To be used introducing in the macro the parameters for // testing, optimization ... // cout << "WARNING: Introducing automatically Geometrical" << endl; cout << "Parameters without containers" << endl; if(ind!=0 || eu1!=0 || eu2!=0 || eu3!=0 || tr1!=0 || tr2!=0 || tr3!=0){ fillMetaRotMatrix(ind,eu1,eu2,eu3); fillMetaTranslation(ind,tr1,tr2,tr3); setMetaEulerAngles(ind,eu1,eu2,eu3); if(MA_DEBUG>0){ cout << "MetaTransformation[" << ind << "]:" << endl; cout <<" Euler angles: " << eu1 << ", " << eu2 << ", " << eu3 << endl; cout << " Translation: " << tr1 << ", " << tr2 << ", " << tr3 << endl; } } else{ // Parameters could enter automatically? cout << " !! Sorry, not implemented !! " << endl; } } Int_t HMdcMetaAligner::execute(void) { // // Iteration in the hit category. Fills histograms // and TTree and calculates relevant quantities // if(fNumMods==1) execute1(); if(fNumMods==2) execute2(); if(fNumMods==3) execute3(); if(fNumMods==4) execute4(); return 0; } void HMdcMetaAligner::execute4(void) { // // New execute for four MDCs // Find the projection of Hits in the most external MDCs // on the inners (besides this can be modified, changing the // order of the MDCs in the constructor). Then merge the hits // and continue to find till the last MDC. If successfull, // fills a TClonesArray of Hits in a TTree for further // analysis and minimization // HMdcHit* pHitA; HMdcHit* pHitB; HMdcHit* pHitC; HMdcHit* pHitD; HMdcHit* pHitCD; HMdcHit* pHitBCD; HMdcHit* pHitABCD; HShowerHit* pHitS; HTofHit* pHitT; HGeomVector* lSVec = new HGeomVector; HGeomVector* lTVec = new HGeomVector; Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; Int_t modC = fLoc[3]; Int_t modD = fLoc[4]; HLocation local; local.setNIndex(2); local.set(2,sector,modD); HGeomVector projPoint; HGeomVector projMetaPoint; Float_t* projSlopes = new Float_t[2]; Bool_t usedA = kFALSE; Bool_t usedB = kFALSE; Bool_t usedC = kFALSE; Bool_t invalidA = kFALSE; Bool_t invalidB = kFALSE; Bool_t invalidC = kFALSE; Bool_t SUsed = kFALSE; Bool_t SInvalid = kFALSE; Bool_t TUsed = kFALSE; Bool_t TInvalid = kFALSE; Int_t detnum=-1; Int_t modTof =0; fNEntries++; fIter1->Reset(); fIter1->gotoLocation(local); while ( (pHitD =(HMdcHit*)fIter1->Next()) != 0 ){ fHitsMdc[0]++; //Hits found in MDC D usedA = kFALSE; //reinit control flags usedB = kFALSE; usedC = kFALSE; invalidA = kFALSE; invalidB = kFALSE; invalidC = kFALSE; if(MA_DEBUG>3) { cout << " ----- SECTOR " << sector << " -----"<< endl; cout << "Module " << modD << ", fHitsMdc " << fHitsMdc[0] << endl; } // Calculation of cross point and slopes in MDC C // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(C) = fRotMat[0]-1 * (X(D) - fTranslation[0]) // ==> X(C) = fRotMat[0]-1 * X(D) - fRotMat[0]-1 * fTranslation[0] // because we are calculating the projection of MDC D Hits on MDC C // HGeomRotation rotInv = fRotMat[0].inverse(); HGeomVector trasInv = -(rotInv * fTranslation[0]); projPoint = findProjPoint(pHitD,rotInv,trasInv,projSlopes); //Iteration on the second MDC (MDCC) for each hit in the first (MDCD) local.set(2,sector,modC); fIter2->Reset(); fIter2->gotoLocation(local); while ( (pHitC =(HMdcHit*)fIter2->Next()) != 0 ){ //hits found in MDCC, provided there is a Hit in MDCD fHitsMdc[1]++; if(MA_DEBUG>3) cout << "Module " << modC << ", fHitsMdc " << fHitsMdc[1] << endl; if(isInsideWindow(1,pHitC,projPoint,projSlopes)&& (invalidB!=kTRUE)&&(invalidA!=kTRUE)){ if(usedC == kFALSE){ usedC = kTRUE; // MDCC hits found in window (only used if unique) fHitsFoundInWindow[0]++; fHitsFoundAndUnique[0]++; //merging hits of MDCC and MDCD on MDCD coordinate system pHitCD = mergeHits(pHitD,pHitC,fRotMat[0],fTranslation[0]); // Calculation of cross point and slopes in MDC B // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(B) = fRotMat[1]-1 * (X(D) - fTranslation[1]) // ==> X(B) = fRotMat[1]-1 * X(D) - fRotMat[1]-1 * fTranslation[1] // because we are calculating the projection of MDCD Hits on MDCB // (in this case the merging of MDCD and MDCC, but in MDCD // coordinate system anyway) rotInv = fRotMat[1].inverse(); trasInv = -(rotInv * fTranslation[1]); projPoint = findProjPoint(pHitCD,rotInv,trasInv,projSlopes); //Iteration on the third MDC (MDCB) for each couple MDCD-MDCC local.set(2,sector,modB); fIter3->Reset(); fIter3->gotoLocation(local); while ( (pHitB =(HMdcHit*)fIter3->Next()) != 0 ){ //hits found in MDCB, provided there is a Hit in MDCD //and one in MDC. If there is more than one in MDCC the //number is also increased but the hit will not be used fHitsMdc[2]++; if(MA_DEBUG>3) cout << "Module " << modB << ", fHitsMdc " << fHitsMdc[2] << endl; if(isInsideWindow(0,pHitB,projPoint,projSlopes)&&(invalidA!=kTRUE)){ if(usedB == kFALSE){ usedB = kTRUE; // MDCB hits found in window (only used if unique) fHitsFoundInWindow[1]++; fHitsFoundAndUnique[1]++; //merging hits of MDCB and MDCC-MDCD on MDCD coordinate system pHitBCD = mergeHits(pHitCD,pHitB,fRotMat[1],fTranslation[1]); // Calculation of cross point and slopes in MDC A // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(A) = fRotMat[2]-1 * (X(D) - fTranslation[2]) // ==> X(A)=fRotMat[2]-1*X(D)-fRotMat[2]-1*fTranslation[2] // because we are calculating the projection of MDCD Hits // on MDCA (in this case the merging of MDCD, MDCC and // MDCB, but in MDCD coordinate system anyway) rotInv = fRotMat[2].inverse(); trasInv = -(rotInv * fTranslation[2]); projPoint = findProjPoint(pHitBCD,rotInv,trasInv,projSlopes); //Iteration on the fourth MDC (MDCA) for each MDCD-MDCC-MDCB local.set(2,sector,modA); fIter4->Reset(); fIter4->gotoLocation(local); while ( (pHitA =(HMdcHit*)fIter4->Next()) != 0 ){ //hits found in MDCA, provided there is a Hit in MDCD //MDC and MDCB. If there is more than one in MDCC or MDCB //the number is also increased but the hit will not be used fHitsMdc[3]++; if(MA_DEBUG>3) cout << "Module " << modA << ", fHitsMdc " << fHitsMdc[3] << endl; if(isInsideWindow(0,pHitA,projPoint,projSlopes)){ if(usedA == kFALSE){ usedA = kTRUE; // MDCA hits found in window (only used if unique) fHitsFoundInWindow[2]++; fHitsFoundAndUnique[2]++; //Real number of matched hits fCount++; //merging all HITS (curiosity or can be used??) pHitABCD = mergeHits(pHitBCD,pHitA,fRotMat[2], fTranslation[2]); //Filling the TClonesArray for storage in TTree //Will be used only if //Now,we enter in the META part! //projectin on SHOWER plane projMetaPoint = findProjMetaPoint(pHitC,pHitD,fRotMat[0], fTranslation[0], fMRotMat[8], fMTranslation[8]); detnum =-1; //SHOWER iteration fIter5->Reset(); while ( (pHitS =(HShowerHit*)fIter5->Next()) != 0 ){ if(isInsideMetaWindow(sector,pHitS,projMetaPoint,lSVec)){ //SHOWER Hit found in window if(SUsed == kFALSE){ SUsed = kTRUE; fHitsFoundInShowerWindow++; fHitsFoundInShowerAndUnique++; fSCount++; detnum =8; //... } else{ //that is, if SUsed == kTRUE if(SInvalid == kFALSE){ fSCount--; SInvalid = kTRUE; fHitsFoundInShowerAndUnique--; fSDiscart++; detnum =-1; } } } } //TOF iteration fIter6->Reset(); while ( (pHitT =(HTofHit*)fIter6->Next()) != 0 ){ modTof = pHitT->getModule(); projMetaPoint = findProjMetaPoint(pHitC,pHitD,fRotMat[0], fTranslation[0], fMRotMat[modTof], fMTranslation[modTof]); if(isInsideMetaWindow(sector,pHitT,projMetaPoint,lTVec)){ //TOF Hit found in window if(TUsed == kFALSE){ TUsed = kTRUE; fHitsFoundInTofWindow[modTof]++; fHitsFoundInTofAndUnique[modTof]++; fTCount++; if(detnum == 8) detnum=detnum+1+modTof; else detnum = modTof; //... } else{ //that is, if TUsed == kTRUE //POR AHORA NO PERMITO DOS HITS EN TOF.CAMBIAR!! if(TInvalid == kFALSE){ fTCount--; TInvalid = kTRUE; fHitsFoundInTofAndUnique[modTof]--; fTDiscart[modTof]++; detnum = -1; } } } } fHits->setMdcHitA(pHitA); fHits->setMdcHitB(pHitB); fHits->setMdcHitC(pHitC); fHits->setMdcHitD(pHitD); if(detnum<0) fHits->setNumberOfMetaHits(0); else { fHits->setNumberOfMetaHits(1); if(detnum>7) { fHits->setNumberOfMetaHits(2); fHits->setLocalHitA(lSVec); fHits->setLocalHitB(lTVec); } else if(detnum==8) fHits->setLocalHitA(lSVec); else fHits->setLocalHitA(lTVec); } fHits->setDetector(detnum); } // End of if(usedA == kFALSE){ else{ //that is, if usedA == kTRUE //More than one candidate on window!! Discart if(invalidA == kFALSE){ fCount--; invalidA = kTRUE; fDiscart[2]++; fHitsFoundAndUnique[2]--; } } // End of else{ //that is, if usedA == kTRUE } // End of if(isInsideWindow(pHitA,projPoint,projSlopes)){ } // End of while( (pHitA =(HMdcHit*)fIter4->Next()) != 0 ){ } // End of if(usedB == kFALSE){ else{ //that is, if usedB == kTRUE //More than one candidate on window!! Discart if(invalidB == kFALSE){ invalidB = kTRUE; fDiscart[1]++; fHitsFoundAndUnique[1]--; } } // End of else{ //that is, if usedB == kTRUE } // End of if(isInsideWindow(pHitB,projPoint,projSlopes)){ } // End of while( (pHitB =(HMdcHit*)fIter3->Next()) != 0 ){ } // End of if(usedC == kFALSE){ else{ //that is, if usedC == kTRUE //More than one candidate on window!! Discart if(invalidC == kFALSE){ invalidC = kTRUE; fDiscart[0]++; fHitsFoundAndUnique[0]--; } } // End of else{ //that is, if usedC == kTRUE } // End of if(isInsideWindow(pHitC,projPoint,projSlopes)){ } // End of while( (pHitC =(HMdcHit*)fIter2->Next()) != 0 ){ if(usedC == kTRUE && invalidC != kTRUE && usedB == kTRUE && invalidB != kTRUE && usedA == kTRUE && invalidA != kTRUE){ fAlignMeta->Fill(); if(fCount%100 ==0) cout << "."<< flush; } } // End of while ( (pHitD =(HMdcHit*)fIter1->Next()) != 0 ){ } void HMdcMetaAligner::execute3(void) { // // New execute for more than two MDCs // HMdcHit* pHitA; HMdcHit* pHitB; HMdcHit* pHitC; HMdcHit* pHitBC; HMdcHit* pHitABC; HShowerHit* pHitS; HTofHit* pHitT; HGeomVector* lSVec = new HGeomVector; HGeomVector* lTVec = new HGeomVector; Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; Int_t modC = fLoc[3]; HLocation local; local.setNIndex(2); local.set(2,sector,modC); HGeomVector projPoint; HGeomVector projMetaPoint; Float_t* projSlopes = new Float_t[2]; Bool_t usedA = kFALSE; Bool_t usedB = kFALSE; Bool_t invalidA = kFALSE; Bool_t invalidB = kFALSE; Bool_t SUsed = kFALSE; Bool_t SInvalid = kFALSE; Bool_t TUsed = kFALSE; Bool_t TInvalid = kFALSE; Int_t detnum=-1; Int_t modTof =0; fNEntries++; fIter1->Reset(); fIter1->gotoLocation(local); while ( (pHitC =(HMdcHit*)fIter1->Next()) != 0 ){ fHitsMdc[0]++; //Hits found in MDC C usedA = kFALSE; //reinit control flags usedB = kFALSE; invalidA = kFALSE; invalidB = kFALSE; if(MA_DEBUG>3) { cout << " ----- SECTOR " << sector << " -----"<< endl; cout << "Module " << modC << ", fHitsMdc " << fHitsMdc[0] << endl; } // Calculation of cross point and slopes in MDC B // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(B) = fRotMat[0]-1 * (X(C) - fTranslation[0]) // ==> X(B) = fRotMat[0]-1 * X(C) - fRotMat[0]-1 * fTranslation[0] // because we are calculating the projection of MDC C Hits on MDC B // HGeomRotation rotInv = fRotMat[0].inverse(); HGeomVector trasInv = -(rotInv * fTranslation[0]); projPoint = findProjPoint(pHitC,rotInv,trasInv,projSlopes); //Iteration on the second MDC (MDCB) for each hit in the first (MDCC) local.set(2,sector,modB); fIter2->Reset(); fIter2->gotoLocation(local); while ( (pHitB =(HMdcHit*)fIter2->Next()) != 0 ){ //hits found in MDCB, provided there is a Hit in MDCC fHitsMdc[1]++; if(MA_DEBUG>3) cout << "Module " << modB << ", fHitsMdc " << fHitsMdc[1] << endl; if(isInsideWindow(1,pHitB,projPoint,projSlopes)&&(invalidA!=kTRUE)){ if(usedB == kFALSE){ usedB = kTRUE; // MDCB hits found in window (only used if unique) fHitsFoundInWindow[0]++; fHitsFoundAndUnique[0]++; //merging hits of MDCB and MDCC on MDCC coordinate system pHitBC = mergeHits(pHitC,pHitB,fRotMat[0],fTranslation[0]); // Calculation of cross point and slopes in MDC A // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(A) = fRotMat[1]-1 * (X(C) - fTranslation[1]) // ==> X(A) = fRotMat[1]-1 * X(C) - fRotMat[1]-1 * fTranslation[1] // because we are calculating the projection of MDCC Hits on MDCA // (in this case the merging of MDCC and MDCB, but in MDCC // coordinate system anyway) rotInv = fRotMat[1].inverse(); trasInv = -(rotInv * fTranslation[1]); projPoint = findProjPoint(pHitBC,rotInv,trasInv,projSlopes); //Iteration on the third MDC (MDCA) for each couple MDCC-MDCB local.set(2,sector,modA); fIter3->Reset(); fIter3->gotoLocation(local); while ( (pHitA =(HMdcHit*)fIter3->Next()) != 0 ){ //hits found in MDCA, provided there is a Hit in MDCC //and one in MDB. If there is more than one in MDCB the //number is also increased but the hit will not be used fHitsMdc[2]++; if(MA_DEBUG>3) cout << "Module " << modB << ", fHitsMdc " << fHitsMdc[1] << endl; if(isInsideWindow(0,pHitA,projPoint,projSlopes)){ if(usedA == kFALSE){ usedA = kTRUE; // MDCA hits found in window (only used if unique) fHitsFoundInWindow[1]++; fHitsFoundAndUnique[1]++; //Real number of matched hits fCount++; //merging all HITS (curiosity or can be used??) pHitABC = mergeHits(pHitBC,pHitA,fRotMat[1], fTranslation[1]); //Now,we enter in the META part! //projectin on SHOWER plane projMetaPoint = findProjMetaPoint(pHitB,pHitC,fRotMat[0], fTranslation[0], fMRotMat[8], fMTranslation[8]); detnum =-1; //SHOWER iteration fIter5->Reset(); while ( (pHitS =(HShowerHit*)fIter5->Next()) != 0 ){ if(isInsideMetaWindow(sector,pHitS,projMetaPoint,lSVec)){ //SHOWER Hit found in window if(SUsed == kFALSE){ SUsed = kTRUE; fHitsFoundInShowerWindow++; fHitsFoundInShowerAndUnique++; fSCount++; detnum =8; //... } else{ //that is, if SUsed == kTRUE if(SInvalid == kFALSE){ fSCount--; SInvalid = kTRUE; fHitsFoundInShowerAndUnique--; fSDiscart++; detnum =-1; } } } } //TOF iteration fIter6->Reset(); while ( (pHitT =(HTofHit*)fIter6->Next()) != 0 ){ modTof = pHitT->getModule(); projMetaPoint = findProjMetaPoint(pHitB,pHitC,fRotMat[0], fTranslation[0], fMRotMat[modTof], fMTranslation[modTof]); if(isInsideMetaWindow(sector,pHitT,projMetaPoint,lTVec)){ //TOF Hit found in window if(TUsed == kFALSE){ TUsed = kTRUE; fHitsFoundInTofWindow[modTof]++; fHitsFoundInTofAndUnique[modTof]++; fTCount++; if(detnum == 8) detnum=detnum+1+modTof; else detnum = modTof; //... } else{ //that is, if TUsed == kTRUE //POR AHORA NO PERMITO DOS HITS EN TOF.CAMBIAR!! if(TInvalid == kFALSE){ fTCount--; TInvalid = kTRUE; fHitsFoundInTofAndUnique[modTof]--; fTDiscart[modTof]++; detnum = -1; } } } } fHits->setMdcHitA(pHitA); fHits->setMdcHitB(pHitB); fHits->setMdcHitC(pHitC); if(detnum<0) fHits->setNumberOfMetaHits(0); else { fHits->setNumberOfMetaHits(1); if(detnum>7) { fHits->setNumberOfMetaHits(2); fHits->setLocalHitA(lSVec); fHits->setLocalHitB(lTVec); } else if(detnum==8) fHits->setLocalHitA(lSVec); else fHits->setLocalHitA(lTVec); } fHits->setDetector(detnum); } // End of if(usedA == kFALSE){ else{ //that is, if usedA == kTRUE //More than one candidate on window!! Discart if(invalidA == kFALSE){ fCount--; invalidA = kTRUE; fDiscart[1]++; fHitsFoundAndUnique[1]--; } } // End of else{ //that is, if usedA == kTRUE } // End of if(isInsideWindow(pHitA,projPoint,projSlopes)){ } // End of while( (pHitA =(HMdcHit*)fIter3->Next()) != 0 ){ } // End of if(usedB == kFALSE){ else{ //that is, if usedB == kTRUE //More than one candidate on window!! Discart if(invalidB == kFALSE){ invalidB = kTRUE; fDiscart[0]++; fHitsFoundAndUnique[0]--; } } // End of else{ //that is, if usedB == kTRUE } // End of if(isInsideWindow(pHitB,projPoint,projSlopes)&&...){ } // End of while( (pHitB =(HMdcHit*)fIter2->Next()) != 0 ){ if(usedB == kTRUE && invalidB != kTRUE && usedA == kTRUE && invalidA != kTRUE){ fAlignMeta->Fill(); if(fCount%100 ==0) cout << "."<< flush; } } // End of while ( (pHitC =(HMdcHit*)fIter1->Next()) != 0 ){ } void HMdcMetaAligner::execute2(void) { // // Adapting the old execute // HMdcHit* pHitA; HMdcHit* pHitB; HMdcHit* pHitAB; HShowerHit* pHitS; HTofHit* pHitT; HGeomVector* lSVec = new HGeomVector; HGeomVector* lTVec = new HGeomVector; Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; HLocation local; local.setNIndex(2); local.set(2,sector,modB); HGeomVector projPoint; HGeomVector projMetaPoint; Float_t* projSlopes = new Float_t[2]; Bool_t usedA = kFALSE; Bool_t invalidA = kFALSE; Bool_t SUsed = kFALSE; Bool_t SInvalid = kFALSE; Bool_t TUsed = kFALSE; Bool_t TInvalid = kFALSE; Int_t detnum=-1; Int_t modTof =0; fNEntries++; fIter1->Reset(); fIter1->gotoLocation(local); while ( (pHitB =(HMdcHit*)fIter1->Next()) != 0 ){ fHitsMdc[0]++; //Hits found in MDC B usedA = kFALSE; //reinit control flag invalidA = kFALSE; SUsed = kFALSE; SInvalid = kFALSE; TUsed = kFALSE; TInvalid = kFALSE; if(MA_DEBUG>3) { cout << " ----- SECTOR " << sector << " -----"<< endl; cout << "Module " << modB << ", fHitsMdc " << fHitsMdc[0] << endl; } // Calculation of cross point and slopes in MDC A // The rotation matrix here used is inverse of fRotMat // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // That is: X(A) = fRotMat[0]-1 * (X(B) - fTranslation[0]) // ==> X(A) = fRotMat[0]-1 * X(B) - fRotMat[0]-1 * fTranslation[0] // because we are calculating the projection of MDC B Hits on MDC A // HGeomRotation rotInv = fRotMat[0].inverse(); HGeomVector trasInv = -(rotInv * fTranslation[0]); projPoint = findProjPoint(pHitB,rotInv,trasInv,projSlopes); //Iteration on the second MDC (MDCA) for each hit in the first (MDCB) local.set(2,sector,modA); fIter2->Reset(); fIter2->gotoLocation(local); while ( (pHitA =(HMdcHit*)fIter2->Next()) != 0 ){ //hits found in MDCA, provided there is a Hit in MDCB fHitsMdc[1]++; if(MA_DEBUG>3) cout << "Module " << modA << ", fHitsMdc " << fHitsMdc[1] << endl; if(isInsideWindow(1,pHitA,projPoint,projSlopes)){ if(usedA == kFALSE){ usedA = kTRUE; // MDCA hits found in window (only used if unique) fHitsFoundInWindow[0]++; fHitsFoundAndUnique[0]++; //Real number of matched hits fCount++; //merging all HITS pHitAB = mergeHits(pHitB,pHitA,fRotMat[0], fTranslation[0]); //Now,we enter in the META part! //projectin on SHOWER plane projMetaPoint = findProjMetaPoint(pHitA,pHitB,fRotMat[0], fTranslation[0], fMRotMat[8], fMTranslation[8]); detnum =-1; //SHOWER iteration fIter5->Reset(); while ( (pHitS =(HShowerHit*)fIter5->Next()) != 0 ){ if(isInsideMetaWindow(sector,pHitS,projMetaPoint,lSVec)){ //SHOWER Hit found in window if(SUsed == kFALSE){ SUsed = kTRUE; fHitsFoundInShowerWindow++; fHitsFoundInShowerAndUnique++; fSCount++; detnum =8; //... } else{ //that is, if SUsed == kTRUE if(SInvalid == kFALSE){ fSCount--; SInvalid = kTRUE; fHitsFoundInShowerAndUnique--; fSDiscart++; detnum = -1; } } } } //TOF iteration fIter6->Reset(); while ( (pHitT =(HTofHit*)fIter6->Next()) != 0 ){ modTof = pHitT->getModule(); projMetaPoint = findProjMetaPoint(pHitA,pHitB,fRotMat[0], fTranslation[0], fMRotMat[modTof], fMTranslation[modTof]); if(isInsideMetaWindow(sector,pHitT,projMetaPoint,lTVec)){ //TOF Hit found in window if(TUsed == kFALSE){ TUsed = kTRUE; fHitsFoundInTofWindow[modTof]++; fHitsFoundInTofAndUnique[modTof]++; fTCount++; if(detnum == 8) detnum=detnum+1+modTof; else detnum = modTof; //... } else{ //that is, if TUsed == kTRUE //POR AHORA NO PERMITO DOS HITS EN TOF.CAMBIAR!! if(TInvalid == kFALSE){ fTCount--; TInvalid = kTRUE; fHitsFoundInTofAndUnique[modTof]--; fTDiscart[modTof]++; detnum = -1; } } } } fHits->setMdcHitA(pHitA); fHits->setMdcHitB(pHitB); lSVec->setZ(detnum); lTVec->setZ(detnum); fHits->setNumberOfMDCs(2); fHits->setDetector(detnum); if(detnum<0) fHits->setNumberOfMetaHits(0); else { fHits->setNumberOfMetaHits(1); if(detnum>8) { fHits->setNumberOfMetaHits(2); fHits->setLocalHitA(lSVec); fHits->setLocalHitB(lTVec); } else if(detnum==8) fHits->setLocalHitA(lSVec); else fHits->setLocalHitA(lTVec); } fHits->setDetector(detnum); } // End of if(usedA == kFALSE){ else{ //that is, if usedA == kTRUE //More than one candidate on window!! Discart if(invalidA == kFALSE){ fCount--; invalidA = kTRUE; fDiscart[0]++; fHitsFoundAndUnique[0]--; } } // End of else{ //that is, if usedA == kTRUE } // End of if(isInsideWindow(pHitA,projPoint,projSlopes)){ } // End of while( (pHitA =(HMdcHit*)fIter2->Next()) != 0 ){ if(MA_DEBUG>3) cout << "Logic values: " << (Int_t)usedA << (Int_t)invalidA << (Int_t)SUsed << (Int_t)SInvalid << (Int_t)TUsed << (Int_t)TInvalid << endl; if(usedA == kTRUE && invalidA != kTRUE && ((SUsed == kTRUE && SInvalid != kTRUE) || (TUsed == kTRUE && TInvalid != kTRUE)) ){ if(MA_DEBUG>3) fHits->print(); fAlignMeta->Fill(); if(fCount%100 ==0) cout << "."<< flush; } } // End of while ( (pHitB =(HMdcHit*)fIter1->Next()) != 0 ){ delete lSVec; delete lTVec; } void HMdcMetaAligner::execute1(void) { // // Only one MDC. // HMdcHit* pHitA; HShowerHit* pHitS; HTofHit* pHitT; HGeomVector* lSVec = new HGeomVector; HGeomVector* lTVec = new HGeomVector; Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; HLocation local; local.setNIndex(2); local.set(2,sector,modA); HGeomVector projPoint; HGeomVector projMetaPoint; Bool_t SUsed = kFALSE; Bool_t SInvalid = kFALSE; Bool_t TUsed = kFALSE; Bool_t TInvalid = kFALSE; Int_t detnum=-1; Int_t modTof =0; fIter1->Reset(); fIter1->gotoLocation(local); while ( (pHitA =(HMdcHit*)fIter1->Next()) != 0 ){ fHitsMdc[0]++; fCount++; //Now,we enter in the META part! //projectin on SHOWER plane projMetaPoint = findProjMetaPoint(pHitA,fMRotMat[8], fMTranslation[8]); detnum =-1; //SHOWER iteration fIter5->Reset(); while ( (pHitS =(HShowerHit*)fIter5->Next()) != 0 ){ if(isInsideMetaWindow(sector,pHitS,projMetaPoint,lSVec)){ //SHOWER Hit found in window if(SUsed == kFALSE){ SUsed = kTRUE; fHitsFoundInShowerWindow++; fHitsFoundInShowerAndUnique++; fSCount++; detnum =8; //... } else{ //that is, if SUsed == kTRUE if(SInvalid == kFALSE){ fSCount--; SInvalid = kTRUE; fHitsFoundInShowerAndUnique--; fSDiscart++; detnum = -1; } } } } //TOF iteration fIter6->Reset(); while ( (pHitT =(HTofHit*)fIter6->Next()) != 0 ){ modTof = pHitT->getModule(); projMetaPoint = findProjMetaPoint(pHitA,fMRotMat[modTof], fMTranslation[modTof]); if(isInsideMetaWindow(sector,pHitT,projMetaPoint,lTVec)){ //TOF Hit found in window if(TUsed == kFALSE){ TUsed = kTRUE; fHitsFoundInTofWindow[modTof]++; fHitsFoundInTofAndUnique[modTof]++; fTCount++; if(detnum == 8) detnum=detnum+1+modTof; else detnum = modTof; //... } else{ //that is, if TUsed == kTRUE //POR AHORA NO PERMITO DOS HITS EN TOF.CAMBIAR!! if(TInvalid == kFALSE){ fTCount--; TInvalid = kTRUE; fHitsFoundInTofAndUnique[modTof]--; fTDiscart[modTof]++; detnum = -1; } } } } fHits->setMdcHitA(pHitA); lSVec->setZ(detnum); lTVec->setZ(detnum); if(detnum<0) fHits->setNumberOfMetaHits(0); else { fHits->setNumberOfMetaHits(1); if(detnum>8) { fHits->setNumberOfMetaHits(2); fHits->setLocalHitA(lSVec); fHits->setLocalHitB(lTVec); } else if(detnum==8) fHits->setLocalHitA(lSVec); else fHits->setLocalHitA(lTVec); } fHits->setDetector(detnum); if( ((SUsed == kTRUE && SInvalid != kTRUE) || (TUsed == kTRUE && TInvalid != kTRUE)) ){ if(MA_DEBUG>3) fHits->print(); fAlignMeta->Fill(); // // It is possible to fill also fAlignMetaCut->Fill(); //with a condition on the Hits // if(fCount%100 ==0) cout << "."<< flush; } } // End of while ( (pHitA =(HMdcHit*)fIter1->Next()) != 0 ){ delete lSVec; delete lTVec; } void HMdcMetaAligner::transfEuler(HGeomRotation eulrot,HGeomVector eulvec, HGeomVector oldV, HGeomVector newV){ // // Transformation from one coordinate system to a new one // newV = eulrot * oldV + eulvec // newV = eulrot * oldV + eulvec; } void HMdcMetaAligner::transfEulerInv(HGeomRotation eulrot,HGeomVector eulvec, HGeomVector oldV, HGeomVector newV){ // // Transformation from one coordinate system to a new one // newV = eulrot * (oldV - eulvec) // newV = eulrot * (oldV - eulvec); } HGeomVector HMdcMetaAligner::findProjMetaPoint(HMdcHit* pHitA, HGeomRotation rotM, HGeomVector traM) { // // Find the projection point of a Hit on META plane // Uses a Hit to project on the META plane // defined by rotM and traM // // Two methods are defined. One suposes target at (0,0,0) // (fTar =kTRUE) and the other uses the hit slopes. // HGeomVector point,direct,result; HGeomVector pointInMeta,directInMeta; Float_t third; third = sqrt(1 - pHitA->getXDir()*pHitA->getXDir() - pHitA->getYDir()*pHitA->getYDir()); point.setX(pHitA->getX()); //getting the hit info point.setY(pHitA->getY()); point.setZ(0.); direct.setX(pHitA->getXDir()); direct.setY(pHitA->getYDir()); direct.setZ(third); //transfor to the new system pointInMeta = rotM * point + traM; directInMeta = rotM * direct; // x-pM.getX() y-pM.getY() -pM.getZ() // ------------------- = ------------------- = ------------------- // direcInMeta.getX() directInMeta.getY() directInMeta.getZ() // result.setX(pointInMeta.getX()- pointInMeta.getZ()*directInMeta.getX()/directInMeta.getZ()); result.setY(pointInMeta.getY()- pointInMeta.getZ()*directInMeta.getY()/directInMeta.getZ()); result.setZ(0.); return result; } HGeomVector HMdcMetaAligner::findProjMetaPoint(HMdcHit* pHitA, HMdcHit* pHitB, HGeomRotation rot, HGeomVector tra, HGeomRotation rotM, HGeomVector traM) { // // Find the projection point of a Hit pair on META plane // // Given a relative rotation and translation from MDC A to MDC B // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // X(B) = rot X(A) + tra // this function obtains the MDC A coordinates in // MDC B coordinate system using rot and tra. // Then uses both Hits to project on the META plane // defined by rotM and traM // HGeomVector a,b,anew,ainmeta,binmeta; HGeomVector newvec; a.setX(pHitA->getX()); //getting the hit info a.setY(pHitA->getY()); a.setZ(0.); b.setX(pHitB->getX()); //getting the hit info b.setY(pHitB->getY()); b.setZ(0.); //coordinates of MDC A in MDC B anew = rot * a + tra; ainmeta = rotM * anew + traM; binmeta = rotM * b + traM; //straigth line from both hits and cross point with plane ZMETA=0 // x-aM.getX() y-aM.getY() z-aM.getZ() // ------------------- = ------------------- = ------------------- // bM.getX()-aM.getX() bM.getY()-aM.getY() bM.getZ()-aM.getZ() // // then for z=0: // x=aM.getX()-aM.getZ()*(bM.getX()-aM.getX())/(bM.getZ()-aM.getZ()) // y=aM.getY()-aM.getZ()*(bM.getY()-aM.getY())/(bM.getZ()-aM.getZ()) // Float_t consta = ainmeta.getZ()/(binmeta.getZ()-ainmeta.getZ()); Float_t xproj = ainmeta.getX()-(binmeta.getX()-ainmeta.getX())*consta; Float_t yproj = ainmeta.getY()-(binmeta.getY()-ainmeta.getY())*consta; if(MA_DEBUG>3)cout << "Projected on META: " << xproj << " " << yproj << endl; newvec.setX(xproj); newvec.setY(yproj); newvec.setZ(0.); return newvec; } HGeomVector HMdcMetaAligner::findProjPoint(HMdcHit* pHit, HGeomRotation rot, HGeomVector tra, Float_t* slopes) { // // Find the projection point of a Hit on another MDC // // Given a relative rotation and translation from MDC A to MDC B // (see CONVENTION in HMdcMetaAligner::setGeomAuxPar()) // X(B) = rot X(A) + tra // this function obtains the projection in MDC B of a Hit in MDC A // // When the user wants to obtain the projection in MDC A of a Hit in // MDC B, should provide the inverse transformations // Example: X(A) = rot-1 * X(B) - rot-1 * tra // = newrot * X(B) + newtra // where: // newrot = rot-1 // newtra = - (rot-1 * tra) // The function returns also the slopes in the new coordinate system HGeomVector newvec; Float_t x, y, s0, s1; Float_t xDir, yDir, aux, den; Float_t zsearch, xsearch, ysearch; x = pHit->getX(); //getting the hit info y = pHit->getY(); xDir = pHit->getXDir(); yDir = pHit->getYDir(); //if using hit! aux = sqrt(1 - xDir * xDir - yDir * yDir); if(aux == 0.){ s0=1; s1=1;} //non-sense values else{ s0 = xDir/aux; s1 = yDir/aux; } if(MA_DEBUG>3){ cout << "VALID MDC HIT: " << x << " " << y << " " << s0 << " " << s1 << endl; } zsearch = -(x*rot(6) + y*rot(7) + tra(2)) / (rot(8) + s0*rot(6) + s1*rot(7)); xsearch = x*rot(0) + y*rot(1) + tra(0) + zsearch*(s0*rot(0) + s1*rot(1) + rot(2)); ysearch = x*rot(3) + y*rot(4) + tra(1) + zsearch*(s0*rot(3) + s1*rot(4) + rot(5)); //For getting the histograms in slopes and also for new cuts!! den = s0*rot(6) + s1*rot(7) + rot(8); if (den == 0) { cout << "ERROR in HMdcMetaAligner::findProjPoint()" << endl; return newvec; } slopes[0] = (s0*rot(0) + s1*rot(1) + rot(2)) / den; slopes[1] = (s0*rot(3) + s1*rot(4) + rot(5)) / den; if(MA_DEBUG>3){ cout << "Projected MDC HIT: " << xsearch << " " << ysearch << " " << slopes[0] << " " << slopes[1] << endl; } newvec.setX(xsearch); newvec.setY(ysearch); newvec.setZ(zsearch); return newvec; } Int_t HMdcMetaAligner::isInsideMetaWindow(Int_t sector, HTofHit* pHitT, HGeomVector proj, HGeomVector* vec){ // //Check if the TOF hit is inside a window around proj point //If it is, include in vec //montones de parametros que deberian estar en contenedores estan a mano! Float_t largeCellWidth = 30.0; Float_t smallCellWidth = 20.0; Float_t xTof = pHitT->getXpos(); // /10?? Other units?? Int_t tofCell = (Int_t)pHitT->getCell(); Int_t tofModule = (Int_t)pHitT->getModule(); Float_t yTof; if(tofModule<4) yTof = largeCellWidth * (3.5 - tofCell); else yTof = smallCellWidth * (3.5 - tofCell); if (pHitT->getSector() == sector) { if (tofModule == 0){ fAllResTof0X->Fill(xTof-proj.getX()); fAllResTof0Y->Fill(yTof-proj.getY()); } if (tofModule == 1){ fAllResTof1X->Fill(xTof-proj.getX()); fAllResTof1Y->Fill(yTof-proj.getY()); } if (tofModule == 2){ fAllResTof2X->Fill(xTof-proj.getX()); fAllResTof2Y->Fill(yTof-proj.getY()); } if (tofModule == 3){ fAllResTof3X->Fill(xTof-proj.getX()); fAllResTof3Y->Fill(yTof-proj.getY()); } if (tofModule == 4){ fAllResTof4X->Fill(xTof-proj.getX()); fAllResTof4Y->Fill(yTof-proj.getY()); } if (tofModule == 5){ fAllResTof5X->Fill(xTof-proj.getX()); fAllResTof5Y->Fill(yTof-proj.getY()); } if (tofModule == 6){ fAllResTof6X->Fill(xTof-proj.getX()); fAllResTof6Y->Fill(yTof-proj.getY()); } if (tofModule == 7){ fAllResTof7X->Fill(xTof-proj.getX()); fAllResTof7Y->Fill(yTof-proj.getY()); } if(MA_DEBUG>3) cout << "VALID TOF HIT: " << sector << " " << tofModule << " " << tofCell << " " << xTof << " " << yTof << endl; if((fabs(xTof-proj.getX())3) cout << " inside window" << endl; vec->setX(xTof); vec->setY(yTof); vec->setZ(0); return (tofModule+1); } else return 0; } else return 0; } Int_t HMdcMetaAligner::isInsideMetaWindow(Int_t sector, HShowerHit* pHitS, HGeomVector proj, HGeomVector* vec){ // //Check if the SHOWER hit is inside a window around proj point //If it is, include in vec //SHOWER //getting info from the address Int_t SAddress = pHitS->getAddress(); Int_t SCol = SAddress%100; Int_t SRow = ((SAddress-SCol)/100)%100; Int_t SMod = ((SAddress-SRow-SCol)/10000)%10; //Int_t SSec = ((SAddress-SRow-SCol-SMod)/100000)%10; //NO. SECTOR 0 = SECTOR 6 arriba!!! Int_t SSec = (Int_t)pHitS->getSector(); Float_t hitS[3]; //getting the pad coordinates fShowerGeometry->getPadParam(0)-> getPad(SRow, SCol)->getPadCenter(&hitS[0],&hitS[1]); hitS[2] = 0.; //cm to mm hitS[0] = hitS[0]*10.; hitS[1] = hitS[1]*10.; if(MA_DEBUG>3){ cout << "VALID SHOWER HIT: " << SAddress << " " << SSec << " " << SMod << " " << SRow << " " << SCol << " " << hitS[0] << " " << hitS[1] << endl; //cout << "TESTING SHOWER HIT: " << (Int_t)pHitS->getSector() // << " " << (Int_t)pHitS->getModule() << " " // << (Int_t)pHitS->getRow() << " " // << (Int_t)pHitS->getCol() << endl; } if(SMod == 0 && SSec == sector){ ShowerPlaneProj->Fill(hitS[0],hitS[1]); fAllResShowerX->Fill(hitS[0]-proj.getX()); fAllResShowerY->Fill(hitS[1]-proj.getY()); if((fabs(hitS[0]-proj.getX())3) cout << " inside window" << endl; vec->setVector(hitS); return 8; } else return 0; } else return 0; } Bool_t HMdcMetaAligner::isInsideWindow(Int_t plot, HMdcHit* pHit, HGeomVector point, Float_t* slope){ // //Check if the hit is inside a window around point and slope // //old check. Implement new one PLEASE! //New one base on equiprobability elipse (Beatriz paper) Float_t xlolimit,xuplimit,ylolimit,yuplimit; xlolimit = point.getX() - fXArea; xuplimit = point.getX() + fXArea; ylolimit = point.getY() - fYArea; yuplimit = point.getY() + fYArea; if(plot){ Float_t x, y, s0, s1; Float_t xDir, yDir, aux; x = pHit->getX(); //getting the hit info y = pHit->getY(); xDir = pHit->getXDir(); yDir = pHit->getYDir(); //if using hit! aux = sqrt(1 - xDir * xDir - yDir * yDir); if(aux == 0.){ s0=1; s1=1;} //non-sense values else{ s0 = xDir/aux; s1 = yDir/aux; } fResX->Fill(x - point.getX()); fResY->Fill(y - point.getY()); fResS0->Fill(s0 - slope[0]); fResS1->Fill(s1 - slope[1]); } if(MA_DEBUG>3) cout << "MDC HIT: " << pHit->getX() << " " << pHit->getY(); if( (pHit->getX()>xlolimit) && (pHit->getX()getY()>ylolimit) && (pHit->getY()3) cout << " inside window" << endl; return kTRUE; } else { if(MA_DEBUG>3) cout << " outside window" << endl; return kFALSE; } /* //New check Float_t corr2X = 0.0484; //(corr x-x')^2 from Bea paper Float_t corr2Y = 0.0841; //corr y-y')^2 from Bea paper Float_t sigma2X = 0.04; //(sigma(X))^2 from Bea paper Float_t sigma2Y = 0.01; //(sigma(Y))^2 from Bea paper Float_t sigma2S0 = 0.01; //(sigma(S0))^2 from Bea paper Float_t sigma2S1 = 0.01; //(sigma(S1))^2 from Bea paper Float_t alpha = 3.29; //0.999% Float_t resX2 = (x - point.getX())(x - point.getX()); Float_t resY2 = (y - point.getY())(y - point.getY()); Float_t resS02 = (s0 - slope[0]) * (s0 - slope[0]); Float_t resS12 = (s1 - slope[1]) * (s1 - slope[1]); Float_t partX = (1/(1-corr2X)) * ( (resX2/sigma2X) + (resS02/sigma2S0) - (2*(x - point.getX())*(s0 - slope[0])* corr2X/(sigma2X*sigma2S0)) ); Float_t partY = (1/(1-corr2Y)) * ( (resY2/sigma2Y) + (resS12/sigma2S1) - (2*(y - point.getY())*(s1 - slope[1])* corr2Y/(sigma2Y*sigma2S1)) ); if(partX+partY < alpha) { if(MA_DEBUG>3) cout << " inside window" << endl; return kTRUE; } else { if(MA_DEBUG>3) cout << " outside window" << endl; return kFALSE; } */ } HMdcHit* HMdcMetaAligner::mergeHits(HMdcHit* hitB, HMdcHit* hitA, HGeomRotation rot,HGeomVector tra){ // // Propagates hitA in hitB coordinate system and merges // the information in a new hit in hitB coordinate system // The rot and tra verifies the CONVENTION // (see HMdcMetaAligner::setGeomAuxPar()), that is: // X(B) = rot X(A) + tra // Normally B is reference MDC, which is farther from target //HMdcHit* pHit = new HMdcHit(); //Now, only mean values. Later, propagate the covariance and add the //multiple scattering effect on hitA /* pHit->setX(,); pHit->setY(,); pHit->setXDir(,); pHit->setYDir(,); return pHit; */ //Only for test purposes, return hitB return hitB; } Float_t* HMdcMetaAligner::transformToSlopes(HMdcHit* pHit){ // //Transform hit angular components in slopes // Float_t* slopes = new Float_t[2]; Float_t xDir,yDir; Float_t aux; xDir = pHit->getXDir(); yDir = pHit->getYDir(); aux = sqrt(1 - xDir * xDir - yDir * yDir); if(aux == 0.){ slopes[0]=1; slopes[1]=1;} //non-sense values else{ slopes[0] = xDir/aux; slopes[1] = yDir/aux; } return slopes; } Bool_t HMdcMetaAligner::finalize(void) { // // Statistical information and Minimization procedure // Int_t sector = fLoc[0]; Int_t modA = fLoc[1]; Int_t modB = fLoc[2]; Int_t modC=-1; Int_t modD=-1; if(fNumMods>2) modC = fLoc[3]; if(fNumMods>3) modD = fLoc[4]; //first statistical information ofstream *fout=NULL; if (fNameAscii) fout = new ofstream(fNameAscii.Data(), ios::app); if (*fout){ *fout << endl << "Sector: " << sector << endl; *fout << "Module A: " << modA; if(fNumMods>1) *fout << " Module B: " << modB; if(fNumMods>2) *fout << " Module C: " << modC; if(fNumMods>3) *fout << " Module D: " << modD; *fout << endl << endl; *fout << "Number of events: " << fNEntries << endl; *fout << "Window (mm): " << fXArea <<"," << fYArea << endl; *fout << "Interpret smaller MDC index as farther away from target" << endl << endl; for(Int_t i=0;i1){ for(Int_t i=0;i" << fSCount << " and TOF->" << fTCount << endl; } if(MA_DEBUG>0){ cout << endl << "Sector: " << sector << endl; cout << "Module A: " << modA; if(fNumMods>1) cout << " Module B: " << modB; if(fNumMods>2) cout << " Module C: " << modC; if(fNumMods>3) cout << " Module D: " << modD; cout << endl << endl; cout << "Number of events: " << fNEntries << endl; cout << "Window (mm): " << fXArea <<"," << fYArea << endl; cout << "Interpret smaller MDC index as farther away from target" << endl << endl; for(Int_t i=0;i1){ for(Int_t i=0;i" << fSCount << " and TOF->" << fTCount << endl; } if(fMinOff==kTRUE){ //if no minimization is performed fillHistograms(0); fillMetaHistograms(0); fitHistograms(0); fillHistograms(1,1); fillMetaHistograms(1,1); storeInFile(); fout->close(); delete fout; return kTRUE; } //filling histograms before minimization fillHistograms(0); fillMetaHistograms(0); fitHistograms(0); fillHistograms(1,1); fillMetaHistograms(1,1); if (*fout) *fout << "Valid hits for META alignment after cuts: " << fCountCut << endl << endl; if(MA_DEBUG>0) cout << "Valid hits for META alignment after cuts: " << fCountCut << endl << endl; if (*fout){ *fout << "Transformation before minimization (init values): " << endl; for(Int_t i=0;i<9;i++){ *fout << "Detector " << i << ", "<< fMEuler[i].getX() << ", " << fMEuler[i].getY() << ", " << fMEuler[i].getZ() << ", " << fMTranslation[i].getX() << ", " << fMTranslation[i].getY() << ", " << fMTranslation[i].getZ() << endl; } } HGeomVector OldEuler; HGeomVector OldTranslation; Int_t IterCounter = 0; Float_t IterCri; for(Int_t i=0;i<9;i++){ fDetector = i; do{ OldEuler = fMEuler[i]; OldTranslation = fMTranslation[i]; IterCounter = 0; IterCri = 0; minfit(fFix,fMEuler[i],fMTranslation[i]); if (*fout){ *fout << "Parameters after minimization " << endl; *fout << "Detector" << i << "> " <10) { cout << "WARNING in HMdcMetaAligner :: finalize" << endl; cout << "Sector: sector. " << "Detector: " <fIterCriteria); fillMetaRotMatrix(i,fMEuler[i].getX(),fMEuler[i].getY(),fMEuler[i].getZ()); } //Recalculate all histograms for the new parameters!!! fillMetaHistograms(2); fillMetaHistograms(3,1); storeInFile(); // close ascii file fout->close(); delete fout; return kTRUE; } void HMdcMetaAligner::fillMetaHistograms (Int_t Mindex, Int_t select = 0){ // // Performs the graphical output from obtained parameters // HMdcHit* hitA; HMdcHit* hitB = NULL; HMdcHit* hitC = NULL; HMdcHit* hitD = NULL; HGeomVector* lHit; HGeomVector* lHitB; HGeomVector projMP; Stat_t entries; if(select != 0) entries = fAlignMetaCut->GetEntries(); else entries = fAlignMeta->GetEntries(); for (Int_t i=0;iGetEntry(i); else fAlignMeta->GetEntry(i); if(MA_DEBUG>0) fHits->print(); hitA = (HMdcHit*) fHits->getMdcHitA(); if(fNumMods>1)hitB = (HMdcHit*) fHits->getMdcHitB(); if(fNumMods>2) hitC = (HMdcHit*) fHits->getMdcHitC(); if(fNumMods>3) hitD = (HMdcHit*) fHits->getMdcHitD(); lHit = fHits->getLocalHitA(); //if(lHit->getZ()>8)lHitB = fHits->getLocalHitB(); if(fHits->getDetector()>8)lHitB = fHits->getLocalHitB(); //Algo va mal con los Int_t de HMdcMetaHit en el TTree!! Int_t det = fHits->getDetector(); //Int_t det = lHit->getZ(); if(MA_DEBUG>0) cout << "detector: " << det << endl; //Constructing all possible projections //The histos represent (value of local hit - value of projected hit) //In any case (Meta part) if(fNumMods>1){ if(det<9) projMP = findProjMetaPoint(hitA,hitB,fRotMat[0], fTranslation[0], fMRotMat[det], fMTranslation[det]); else projMP = findProjMetaPoint(hitA,hitB,fRotMat[0], fTranslation[0], fMRotMat[det-9], fMTranslation[det-9]); } else { if(det<9) projMP = findProjMetaPoint(hitA,fMRotMat[det], fMTranslation[det]); else projMP = findProjMetaPoint(hitA,fMRotMat[det-9], fMTranslation[det-9]); } if(det>7) fResShowerX[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det>7) fResShowerY[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==0) fResTof0X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==0) fResTof0Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==1) fResTof1X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==1) fResTof1Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==2) fResTof2X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==2) fResTof2Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==3) fResTof3X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==3) fResTof3Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==4) fResTof4X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==4) fResTof4Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==5) fResTof5X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==5) fResTof5Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==6) fResTof6X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==6) fResTof6Y[Mindex]->Fill(lHit->getY()-projMP.getY()); if(det==7) fResTof7X[Mindex]->Fill(lHit->getX()-projMP.getX()); if(det==7) fResTof7Y[Mindex]->Fill(lHit->getY()-projMP.getY()); } } void HMdcMetaAligner::fillHistograms (Int_t index, Int_t select = 0){ // // Performs the graphical output from obtained parameters // HMdcHit* hitA; HMdcHit* hitB=NULL; HMdcHit* hitC=NULL; HMdcHit* hitD=NULL; HGeomVector projPoint; Float_t* projSlopes = new Float_t[2]; Float_t* origSlopes = new Float_t[2]; HGeomRotation rotaux; HGeomVector transaux; HGeomVector transf[4]; HGeomVector a,b,c,d; Float_t errorx[4]; Float_t errory[4]; Stat_t entries; if(select != 0) entries = fAlignMetaCut->GetEntries(); else entries = fAlignMeta->GetEntries(); for (Int_t i=0;iGetEntry(i); else fAlignMeta->GetEntry(i); hitA = (HMdcHit*) fHits->getMdcHitA(); if(fNumMods>1) hitB = (HMdcHit*) fHits->getMdcHitB(); if(fNumMods>2) hitC = (HMdcHit*) fHits->getMdcHitC(); if(fNumMods>3) hitD = (HMdcHit*) fHits->getMdcHitD(); //Constructing all possible projections //The histos represent (value of local hit - value of projected hit) if(fNumMods==4) { //first projecting on MDC D projPoint = findProjPoint(hitC,fRotMat[0],fTranslation[0],projSlopes); CvsDinDCS_X[index]->Fill(hitD->getX() - projPoint.getX()); CvsDinDCS_Y[index]->Fill(hitD->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitD); CvsDinDCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); CvsDinDCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); projPoint = findProjPoint(hitB,fRotMat[1],fTranslation[1],projSlopes); BvsDinDCS_X[index]->Fill(hitD->getX() - projPoint.getX()); BvsDinDCS_Y[index]->Fill(hitD->getY() - projPoint.getY()); BvsDinDCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsDinDCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); projPoint = findProjPoint(hitA,fRotMat[2],fTranslation[0],projSlopes); AvsDinDCS_X[index]->Fill(hitD->getX() - projPoint.getX()); AvsDinDCS_Y[index]->Fill(hitD->getY() - projPoint.getY()); AvsDinDCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsDinDCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //second projecting on MDC C rotaux = fRotMat[0].inverse(); transaux = -(rotaux * fTranslation[0]); projPoint = findProjPoint(hitD,rotaux,transaux,projSlopes); DvsCinCCS_X[index]->Fill(hitC->getX() - projPoint.getX()); DvsCinCCS_Y[index]->Fill(hitC->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitC); DvsCinCCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); DvsCinCCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[0].inverse())*fRotMat[1]; transaux = (fRotMat[0].inverse())*(fTranslation[1]-fTranslation[0]); projPoint = findProjPoint(hitB,rotaux,transaux,projSlopes); BvsCinCCS_X[index]->Fill(hitC->getX() - projPoint.getX()); BvsCinCCS_Y[index]->Fill(hitC->getY() - projPoint.getY()); BvsCinCCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsCinCCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[0].inverse())*fRotMat[2]; transaux = (fRotMat[0].inverse())*(fTranslation[2]-fTranslation[0]); projPoint = findProjPoint(hitA,rotaux,transaux,projSlopes); AvsCinCCS_X[index]->Fill(hitC->getX() - projPoint.getX()); AvsCinCCS_Y[index]->Fill(hitC->getY() - projPoint.getY()); AvsCinCCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsCinCCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //third projecting on MDC B rotaux = fRotMat[1].inverse(); transaux = -(rotaux * fTranslation[1]); projPoint = findProjPoint(hitD,rotaux,transaux,projSlopes); DvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); DvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitB); DvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); DvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[1].inverse())*fRotMat[0]; transaux = (fRotMat[1].inverse())*(fTranslation[0]-fTranslation[1]); projPoint = findProjPoint(hitC,rotaux,transaux,projSlopes); CvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); CvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); CvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); CvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[1].inverse())*fRotMat[2]; transaux = (fRotMat[1].inverse())*(fTranslation[2]-fTranslation[1]); projPoint = findProjPoint(hitA,rotaux,transaux,projSlopes); AvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); AvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); AvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //last, projecting on MDC A rotaux = fRotMat[2].inverse(); transaux = -(rotaux * fTranslation[2]); projPoint = findProjPoint(hitD,rotaux,transaux,projSlopes); DvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); DvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitA); DvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); DvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[2].inverse())*fRotMat[0]; transaux = (fRotMat[2].inverse())*(fTranslation[0]-fTranslation[2]); projPoint = findProjPoint(hitC,rotaux,transaux,projSlopes); CvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); CvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); CvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); CvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[2].inverse())*fRotMat[1]; transaux = (fRotMat[2].inverse())*(fTranslation[1]-fTranslation[2]); projPoint = findProjPoint(hitB,rotaux,transaux,projSlopes); BvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); BvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); BvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); } if(fNumMods==3){ //first projecting on MDC C projPoint = findProjPoint(hitB,fRotMat[0],fTranslation[0],projSlopes); BvsCinCCS_X[index]->Fill(hitC->getX() - projPoint.getX()); BvsCinCCS_Y[index]->Fill(hitC->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitC); BvsCinCCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsCinCCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); projPoint = findProjPoint(hitA,fRotMat[1],fTranslation[1],projSlopes); AvsCinCCS_X[index]->Fill(hitC->getX() - projPoint.getX()); AvsCinCCS_Y[index]->Fill(hitC->getY() - projPoint.getY()); AvsCinCCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsCinCCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //then projecting on MDC B rotaux = fRotMat[0].inverse(); transaux = -(rotaux * fTranslation[0]); projPoint = findProjPoint(hitC,rotaux,transaux,projSlopes); CvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); CvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitB); CvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); CvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[0].inverse())*fRotMat[1]; transaux = (fRotMat[0].inverse())*(fTranslation[1]-fTranslation[0]); projPoint = findProjPoint(hitA,rotaux,transaux,projSlopes); AvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); AvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); AvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //last, projecting on MDC A rotaux = fRotMat[1].inverse(); transaux = -(rotaux * fTranslation[1]); projPoint = findProjPoint(hitC,rotaux,transaux,projSlopes); CvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); CvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitA); CvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); CvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); rotaux = (fRotMat[1].inverse())*fRotMat[0]; transaux = (fRotMat[1].inverse())*(fTranslation[0]-fTranslation[1]); projPoint = findProjPoint(hitB,rotaux,transaux,projSlopes); BvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); BvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); BvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); } if(fNumMods==2) { //first projecting on MDC B projPoint = findProjPoint(hitA,fRotMat[0],fTranslation[0],projSlopes); AvsBinBCS_X[index]->Fill(hitB->getX() - projPoint.getX()); AvsBinBCS_Y[index]->Fill(hitB->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitB); AvsBinBCS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); AvsBinBCS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); //then projecting on MDC A rotaux = fRotMat[0].inverse(); transaux = -(rotaux * fTranslation[0]); projPoint = findProjPoint(hitB,rotaux,transaux,projSlopes); BvsAinACS_X[index]->Fill(hitA->getX() - projPoint.getX()); BvsAinACS_Y[index]->Fill(hitA->getY() - projPoint.getY()); origSlopes = transformToSlopes(hitA); BvsAinACS_XSlope[index]->Fill(origSlopes[0] - projSlopes[0]); BvsAinACS_YSlope[index]->Fill(origSlopes[1] - projSlopes[1]); } //more histograms if(fNumMods>2) { //First, calculate the straigth line //and all insteresting quantities if(fNumMods>3){ d.setX(hitD->getX()); d.setY(hitD->getY()); d.setZ(0.); } c.setX(hitC->getX()); c.setY(hitC->getY()); c.setZ(0.); b.setX(hitB->getX()); b.setY(hitB->getY()); b.setZ(0.); a.setX(hitA->getX()); a.setY(hitA->getY()); a.setZ(0.); //ERRORS if(fUseUnitErrors==kFALSE){ errorx[0] = (hitA->getErrX()!=0)? hitA->getErrX() : 0.2; errorx[1] = (hitB->getErrX()!=0)? hitB->getErrX() : 0.2; errorx[2] = (hitC->getErrX()!=0)? hitC->getErrX() : 0.2; if(fNumMods>3) errorx[3] = (hitD->getErrX()!=0)? hitD->getErrX():0.2; errory[0] = (hitA->getErrY()!=0)? hitA->getErrY() : 0.1; errory[1] = (hitB->getErrY()!=0)? hitB->getErrY() : 0.1; errory[2] = (hitC->getErrY()!=0)? hitC->getErrY() : 0.1; if(fNumMods>3) errory[3] = (hitD->getErrY()!=0)? hitD->getErrY():0.1; } else for(Int_t i=0;i3){ for(Int_t i=0;i3) cout << "Xwt=" << Xwt << " Xss=" << Xss << " Xsx=" << Xsx << " Xsy=" << Xsy << "Ywt=" << Ywt << " Yss=" << Yss << " Ysx=" << Ysx << " Ysy=" << Ysy << endl; } Xsxoss = Xsx/Xss; Ysxoss = Ysx/Yss; if(MA_DEBUG>3) cout << "Xsxoss=" << Xsxoss << " Ysxoss=" << Ysxoss << endl; for(Int_t i=0;i3) cout << "Xt=" << Xt << " Xst2=" << Xst2 << " b1 (partial)=" << b1 << "Yt=" << Yt << " Yst2=" << Yst2 << " b2 (partial)=" << b2 << endl; } b1 /= Xst2; a1 = (Xsy-(Xsx*b1))/Xss; b2 /= Yst2; a2 = (Ysy-(Ysx*b2))/Yss; if(MA_DEBUG>3) cout << "b1=" << b1 << " a1=" << a1 << "b2=" << b2 << " a2=" << a2 << endl; siga1 = sqrt((1.0+Xsx*Xsx/(Xss*Xst2))/Xss); sigb1 = sqrt(1.0/Xst2); siga2 = sqrt((1.0+Ysx*Ysx/(Yss*Yst2))/Yss); sigb2 = sqrt(1.0/Yst2); if(MA_DEBUG>3) cout << "sigb1=" << sigb1 << " siga1=" << siga1 << "sigb2=" << sigb2 << " siga2=" << siga2 << endl; //Aqui falta evaluar la calidad del ajuste o bien encontrar //cuales son los errores que se esperan en los datos para un buen ajuste for(Int_t i=0;i3) cout << "XChi2=" << XChi2 << " YChi2=" << YChi2 << endl; } XChi2Hist[index]->Fill(XChi2); YChi2Hist[index]->Fill(YChi2); TotalChi2[index]->Fill(XChi2+YChi2); // Also can be defined by a vector V and a point P: // V=(b1,b2,1) P=(a1,a2,0) // Let us calculate the square distance from the straigth // line to the points (fit residuals) Float_t part1=0.0,part2=0.0,part3=0.0,sdistance=0.0,sdist=0.0; for(Int_t i=0;i3) cout << "Square Distance point " << i << " - line: " << sdistance << endl; if(i==0)SqrDistToA[index]->Fill(sdist); if(i==1)SqrDistToB[index]->Fill(sdist); if(i==2)SqrDistToC[index]->Fill(sdist); if(i==3)SqrDistToD[index]->Fill(sdist); } SqrDist[index]->Fill(sdistance); } else{ //2 MDCs } } // end of for (Int_t i=0; i2) modC = fLoc[3]; if(fNumMods>3) modD = fLoc[4]; fAlignMeta->Write(); fAlignMetaCut->Write(); // const char* oldDirName = gDirectory->GetPath(); static Char_t newDirName[255]; if(fNumMods == 1) sprintf(newDirName,"%s%s%i%s%i","MetaAligner_", "S_",sector,"_ModA_",modA); if(fNumMods == 2) sprintf(newDirName,"%s%s%i%s%i%s%i","MetaAligner_", "S_",sector,"_ModA_",modA,"_ModB_",modB); if(fNumMods == 3) sprintf(newDirName,"%s%s%i%s%i%s%i%s%i","MetaAligner_", "S_",sector,"_ModA_",modA,"_ModB_",modB, "_ModC_",modC); if(fNumMods == 4) sprintf(newDirName,"%s%s%i%s%i%s%i%s%i%s%i", "MetaAligner_","S_",sector,"_ModA_",modA, "_ModB_",modB,"_ModC_",modC,"_ModD_",modD); fOutRoot->cd(newDirName); //Enter in the file the HMdcGeomRotation HMdcGeomVector if(fNumMods>1){ for(Int_t i=0;icd(); fRecCount--; if(!fRecCount) { fOutRoot->Write(); fOutRoot->Close(); } } void HMdcMetaAligner::fillMetaRotMatrix(Int_t ind, Float_t prim, Float_t segu, Float_t terc){ // // Fill a matrix using the Euler angles of the relative transformation // /*OLD fRotMat[0][0]=(cos(prim) * cos(segu) * cos(terc)) - (sin(prim) * sin(terc)); fRotMat[1][0]=( - cos(prim) * cos(segu) * sin(terc)) - (sin(prim) * cos(terc)); fRotMat[2][0]=(cos(prim) * sin(segu)); fRotMat[0][1]=(sin(prim) * cos(segu) * cos(terc)) + (cos(prim) * sin(terc)); fRotMat[1][1]=( - sin(prim) * cos(segu) * sin(terc)) + (cos(prim) * cos(terc)); fRotMat[2][1]=(sin(prim) * sin(segu)); fRotMat[0][2]=( - sin(segu) * cos(terc)); fRotMat[1][2]=(sin(segu) * sin(terc)); fRotMat[2][2]=(cos(segu)); */ //ATT! Correspondencia entre HGeomRotation y la antigua transf. arriba // // rot(0) rot(1) rot(2) fRotMat[0][0] fRotMat[1][0] fRotMat[2][0] // rot(3) rot(4) rot(5) = fRotMat[0][1] fRotMat[1][1] fRotMat[2][1] // rot(6) rot(7) rot(8) fRotMat[0][2] fRotMat[1][2] fRotMat[2][2] const double rad2deg = 57.29577951308; fMRotMat[ind].setEulerAngles(prim*rad2deg,segu*rad2deg,terc*rad2deg); } void HMdcMetaAligner::fillMetaTranslation(Int_t ind,Float_t x, Float_t y, Float_t z){ // // Translation from MDC A to MDC B // fMTranslation[ind].setX(x); fMTranslation[ind].setY(y); fMTranslation[ind].setZ(z); } void HMdcMetaAligner::fillRotMatrix(Int_t ind, Float_t prim, Float_t segu, Float_t terc){ // // Fill a matrix using the Euler angles of the relative transformation // /*OLD fRotMat[0][0]=(cos(prim) * cos(segu) * cos(terc)) - (sin(prim) * sin(terc)); fRotMat[1][0]=( - cos(prim) * cos(segu) * sin(terc)) - (sin(prim) * cos(terc)); fRotMat[2][0]=(cos(prim) * sin(segu)); fRotMat[0][1]=(sin(prim) * cos(segu) * cos(terc)) + (cos(prim) * sin(terc)); fRotMat[1][1]=( - sin(prim) * cos(segu) * sin(terc)) + (cos(prim) * cos(terc)); fRotMat[2][1]=(sin(prim) * sin(segu)); fRotMat[0][2]=( - sin(segu) * cos(terc)); fRotMat[1][2]=(sin(segu) * sin(terc)); fRotMat[2][2]=(cos(segu)); */ //ATT! Correspondencia entre HGeomRotation y la antigua transf. arriba // // rot(0) rot(1) rot(2) fRotMat[0][0] fRotMat[1][0] fRotMat[2][0] // rot(3) rot(4) rot(5) = fRotMat[0][1] fRotMat[1][1] fRotMat[2][1] // rot(6) rot(7) rot(8) fRotMat[0][2] fRotMat[1][2] fRotMat[2][2] const double rad2deg = 57.29577951308; fRotMat[ind].setEulerAngles(prim*rad2deg,segu*rad2deg,terc*rad2deg); } void HMdcMetaAligner::fillTranslation(Int_t ind,Float_t x, Float_t y, Float_t z){ // // Translation from MDC A to MDC B // fTranslation[ind].setX(x); fTranslation[ind].setY(y); fTranslation[ind].setZ(z); } void HMdcMetaAligner::setMetaEulerAngles(Int_t ind,Float_t f, Float_t s, Float_t t){ // // Euler angles in transformation MDC B - META // fMEuler[ind].setX(f); fMEuler[ind].setY(s); fMEuler[ind].setZ(t); } void HMdcMetaAligner::setEulerAngles(Int_t ind,Float_t f, Float_t s, Float_t t){ // // Euler angles in transformation from MDC A to MDC B // fEuler[ind].setX(f); fEuler[ind].setY(s); fEuler[ind].setZ(t); } void HMdcMetaAligner::setSearchLimits(Float_t x, Float_t y){ // // Limits of the square defined in the search procedure // fXArea = x; fYArea = y; } void HMdcMetaAligner::setIterCriteria(Float_t cri){ // // Set the criteria for iteration in the minimization (see finalize()) // fIterCriteria = cri; } void HMdcMetaAligner::setWeights(Float_t w0,Float_t w1,Float_t w2,Float_t w3){ // // Set the weights in the fcn() // fWeights[0]= w0; fWeights[1]= w1; fWeights[2]= w2; fWeights[3]= w3; } void HMdcMetaAligner::setSteps(Float_t s0,Float_t s1,Float_t s2, Float_t s3, Float_t s4, Float_t s5){ // // Set the steps in the minimization // fSteps[0]= s0; fSteps[1]= s1; fSteps[2]= s2; fSteps[3]= s3; fSteps[4]= s4; fSteps[5]= s5; if(MA_DEBUG>3) cout << "Steps in the minimization adjusted to " << s0 << ", " << s1 << ", " << s2 << ", " << s3 << ", " << s4 << ", " << s5 << endl; } void HMdcMetaAligner::setLimits(Float_t l0,Float_t l1,Float_t l2, Float_t l3, Float_t l4, Float_t l5){ // // Set the criteria for iteration in the minimization (see finalize()) // fLimits[0]= l0; fLimits[1]= l1; fLimits[2]= l2; fLimits[3]= l3; fLimits[4]= l4; fLimits[5]= l5; if(MA_DEBUG>3) cout << "Limits in the minimization adjusted to " << l0 << ", " << l1 << ", " << l2 << ", " << l3 << ", " << l4 << ", " << l5 << endl; } void HMdcMetaAligner::minfit(Int_t fix, HGeomVector fE, HGeomVector fT){ // // Minuit menu // Put first argument to 1 to fix angular values // Put first argument to 2 to fix translation values // Other arguments are init values for the parameters! // Double_t args[8]; Int_t err = 0; Float_t* limitL; Float_t* limitH; limitL = new Float_t[6]; limitH = new Float_t[6]; //This depends on MDCA and MDCB //DANGEROUS! Can read memory out of scope Double_t start_val[]={fE.getX(),fE.getY(),fE.getZ(), fT.getX(),fT.getY(),fT.getZ()}; //seting limits for(int i=0;i<6;i++){ if(fLimits[i]==0){ limitL[i]=0; limitH[i]=0; } else { limitL[i]=start_val[i]-fLimits[i]; limitH[i]=start_val[i]+fLimits[i]; cout << " LIMITS IN THE MINIMIZATION " << endl; cout << " (from 0 to 5) Par " << i << " between " << limitL[i] << " and " << limitH[i] << endl; } } TMinuit *minuit=new TMinuit(6); //setting the object to minimize minuit->SetFCN(fcnMA); minuit->SetObjectFit(this); if(MA_DEBUG>0){ cout << "HMdcMetaAligner::minfit()" <mnparm(i,pname,start_val[i],fSteps[i],limitL[i],limitH[i],err); } //FIXING if(fix == 1) { //fixing angles minuit->FixParameter(0); minuit->FixParameter(1); minuit->FixParameter(2); } if(fix == 2) { //fixing translations minuit->FixParameter(3); minuit->FixParameter(4); minuit->FixParameter(5); } //FIXING if(fUseSlopeCut) //when cut in slopes, fixing the Z translation minuit->FixParameter(5); //To remove all printout if(MA_DEBUG<3){ minuit->SetPrintLevel(-1); } //args is the array of the numerical arguments of the command //the third field is the number of arguments esp3cified //For MIGRAD arguments are maxcalls and tolerance, both optionals minuit->mnexcm("MIGRAD",args,0,err); //minuit->mnexcm("SIMPLEX",args,0,err); Double_t parresult[6]; for(Int_t i=0;i<6;i++){ minuit->GetParameter(i,parresult[i],fError[i]); } fMEuler[fDetector].setX(parresult[0]); fMEuler[fDetector].setY(parresult[1]); fMEuler[fDetector].setZ(parresult[2]); fMTranslation[fDetector].setX(parresult[3]); fMTranslation[fDetector].setY(parresult[4]); fMTranslation[fDetector].setZ(parresult[5]); if (err != 0) cout << endl <<"MINIMIZATION EXITED WITH ERROR " << err << endl; } void fcnMA(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){ // // This function contains the functional to minimize // Double_t chisq = 0.; HGeomRotation rotM; HGeomRotation rot; HGeomVector tra; HGeomVector a,b,anew,ainmeta,binmeta; HGeomVector point,direct; HGeomVector pointInMeta,directInMeta; Double_t prim, segu, terc; Float_t xproj,yproj,consta; HGeomVector* vec; Int_t ldet; Float_t third; prim = par[0]; segu = par[1]; terc = par[2]; HGeomVector traM(par[3],par[4],par[5]); if (MA_DEBUG>2){ cout << "HMdcMetaAligner::fcnalign() Parameters: " << par[0] << "," << par[1] << "," << par[2] << "," << par[3] << "," << par[4] << "," << par[5] << endl; } rotM.setEulerAngles(prim*180./TMath::Pi(), segu*180./TMath::Pi(), terc*180./TMath::Pi()); /* 0 rotmat[0][0]=(cos(prim)*cos(segu)*cos(terc))-(sin(prim)*sin(terc)); 1 rotmat[1][0]=(-cos(prim)*cos(segu)*sin(terc))-(sin(prim)*cos(terc)); 2 rotmat[2][0]=(cos(prim) * sin(segu)); 3 rotmat[0][1]=(sin(prim)*cos(segu)*cos(terc))+(cos(prim)*sin(terc)); 4 rotmat[1][1]=(-sin(prim)*cos(segu)*sin(terc))+(cos(prim)*cos(terc)); 5 rotmat[2][1]=(sin(prim) * sin(segu)); 6 rotmat[0][2]=( - sin(segu) * cos(terc)); 7 rotmat[1][2]=(sin(segu) * sin(terc)); 8 rotmat[2][2]=(cos(segu)); */ Int_t numMods; TTree* pData; HMdcMetaAligner* pAlign; HMdcMetaHit* theHits; pAlign = (HMdcMetaAligner*)(gMinuit->GetObjectFit()); pData = pAlign->getTree(); numMods = pAlign->getNumMods(); if(numMods>1){ rot = pAlign->getRotMatrix(); tra = pAlign->getTranslation(); } Int_t detector = pAlign->getDetector(); // theHits = pAlign->getHits(); Stat_t entries = pData->GetEntries(); Float_t* weights; weights = new Float_t[4]; weights = pAlign->getWeights(); HMdcHit* hitA; HMdcHit* hitB=NULL; //loop on hits for (Int_t i=0;iGetEntry(i); theHits = pAlign->getHits(); //data = pData->GetArgs(); //filling from the Tree vec = (HGeomVector*) theHits->getLocalHitA(); //ldet = (Int_t) vec->getZ(); //CUando arregle el problema con el tree de HMdcMetaHit, cambiar ldet = (Int_t) theHits->getDetector(); if(ldet==detector||(ldet-9)==detector) { hitA = (HMdcHit*) theHits->getMdcHitA(); if(numMods>1) hitB = (HMdcHit*) theHits->getMdcHitB(); if(ldet>8) vec = (HGeomVector*) theHits->getLocalHitB(); if(numMods>1){ //filling from the ntuple a.setX(hitA->getX()); //getting the hit info a.setY(hitA->getY()); a.setZ(0.); b.setX(hitB->getX()); //getting the hit info b.setY(hitB->getY()); b.setZ(0.); anew = rot * a + tra; ainmeta = rotM * anew + traM; binmeta = rotM * b + traM; //straigth line from both hits and cross point with plane ZMETA=0 // x-aM.getX() y-aM.getY() z-aM.getZ() // ------------------- = ------------------- = ------------------- // bM.getX()-aM.getX() bM.getY()-aM.getY() bM.getZ()-aM.getZ() // // then for z=0: // x=aM.getX()-aM.getZ()*(bM.getX()-aM.getX())/(bM.getZ()-aM.getZ()) // y=aM.getY()-aM.getZ()*(bM.getY()-aM.getY())/(bM.getZ()-aM.getZ()) // consta = ainmeta.getZ()/(binmeta.getZ()-ainmeta.getZ()); xproj = ainmeta.getX()-(binmeta.getX()-ainmeta.getX())*consta; yproj = ainmeta.getY()-(binmeta.getY()-ainmeta.getY())*consta; //Some extra debbuging (always needed ;-) ) if(MA_DEBUG >3){ cout << " ++++++++ VALUES IN fcnMA() +++++++++ " << endl; cout << "HITA: " ; a.print(); cout << "HITB: " ; b.print(); cout << "HITA IN MDC B CS: "; anew.print(); cout << "HITA IN META CS: "; ainmeta.print(); cout << "HITB IN META CS: "; binmeta.print(); cout << "Dif X: " << xproj-(vec->getX()) << " Dif Y: " << yproj-(vec->getY())<< endl; } } else{ //1MDC third = sqrt(1 - hitA->getXDir()*hitA->getXDir() - hitA->getYDir()*hitA->getYDir()); point.setX(hitA->getX()); //getting the hit info point.setY(hitA->getY()); point.setZ(0.); direct.setX(hitA->getXDir()); direct.setY(hitA->getYDir()); direct.setZ(third); //transfor to the new system pointInMeta = rotM * point + traM; directInMeta = rotM * direct; // x-pM.getX() y-pM.getY() -pM.getZ() // ------------------- = ------------------- = ------------------- // direcInMeta.getX() directInMeta.getY() directInMeta.getZ() // xproj = pointInMeta.getX()- pointInMeta.getZ()*directInMeta.getX()/directInMeta.getZ(); yproj = pointInMeta.getY()- pointInMeta.getZ()*directInMeta.getY()/directInMeta.getZ(); } chisq += (xproj-(vec->getX()))*(xproj-(vec->getX()))/weights[0] + (yproj-(vec->getY()))*(yproj-(vec->getY()))/weights[1]; } //end of if(... } //end of for (Int_t i=0;... f = chisq; if(MA_DEBUG>2){ cout << "chisqr= " << chisq << " out of " << entries << " combinations."<< endl; } }