//*--- Author : V.Pechenov //*--- Modified: 26.06.07 V.Pechenov //*--- $Id: hmdcmille.cc,v 1.61 2007/07/09 12:38:36 pechenov Exp $ using namespace std; #include "hmdcmille.h" #include "Mille.h" #include "hades.h" #include "hruntimedb.h" #include "hcategory.h" #include "hiterator.h" #include "hmdcsizescells.h" #include "hmdcclussim.h" #include "hmdctrackdset.h" #include "hmdctrackparam.h" #include "hmdcwiredata.h" #include "hgeantkine.h" #include "hgeantmdc.h" #include "hmdcgetcontainers.h" #include "hmdclayergeompar.h" #include "hmdclayercorrpar.h" #include "hmdcgeompar.h" #include "hspecgeompar.h" #include "hmdcgeomstruct.h" #include "hparasciifileio.h" #include "hmdcdetector.h" #include "hmdclistcells.h" #include #include #include //_HADES_CLASS_DESCRIPTION ////////////////////////////////////////////////////////////////////////////// // // HMdcMille // // Interface class for Millipede alignment. // ////////////////////////////////////////////////////////////////////////////// ClassImp(HMdcMille) HMdcMille::HMdcMille() : HMdc12Fit("mdc.Mille","mdc.Mille") { setDef(); } HMdcMille::HMdcMille(const Char_t* milleOutFName,const Char_t* geomPFName, Int_t bmTmID,const Char_t* sumOfShFName) : HMdc12Fit("mdc.Mille","mdc.Mille") { setDef(); if(bmTmID<0 || bmTmID>9) printErrorAndExit("HMdcMille","beamTimeId must be >= 0 and <=9 !!!"); beamTimeId = bmTmID; //! if(milleOutFName) mille = new Mille(milleOutFName); //! else printErrorAndExit("HMdcMille", //! "Mille output file name is not specifaed !!!"); milleFileName = milleOutFName; openNewBinaryFile(); if(sumOfShFName != 0) sumOfShiftsFName = sumOfShFName; if(geomPFName != 0) geomParFileName = geomPFName; readSumOfShiftsFile(); readPedeResFile(pedeResFName + ".res"); if(!loadGeometryPar()) printErrorAndExit("HMdcMille","Can't load mdc geometry!!! Stop!"); if(!makeShifts()) printErrorAndExit("HMdcMille","Can't change mdc geometry!"); char buf[50]; sprintf(buf,".st%02i",iteration); stepIter = buf; } void HMdcMille::setDef(void) { HMdcTrackDSet::setAnotherFit(this); nParMax = 9; iteration = 0; mille = 0; shiftType = 1; doCopyGeomFile = kFALSE; doCopySumShFile = kTRUE; doCopyResFile = kTRUE; doCopyLogFile = kFALSE; doCopyHisFile = kFALSE; createPedeInFile = kFALSE; // kTRUE - create par.file isGeomChanged = kFALSE; isGeomFileExist = kFALSE; nTracks = 0; nWiresTot = 0; nLayersTot = 0; parSteps[0] = 0.01; // shift along X parSteps[1] = 0.01; // shift along Y parSteps[2] = 0.01; // shift along Z parSteps[3] = 0.003; // rotation around X (in degree) parSteps[4] = 0.003; // rotation around Y (in degree) parSteps[5] = 0.003; // rotation around Z (in degree) parSteps[6] = 0.004; // cell thickness cellThicknFree = kFALSE; nHists = 0; fixFitTOffset = kTRUE; useMdcShParOnly = kFALSE; doLayP2Align = kTRUE; nWiresMin = 15; thetaCut = -1.; pGetCont = HMdcGetContainers::getObject(); pMdcDet = pGetCont->getMdcDetector(); nBinaryFiles = 0; milleFileSize = 0; crPedeTaskFile = kFALSE; printDebugFlag = kFALSE; doHists = kTRUE; nClustersCut = 5; //??? filteringFlag = 0; for(Int_t s=0;s<6;s++) { useSector[s] = kTRUE; nWiresCut[s] = 320; } // Clean arrais: for(Int_t im=0;im<24;im++) pSCMod[im] = NULL; for(Int_t il=0;il<144;il++) pSCLayer[il] = NULL; for(Int_t il=0;il<144;il++) nLayerParts[il] = 0; for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) { histInd[s] = -1; mdcMods[s][m] = 0; shiftSysMdc[s][m] = -1; for(Int_t p=0;pdeleteCont(); } Bool_t HMdcMille::init(void) { return HMdc12Fit::init(); } Bool_t HMdcMille::reinit(void) { if( !HMdc12Fit::reinit() ) return kFALSE; HMdcSizesCells* pSizesCells=HMdcSizesCells::getExObject(); if(pSizesCells == 0) return kFALSE; setShiftTransformation(); for(Int_t sec=0;sec<6;sec++) { nMods[sec] = 0; HMdcSizesCellsSec& pSizesCellsSec = (*pSizesCells)[sec]; if(&pSizesCellsSec == 0) continue; if(doHists) createHists(sec); for(module=0;module<4;module++) { HMdcSizesCellsMod& pSizesCellsMod = pSizesCellsSec[module]; if(&pSizesCellsMod == 0) continue; nMods[sec]++; setIMod(sec,module); // calculate iMod pSCMod[iMod] = &pSizesCellsMod; for(layer=0;layer<6;layer++) { HMdcSizesCellsLayer& rSizesCellsLay = pSizesCellsMod[layer]; setILay(layer); // calculate iLay if(&rSizesCellsLay == 0) continue; pSCLayer[iLay] = &rSizesCellsLay; nLayerParts[iLay] = rSizesCellsLay.getLayerNParts(); if(shiftType == 0) setTrShiftInSecSys(); else if(shiftType == 1) setTrShiftInModSys(); } } } return kTRUE; } void HMdcMille::createHists(Int_t sec) { if(fHistograms != 0 && nHists>0) return; if(fHistograms == 0) fHistograms = new TList(); char name[100]; char title[100]; histInd[sec] = nHists; sprintf(name,"pThVsPhS%i",sec+1); sprintf(title,"Sec.%i #Theta vs #phi",sec+1); fHistograms->Add(new TH2F(name,title,60,60.,120.,80,10.,90.)); nHists++; sprintf(name,"hChi2S%i",sec+1); sprintf(title,"Sec.%i Chi^{2}/NDF",sec+1); fHistograms->Add(new TH1F(name,title,200,0.,20.)); nHists++; sprintf(name,"hNWiresS%i",sec+1); sprintf(title,"Sec.%i Num.wires per track",sec+1); Int_t nBins = clusFitAlg==4 ? 100 : 50; fHistograms->Add(new TH1F(name,title,nBins,0.,nBins)); nHists++; nBins = clusFitAlg==4 ? 49 : 25; sprintf(name,"hNLayersS%i",sec+1); sprintf(title,"Sec.%i Num.layers per track",sec+1); fHistograms->Add(new TH1F(name,title,nBins,0.,nBins)); nHists++; } void HMdcMille::setShiftTransformation(void) { Double_t steps[6]={0.,0.,0.,0.,0.,0.}; for(Int_t parn=0;parn<6;parn++) { derNorm[parn] = 1./parSteps[parn]; steps[parn] = parSteps[parn]; HMdcSizesCells::setTransform(steps,posShifts[parn]); steps[parn] = -parSteps[parn]; HMdcSizesCells::setTransform(steps,negShifts[parn]); steps[parn] = 0.; } // For layer thickness align. parameter: for(Int_t lay=0;lay<6;lay++) { steps[2] = parSteps[6]*(lay-3+0.5); derLThNorm[lay] = 1./steps[2]; HMdcSizesCells::setTransform(steps,posLThSh[lay]); steps[2] = -steps[2]; HMdcSizesCells::setTransform(steps,negLThSh[lay]); steps[2] = 0.; } } void HMdcMille::setTrShiftInModSys(void) { // mdcSys - mdc<->sec. transfomation // laySysMod - rotated_layer<->mdc transfomation const HGeomTransform *mdcSys = pSCMod[iMod]->getSecTrans(); const HGeomTransform *layMod = pSCLayer[iLay]->getRotLayP1SysRMod(); const HGeomTransform *layP2Mod = NULL; if(nLayerParts[iLay] == 2) layP2Mod = pSCLayer[iLay]->getRotLayP2SysRMod(); for(Int_t parn=0;parn<6;parn++) { calcShiftInModSys(mdcSys,layMod,posShifts[parn],layPosSh[iLay][parn]); calcShiftInModSys(mdcSys,layMod,negShifts[parn],layNegSh[iLay][parn]); if(layP2Mod == NULL) continue; calcShiftInModSys(mdcSys,layP2Mod,posShifts[parn],layPosSh[iLay][parn+20]); calcShiftInModSys(mdcSys,layP2Mod,negShifts[parn],layNegSh[iLay][parn+20]); } // Cell thickness: calcShiftInModSys(mdcSys,layMod,posLThSh[layer],layPosSh[iLay][6]); calcShiftInModSys(mdcSys,layMod,negLThSh[layer],layNegSh[iLay][6]); if(layP2Mod != NULL) { calcShiftInModSys(mdcSys,layP2Mod,posLThSh[layer],layPosSh[iLay][6+20]); calcShiftInModSys(mdcSys,layP2Mod,negLThSh[layer],layNegSh[iLay][6+20]); } } void HMdcMille::setTrShiftInSecSys(void) { // laySysSec - rotated_layer<->sec transfomation const HGeomTransform* laySec = pSCLayer[iLay]->getRotLayP1SysRSec(); const HGeomTransform* layP2Sec = NULL; if(nLayerParts[iLay] == 2) layP2Sec = pSCLayer[iLay]->getRotLayP2SysRSec(); for(Int_t parn=0;parn<6;parn++) { calcShiftInSecSys(laySec,posShifts[parn],layPosSh[iLay][parn]); calcShiftInSecSys(laySec,negShifts[parn],layNegSh[iLay][parn]); if(layP2Sec == NULL) continue; // Layer second part: calcShiftInSecSys(layP2Sec,posShifts[parn],layPosSh[iLay][parn+20]); calcShiftInSecSys(layP2Sec,negShifts[parn],layNegSh[iLay][parn+20]); } // Cell thickness: calcShiftInSecSys(laySec,posLThSh[layer],layPosSh[iLay][6]); calcShiftInSecSys(laySec,negLThSh[layer],layNegSh[iLay][6]); if(layP2Sec == NULL) return; // Layer second part: calcShiftInSecSys(layP2Sec,posLThSh[layer],layPosSh[iLay][6+20]); calcShiftInSecSys(layP2Sec,negLThSh[layer],layNegSh[iLay][6+20]); } void HMdcMille::calcShiftInModSys(const HGeomTransform* mdcSys, const HGeomTransform* laySysMod,const HGeomTransform& shift, HGeomTransform& laySh) { // mdcSys - mdc<->sec. transfomation // laySysMod - rotated_layer<->mdc transfomation // shift - transfomation for shift of one global parameter // laySh - new (shifted) rotated_layer<->sec transfomation HGeomTransform shiftedMdcSys(shift); shiftedMdcSys.transFrom(*mdcSys); // Do shift of mdc<->sec. transfomation laySh = *laySysMod; laySh.transFrom(shiftedMdcSys); // Calc. new lay<->sec. transfomation } void HMdcMille::calcShiftInSecSys( const HGeomTransform* laySysSec, const HGeomTransform& shift,HGeomTransform& laySh) { // laySysSec - rotated_layer<->sec transfomation // shift - transfomation for shift of one global parameter // laySh - new (shifted) rotated_layer<->sec transfomation laySh = *laySysSec; laySh.transFrom(shift); // Do shift of layer<->sec. transfomation } Bool_t HMdcMille::finalize(void) { if(mille) delete mille; mille = 0; creatSumOfShiftsFile(); if(!writeParAsciiFile()) printErrorAndExit("finalize", "Can't write geometry parameters ASCII file!!!"); if(createPedeInFile) creatPedeInParamFile(); if(crPedeTaskFile) writePedeTaskFile(); return kTRUE; } Int_t HMdcMille::execute(void) { if(HMdcTrackDSet::fPrint()) { fitpar.setPrintFlag(gHades->getEventCounter()); if(fitpar.getPrintFlag()) printf( "============================== Event %i =============================\n", gHades->getEventCounter()); } if(clusFitAlg==3) { if(filteringFlag>0) { if( !eventFilter() ) return kSkipEvent; if( filteringFlag==2 ) return 0; } fitAlgorithmForMille(); } else if(clusFitAlg==4) fitAlgForMilleCosmic(); return 0; } Bool_t HMdcMille::eventFilter(void) { // For beam data only! HMdcClus* pClst; UInt_t nClusters[6] = {0,0,0,0,0,0}; iterClus->Reset(); while((pClst=(HMdcClus*)iterClus->Next())) if(pClst->getIndexParent()<0) { UChar_t sec = pClst->getSec(); nClusters[sec]++; } HMdcEvntListCells* pEvLCells = HMdcEvntListCells::getExObject(); Bool_t isGoodSectors = kFALSE; for(Int_t s=0;s<6;s++) { useSector[s] = kFALSE; if(nMods[s] <= 1) continue; if(nClusters[s]>0 && nClusters[s]<=nClustersCut) { if(pEvLCells != NULL) { Int_t nCells = (*pEvLCells)[s].getNCells(); Int_t nWiresC = (nWiresMin*nMods[s])/4; if(nCells>=nWiresC && nCells<=nWiresCut[s]) { useSector[s] = kTRUE; isGoodSectors = kTRUE; } } else { useSector[s] = kTRUE; isGoodSectors = kTRUE; } } } return isGoodSectors; } void HMdcMille::fitAlgorithmForMille(void) { // Sector fit // Magnet off + combined clusters in sector HMdcClus* fClst1; iterClus->Reset(); while((fClst1=(HMdcClus*)iterClus->Next())) if(fClst1->getIndexParent()<0) { UChar_t sec = fClst1->getSec(); if(!useSector[sec]) continue; //--- Angle theta cut: //!!! if(thetaCut>0. && fClst1->getTheta()*TMath::RadToDeg()getTheta()*TMath::RadToDeg(); // Double_t phi = fClst1->getPhi()*TMath::RadToDeg(); // Int_t sec = fClst1->getSec(); /* if(sec==1) { if(theta<35.) { if(phi>=108. && theta<35.) continue; if(phi>=102. && theta<25.) continue; if(phi>=90. && theta<20.) continue; } } else if(sec==2) { if(theta<30.) { if(phi>=96. && theta<30.) continue; if(phi>=78. && theta<25.) continue; } } else if(sec==3) { if(phi>=96. && theta<25.) continue; } else if(sec==4) { if(phi>=108. && theta<25.) continue; } */ fittersArr[0].reset(); Int_t first,last; fClst1->getIndexRegChilds(first,last); HMdcClus* fClst2 = (first<0) ? 0:(HMdcClus*)fClusCat->getObject(first); // Num.wires test: Int_t nWr = fClst1->getNDrTimes(); if(fClst2) nWr += fClst2->getNDrTimes(); nWiresMinTr = (nWiresMin*nMods[fClst1->getSec()])/4; if(nWr < nWiresMinTr) continue; Bool_t flag=kFALSE; if(fClst2==0) { if(fClst1->getNDrTimesMod(0)==0 || fClst1->getNDrTimesMod(1)==0) continue; Int_t typeClFn = fClst1->getTypeClFinder(); if(typeClFn==0) flag = fitSeg(fClst1); // Two mdc per segment else if(typeClFn==2) flag = fitMixedClus(fClst1); if(!flag) continue; fitter=fittersArr[0].getFitter(0); // ??? if(fitter->getSegIndex()<0) continue; fillTrkCandISeg(); } else { flag = fitSec(fClst1,fClst2); //?????????????? } if(flag && fitter->getFitStatus()) sendToMille(); } } void HMdcMille::fitAlgForMilleCosmic(void) { // Two sector fit // Magnet off + combined clusters in sector fitpar.setChi2PerNdfCut(100.); HMdcClus* fClst = NULL; iterClus->Reset(); HMdcClus *fClstArr[6][2]; // [sector][0,1 -IOseg] Int_t nSegs[6]; Int_t nMods[6]; Int_t secList[4]; Int_t nLaySec[6]; Int_t &s1 = secList[0]; Int_t &s2 = secList[1]; Int_t &s3 = secList[2]; Int_t &s4 = secList[3]; // In cosmic one track per event only! First cluster finding: while((fClst=(HMdcClus*)iterClus->Next())) if(fClst->getIndexParent()<0) { Int_t nSectors = 0; Int_t nSegments = 1; Int_t nModules = 0; Int_t nWires = 0; for(Int_t s=0;s<6;s++) { fClstArr[s][0] = NULL; fClstArr[s][1] = NULL; nSegs[s] = 0; nMods[s] = 0; nLaySec[s] = 0; } while(kTRUE) { Int_t sec = fClst->getSec(); Int_t seg = fClst->getIOSeg(); if(fClstArr[sec][seg] != 0) { Error("fitAlgForMilleCosmic","Second time HMdcClus for sec.%i seg.%i",sec,seg); return; } fClstArr[sec][seg] = fClst; if(nSegs[sec] == 0) { secList[nSectors] = sec; nSectors++; } nLaySec[sec] += fClst->getNLayers(); nSegs[sec]++; nSegments++; Int_t nm = fClst->getActiveModule()<0; nMods[sec] += nm; nModules += nm; nWires += fClst->getNDrTimes(); Int_t first,last; fClst->getIndexRegChilds(first,last); if(first < 0 || nSectors == 4) break; fClst = (HMdcClus*)fClusCat->getObject(first); } for(Int_t s=0;s<6;s++) if(fClstArr[s][0]==NULL && fClstArr[s][1]!=NULL) { fClstArr[s][0] = fClstArr[s][1]; fClstArr[s][1] = NULL; } if(nModules==1) continue; if(nModules>=4) nWiresMinTr = nWiresMin; else nWiresMinTr = (nWiresMin*nModules)/4; if(nWires < nWiresMinTr) continue; // Num.wires test fittersArr[0].reset(); Bool_t flag=kFALSE; if(nSegments == 1) { if(nModules < 2) continue; Int_t typeClFn = fClst->getTypeClFinder(); if(typeClFn==0) flag = fitSeg(fClst); else if(typeClFn==2) flag = fitMixedClus(fClst); } else { if (nSectors == 1) flag = fitSec(fClstArr[s1][0],fClstArr[s1][1]); else if(nSectors == 2) flag = fit2Sectors(fClstArr[s1][0],fClstArr[s1][1], fClstArr[s2][0],fClstArr[s2][1]); // else if(nSectors > 2) { //??????????? // Int_t sec1 = 0; // for(Int_t s=0;s<6;s++) if(s!=sec1 && nLaySec[s]>nLaySec[sec1]) sec1 = s; // Int_t sec2 = sec1>0 ? 0 : 1; // for(Int_t s=0;s<6;s++) if(s!=sec1 && s!=sec2 && nLaySec[s]>nLaySec[sec2]) sec2 = s; // flag = fit2Sectors(fClstArr[sec1][0],fClstArr[sec1][1],fClstArr[sec2][0],fClstArr[sec2][1]); // } else if(nSectors == 3) flag = fitNSectors(fClstArr[s1][0],fClstArr[s1][1], fClstArr[s2][0],fClstArr[s2][1], fClstArr[s3][0],fClstArr[s3][1],NULL,NULL); else if(nSectors == 4) flag = fitNSectors(fClstArr[s1][0],fClstArr[s1][1], fClstArr[s2][0],fClstArr[s2][1], fClstArr[s3][0],fClstArr[s3][1], fClstArr[s4][0],fClstArr[s4][1]); } if(flag && fitter->getFitStatus()) sendToMille(); // Int_t first,last; // Int_t sec = fClst1->getSec(); // fClst1->getIndexRegChilds(first,last); // HMdcClus* fClst2 = (first<0) ? NULL : (HMdcClus*)fClusCat->getObject(first); // HMdcClus* fClst3 = NULL; // HMdcClus* fClst4 = NULL; // if(fClst2) { // if(sec == fClst2->getSec()) { // fClst2->getIndexRegChilds(first,last); // fClst3 = (first<0) ? 0:(HMdcClus*)fClusCat->getObject(first); // } else { // fClst3 = fClst2; // fClst2 = NULL; // } // if(fClst3) { // fClst3->getIndexRegChilds(first,last); // fClst4 = (first<0) ? 0:(HMdcClus*)fClusCat->getObject(first); // } // } // Int_t sec2 = -1; // if(fClst3) sec2 = fClst3->getSec(); // else if(fClst4) sec2 = fClst4->getSec(); // Int_t nModules = 0; // nModules += fClst1->getActiveModule()<0 ? 2:1; // if(fClst2) nModules += fClst2->getActiveModule()<0 ? 2:1; // if(fClst3) nModules += fClst3->getActiveModule()<0 ? 2:1; // if(fClst4) nModules += fClst4->getActiveModule()<0 ? 2:1; // if(nModules==1) continue; // if(nModules>=4) nWiresMinTr = nWiresMin; // else nWiresMinTr = (nWiresMin*nModules)/4; // // // Num.wires test: // Int_t nWr = fClst1->getNDrTimes(); // if(fClst2) nWr += fClst2->getNDrTimes(); // if(fClst3) nWr += fClst3->getNDrTimes(); // if(fClst4) nWr += fClst4->getNDrTimes(); // if(nWr < nWiresMinTr) continue; // // Bool_t flag=kFALSE; // if(fClst2==0 && fClst3==0 && fClst4==0) { // if(fClst1->getNDrTimesMod(0)==0 || fClst1->getNDrTimesMod(1)==0) continue; // Int_t typeClFn = fClst1->getTypeClFinder(); // if(typeClFn==0) flag = fitSeg(fClst1); // else if(typeClFn==2) flag = fitMixedClus(fClst1); // } else { // if(fClst3==0&&fClst4==0) flag = fitSec(fClst1,fClst2); // else flag = fit2Sectors(fClst1,fClst2,fClst3,fClst4); // } // if(flag && fitter->getFitStatus()) sendToMille(); } } void HMdcMille::sendToMille(void) { if(mille == 0) return; finalParam = fitter->getFinalParam(); HMdcWiresArr& wiresArr = fitter->getWiresArr(); Int_t nWiresData = wiresArr.getNDriftTimes(); // test chi2 and ...: // if(fitter->getMaxChi2()>9) return; // TEST TEST ???????????????!!!!!!!!!!!!!!!!!!!!!!!!!!! HMdcList24GroupCells& listCells = wiresArr.getOutputListCells(); for(Int_t mod=0;mod<4;mod++) { Int_t nLayersMod = listCells.getNLayersMod(mod); if(nLayersMod == 1) return; //??? if(mod==1 && nLayersMod < 4) return; } if(wiresArr.getSector2()>=0) { HMdcList24GroupCells &listCells2 = wiresArr.getOutputListCells2(); for(Int_t mod=0;mod<4;mod++) if(listCells2.getNLayersMod(mod) == 1) return; if(wiresArr.getSector3()>=0) { HMdcList24GroupCells &listCells3 = wiresArr.getOutputListCells3(); for(Int_t mod=0;mod<4;mod++) if(listCells3.getNLayersMod(mod) == 1) return; if(wiresArr.getSector4()>=0) { HMdcList24GroupCells &listCells4 = wiresArr.getOutputListCells4(); for(Int_t mod=0;mod<4;mod++) if(listCells4.getNLayersMod(mod) == 1) return; } } } for(Int_t wInd=0;wIndgetHitStatus() != 1) continue; if(!wireData->isInCell()) return; } p1.setXYZ(finalParam->x1(),finalParam->y1(),finalParam->z1()); p2.setXYZ(finalParam->x2(),finalParam->y2(),finalParam->z2()); Int_t nWires = finalParam->getNumOfGoodWires(); Int_t nLayers = wiresArr.getOutputNLayers(); if(nWires < nWiresMinTr) return; // Filling histograms: if(doHists) { Int_t hstInd = histInd[wiresArr.getSector()]; if(hstInd>=0) { ((TH2F*)fHistograms->At(hstInd))->Fill(finalParam->getPhiDeg(),finalParam->getThetaDeg()); ((TH1F*)fHistograms->At(hstInd+1))->Fill(finalParam->getChi2()); ((TH1F*)fHistograms->At(hstInd+2))->Fill(nWires+0.1); ((TH1F*)fHistograms->At(hstInd+3))->Fill(nLayers+0.1); } } // Derivatives calculation: if(fixFitTOffset) { for(Int_t wInd=0;wIndgetAnalytDeriv(locDer,0)) continue; int nGlobDer = getGlobalDer(); mille->mille(4,locDer,nGlobDer,globDer,globLabel, -wireData->getDev(),wireData->getErrTdcTime()); } // nWires } else { for(Int_t mi=0;mi<24;mi++) for(Int_t p=0;pgetHitStatus() != 1) continue; setIModILay(); // calculate iMod and iLay + ... for(Int_t p=0;p<6;p++) tofDerCorr[iMod][p] += calcDerTofCorr(p); if(mShFlag[sector][module][6]) tofDerCorr[iMod][6] += calcDerTofCorr(6); } for(Int_t wInd=0;wIndgetAnalytDeriv(locDer,finalParam)) continue; int nGlobDer = getGlobalDerWTof(); mille->mille(4,locDer,nGlobDer,globDer,globLabel, -wireData->getDev(),wireData->getErrTdcTime()); // if(module==0) printf("%iw.Glob: Der %g %g %g %g %g %g %g %g\n", // wInd,globDer[0],globDer[1],globDer[2],globDer[3],globDer[4],globDer[5], // globDer[6],globDer[7]); } // nWires } nWiresTot += nWires; nLayersTot += nLayers; nTracks++; milleFileSize += mille->end(); if(milleFileSize > 1800000000) openNewBinaryFile(); } void HMdcMille::openNewBinaryFile(void) { // New binary file: if(mille) delete mille; TString newFile(milleFileName); newFile += nBinaryFiles; mille = new Mille(newFile); milleFileSize = 0; nBinaryFiles++; } void HMdcMille::setIModILay(void) { sector = wireData->getSector(); module = wireData->getModule(); layer = wireData->getLayer(); setIModILay(sector,module,layer); // calculate iMod and iLay layPosShCurr = layPosSh[iLay]; layNegShCurr = layNegSh[iLay]; layerPart = pSCLayer[iLay]->getLayerPart(wireData->getCell()); if(layerPart == 1) { layPosShCurr += 20; layNegShCurr += 20; } } Int_t HMdcMille::getGlobalDer(void) { setIModILay(); Int_t nParam = 0; //!!! Int_t nMdcPar = mShFlag[sector][module][6] ? 7 : 6; Int_t nMdcPar = cellThicknFree ? 7 : 6; for(Int_t p=0;pgetWtNorm(); Double_t sumWtNorm = finalParam->getSumWtNorm(wireData->getModIndex()); if(sumWtNorm > 0.) return -dDrT*wtNorm/sumWtNorm; return 0.; } int HMdcMille::packLabel(Int_t s,Int_t m,Int_t l,Int_t parn) { // if lay<0 - it is shift of mdc, otherwise shift of layer if(l<0) l = -1; if(createPedeInFile || !isGeomChanged) { if(l<0) mShFlag[s][m][parn] = kTRUE; else lShFlag[s][m][l][parn] = kTRUE; } int label = (((shiftType*10 + s+1)*10 + m+1)*10 + l+1)*10 + parn+1; if(beamTimeId>0 && l<0 && parn!=6) label += beamTimeId*100000; return label; } Bool_t HMdcMille::unpackLabel(Int_t label,Int_t& btId,Int_t& sys, Int_t& sec,Int_t& mod,Int_t& lay,Int_t& parn) { // shSys = 0 shift in respect sector system // shSys = 1 shift in respect mdc system // if lay==-1 - it is shift of mdc, otherwise shift of layer // parn = 0 - shift along X // parn = 1 - shift along Y // parn = 2 - shift along Z // parn = 3 - rotation around X // parn = 4 - rotation around Y // parn = 5 - rotation around Z // parn = 6 - cell thickness // parn = 7 - shift along Y for second layer part // parn = 8 - rotation around Z for second layer part btId = label/100000; // beam time Id sys = (label%100000)/10000; sec = (label%10000)/1000 - 1; mod = (label%1000)/100 - 1; lay = (label%100)/10 - 1; parn = (label%10) - 1; if(sec<0||sec>5 || mod<0||mod>3 || lay<-1||lay>5 || parn<0) return kFALSE; return kTRUE; } Double_t HMdcMille::calcGlobDer(Int_t parNum) { Double_t dDrT = calcDDriftTime(parNum); if(parNum==6) return dDrT*derLThNorm[layer]; else return dDrT*derNorm[parNum]; // derNorm=1/parSteps } Double_t HMdcMille::calcDDriftTime(Int_t p) { Int_t dSign1 = {0}; Int_t dSign2 = {0}; Double_t drTime1 = calcDriftTime(layPosShCurr[p],dSign1); Double_t drTime2 = calcDriftTime(layNegShCurr[p],dSign2); if(dSign1 == dSign2) return (drTime1-drTime2)*0.5; Int_t dSignT = wireData->getDistanceSign(); if(dSignT==dSign1) return drTime1 - wireData->getDrTime(); return wireData->getDrTime() - drTime2; // dSignT==dSign2 } Double_t HMdcMille::calcDriftTime(const HGeomTransform& laySys,Int_t& distSign) { HGeomVector p1l,p2l; Int_t parSec = finalParam->getSec(); if(parSec==sector || parSec<0) { p1l = p1; p2l = p2; } else { HMdcSizesCells* pSCells = HMdcSizesCells::getObject(); pSCells->transLineToOtherSec(*finalParam,sector,p1l,p2l); } p1l = laySys.transTo(p1l); p2l = laySys.transTo(p2l); return wireData->calcDriftTimeForAlign(p1l,p2l,distSign); } void HMdcMille::setPedeInFileName(const char* fname) { pedeInParFName = fname; if(pedeInParFName.Length() > 0) createPedeInFile = kTRUE; else Error("setPedeInFileName","File name is not specified!!!"); } void HMdcMille::creatPedeInParamFile(void) { FILE* file = fopen(pedeInParFName.Data(),"r"); if(file) { printf("******CreatePedeInParamFile******\n"); printf("* File %s exist!\n",pedeInParFName.Data()); printf("* If you want to change this file\n"); printf("* do it by hand or delete file.\n"); printf("*********************************\n"); fclose(file); return; } file = fopen(pedeInParFName.Data(),"w"); if(file == 0) Error("creatPedeInParamFile","Can't open file %s", pedeInParFName.Data()); else { TString sys; if(shiftType == 0) sys = "All shifts done in sector sys."; if(shiftType == 1) sys = "All shifts done in mdc sys."; fprintf(file," Parameter ! %s\n",sys.Data()); for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) for(Int_t l=-1;l<6;l++) { if(l>=0 && useMdcShParOnly) continue; for(Int_t p=0;p=0 && !lShFlag[s][m][l][p]) continue; if( (l<0 && mFixFlag[s][m][p]) || (l>=0 && lFixFlag[s][m][l][p]) ) fprintf(file,"%11i 0.0 -1.0000 !",packLabel(s,m,l,p)); else fprintf(file,"%11i 0.0 0.0000 !",packLabel(s,m,l,p)); if(l<0) fprintf(file," %is. %im.: mdc",s+1,m+1); else { fprintf(file," %is. %im. %il.: layer",s+1,m+1,l+1); if(p==7 || p==8) fprintf(file," part II"); } fprintf(file," %s\n",shitsInfo[p].Data()); // Print description } } fclose(file); } } Bool_t HMdcMille::loadGeometryPar(void) { HRuntimeDb* rtdb = gHades->getRuntimeDb(); if(rtdb == 0) return kFALSE; HParAsciiFileIo* inputFile = new HParAsciiFileIo; if(!inputFile->open(const_cast(geomParFileName.Data()),"in")) return kTRUE; // No file! HMdcGeomStruct* pMdcGeomStruct = (HMdcGeomStruct*)rtdb->getContainer("MdcGeomStruct"); if( !((HParSet*)pMdcGeomStruct)->init(inputFile) ) return kFALSE; pMdcGeomStruct->setStatic(); pMdcGeomStruct->setInputVersion(1,1); HSpecGeomPar* pSpecGeomPar = (HSpecGeomPar*)rtdb->getContainer("SpecGeomPar"); HMdcGeomPar* pMdcGeomPar = (HMdcGeomPar*)rtdb->getContainer("MdcGeomPar"); HMdcLayerGeomPar* pMdcLGeomPar = (HMdcLayerGeomPar*)rtdb->getContainer("MdcLayerGeomPar"); HMdcLayerCorrPar* pMdcLCorrPar = (HMdcLayerCorrPar*)rtdb->getContainer("MdcLayerCorrPar"); if( !((HParSet*)pSpecGeomPar)->init(inputFile) || !((HParSet*)pMdcGeomPar )->init(inputFile) || !((HParSet*)pMdcLGeomPar)->init(inputFile) || !((HParSet*)pMdcLCorrPar)->init(inputFile)) return kFALSE; pSpecGeomPar->setInputVersion(1,1); pMdcGeomPar ->setInputVersion(1,1); pMdcLGeomPar->setInputVersion(1,1); pMdcLCorrPar->setInputVersion(1,1); pSpecGeomPar->setStatic(); pMdcGeomPar ->setStatic(); pMdcLGeomPar->setStatic(); pMdcLCorrPar->setStatic(); inputFile->close(); delete inputFile; isGeomFileExist = kTRUE; return kTRUE; } void HMdcMille::readPedeResFile(const char* fileName) { if(fileName == 0) printErrorAndExit("readPedeResFile","File name of Pede results not specified!"); // Open file: FILE* file = fopen(fileName,"r"); if(file == 0) { if(iteration==0) return; printErrorAndExit("readPedeResFile","No pede result file!"); } // Read file: TString buffer; if( !buffer.Gets(file) ) printErrorAndExit("readPedeResFile","Wrong format of file %s",fileName); Ssiz_t pos = buffer.Index("Parameter"); if(pos<0 || pos>10) printErrorAndExit("readPedeResFile","Wrong format of file %s",fileName); while(buffer.Gets(file) && buffer.Length()>16 && addShifts(buffer)); fclose(file); } Bool_t HMdcMille::makeShifts(void) { if( !isGeomChanged ) return kTRUE; if( !isGeomFileExist ) return kFALSE; HMdcLayerGeomPar* pMdcLGeomPar = pGetCont->getMdcLayerGeomPar(); HMdcLayerCorrPar* pMdcLCorrPar = pGetCont->getMdcLayerCorrPar(); Bool_t reCalcLayerTransf = kFALSE; HGeomTransform shift; HGeomTransform trans; for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) { if(pMdcDet->getModule(s,m) == 0) continue; HMdcLayerGeomParMod& fLayerParMod = (*pMdcLGeomPar)[s][m]; if(!useMdcShParOnly) for(Int_t l=0;l<6;l++) { if(shiftSysLay[s][m][l] != 1 ) continue; //!!!!!!!!!!! // Tol'ko dlya Y seychas!!!!!!!!!!!!!!!!!!!!!!!! //V kakoy posledovatel'nosti???? HMdcLayerGeomParLay& fLayerParLay = fLayerParMod[l]; Double_t pitch = fLayerParLay.getPitch(); Double_t centWrOld = fLayerParLay.getCentralWireNr(); Double_t wireOrOld = fLayerParLay.getWireOrient(); Double_t wireOrNew = wireOrOld + shiftsLay[s][m][l][5]; // Layer rotation Double_t wireOrR = wireOrNew*TMath::DegToRad(); Double_t layShift = shiftsLay[s][m][l][1]*TMath::Cos(wireOrR) - shiftsLay[s][m][l][0]*TMath::Sin(wireOrR); // Layer shift Double_t centWrNew = centWrOld - layShift/pitch; fLayerParLay.setWireOrient(wireOrNew); fLayerParLay.setCentralWireNr(centWrNew); // ------- Layer second part: ------- if(m<2 || !doLayP2Align) continue; Int_t firstWrPII; Float_t layShiftPIIold, wrOrCorrPIIold; if(!pMdcLCorrPar->getLayerCorrPar(s,m,l, firstWrPII,layShiftPIIold,wrOrCorrPIIold)) continue; if(firstWrPII>300) pMdcLCorrPar->addLayerShift(s,m,l,firstWrPII,0.,0.); else { // Double_t wireOrPIIold = wireOrOld - wrOrCorrPIIold; // Double_t wireOrPIInew = wireOrPIIold + shiftsLay[s][m][l][8]; // Layer PII rotation // // Double_t yShiftPIIold = layShiftPIIold/TMath::Cos(wireOrPIIold); // Double_t yShiftPIInew = yShiftPIIold + shiftsLay[s][m][l][7]; // Layer PII y-shift // // Float_t wrOrCorrPIInew = wireOrNew - wireOrPIInew; // Float_t layShiftPIInew = yShiftPIInew*TMath::Cos(wireOrPIInew*TMath::DegToRad()) - layShift; Float_t wrOrCorrPIInew = wrOrCorrPIIold + shiftsLay[s][m][l][5] - shiftsLay[s][m][l][8]; Double_t cosWrOrPIInew = TMath::Cos((wireOrNew - wrOrCorrPIInew)*TMath::DegToRad()); Float_t layShiftPIInew = layShiftPIIold + shiftsLay[s][m][l][7]*cosWrOrPIInew - layShift; pMdcLCorrPar->addLayerShift(s,m,l,firstWrPII,layShiftPIInew,wrOrCorrPIInew); } } if(shiftSysMdc[s][m]<0) continue; // Cell thickness: if(shiftsMdc[s][m][6] != 0.) { for(Int_t l=0;l<6;l++) { HMdcLayerGeomParLay& fLayerParLay = fLayerParMod[l]; Float_t catDist = fLayerParLay.getCatDist(); fLayerParLay.setCatDist(catDist + shiftsMdc[s][m][6]); } reCalcLayerTransf = kTRUE; } // Shift MDC: HMdcSizesCells::setTransform(shiftsMdc[s][m],shift); //first 6 par. are used HModGeomPar* fMod = pGetCont->getModGeomPar(s,m); if(shiftSysMdc[s][m]==0) { // shifts in sector coor.system if( !pGetCont->getSecTransMod(trans,s,m) ) return kFALSE; trans.transFrom(shift); trans.transFrom(pGetCont->getLabTransSec(s)); fMod->getLabTransform() = trans; } else if(shiftSysMdc[s][m]==1) { // shifts in mdc coor.system if(!pGetCont->getLabTransMod(trans,s,m)) return kFALSE; shift.transFrom(trans); fMod->getLabTransform() = shift; } } if(reCalcLayerTransf) pMdcLGeomPar->calcLayerTransformations(); if(isGeomChanged) printf("\nMDC geometry has changed !!!\n\n"); return kTRUE; } Bool_t HMdcMille::writeParAsciiFile(void) { TString outGeomFileName(geomParFileName + ".last"); HParAsciiFileIo* outputFile = new HParAsciiFileIo; if( !outputFile->open(const_cast(outGeomFileName.Data()),"out") ) return kFALSE; HRuntimeDb* rtdb = gHades->getRuntimeDb(); if(rtdb == 0) return kFALSE; if(rtdb->getContainer("MdcGeomStruct")->write(outputFile)<0) return kFALSE; if(rtdb->getContainer("SpecGeomPar")->write(outputFile)<0) return kFALSE; if(rtdb->getContainer("MdcGeomPar")->write(outputFile)<0) return kFALSE; if(rtdb->getContainer("MdcLayerGeomPar")->write(outputFile)<0) return kFALSE; if(rtdb->getContainer("MdcLayerCorrPar")->write(outputFile)<0) return kFALSE; // if(rtdb->getContainer("MdcCalParRaw")->write(outputFile)<0) return kFALSE; outputFile->close(); delete outputFile; if(doCopyGeomFile) copyFile("cp",outGeomFileName); if(!isGeomFileExist && !isGeomChanged) { TString cpFile("cp "+outGeomFileName+" "+geomParFileName); gSystem->Exec(cpFile); } return kTRUE; } void HMdcMille::fixFullMod(Int_t s,Int_t m) { for(Int_t p=0;p<7;p++) fixModPar(s,m,p); } void HMdcMille::fixModPar(Int_t s,Int_t m,Int_t p) { if(s>=0&&s<6 && m>=0&&m<4 && p>=0&&p<7) { mFixFlag[s][m][p] = kTRUE; if(sigmaMdc[s][m][p] == 0.) sigmaMdc[s][m][p] = -1.; } } void HMdcMille::fixFullLay(Int_t s,Int_t m,Int_t l) { for(Int_t p=0;p5 || m<0 || m>3 || l<0 || l>5 || p<0 || p>8) return; lFixFlag[s][m][l][p] = kTRUE; if(sigmaLay[s][m][l][p] == 0.) sigmaLay[s][m][l][p] = -1.; } void HMdcMille::fixPar(Int_t s,Int_t m,Int_t l,Int_t p) { if(l<0) { if(p<0) fixFullMod(s,m); else fixModPar(s,m,p); } else { if(p<0) fixFullLay(s,m,l); else fixLayPar(s,m,l,p); } } void HMdcMille::readSumOfShiftsFile(void) { // Open file: FILE* file = fopen(sumOfShiftsFName.Data(),"r"); if(file == 0) return; // Read file: TString buffer; if( !buffer.Gets(file) ) printErrorAndExit("readSumOfShiftsFile", "Wrong format of file %s",sumOfShiftsFName.Data()); // Stop! Ssiz_t pos = buffer.Index("SumOfShifts:"); if(pos<0 || pos>3) printErrorAndExit("readSumOfShiftsFile", "Wrong format of file %s",sumOfShiftsFName.Data()); // Stop! Int_t sys; Int_t n = sscanf(buffer.Data(),"%*s %i %*s %*s %i",&iteration,&sys); if(n!=2 || iteration<0) printErrorAndExit("readSumOfShiftsFile", "Wrong format of file %s",sumOfShiftsFName.Data()); // Stop! iteration++; if( !buffer.Gets(file) ) printErrorAndExit("readSumOfShiftsFile", "Wrong format of file %s",sumOfShiftsFName.Data()); // Stop! pos = buffer.Index("Statistics:"); if(pos<0 || pos>3) printErrorAndExit("readSumOfShiftsFile", "Wrong format of file %s",sumOfShiftsFName.Data()); // Stop! while(buffer.Gets(file) && buffer.Length()>16 && addShifts(buffer)); fclose(file); } Bool_t HMdcMille::addShifts(TString& buffer) { TObjArray* arr = buffer.Tokenize(" "); Bool_t exitFlag = kFALSE; if(arr->GetEntries()>=3) { Int_t label = ((TObjString*)(arr->At(0)))->GetString().Atoi(); Double_t shift = ((TObjString*)(arr->At(1)))->GetString().Atof(); Double_t coef = ((TObjString*)(arr->At(2)))->GetString().Atof(); Int_t btId,sys,sec,mod,lay,parn; if(!unpackLabel(label,btId,sys,sec,mod,lay,parn)) printErrorAndExit("addShifts","Wrong format of input file!!! Stop!"); exitFlag = kTRUE; if((btId == 0 || btId == beamTimeId) && pMdcDet->getModule(sec,mod)) { mdcMods[sec][mod] = 1; if(lay < 0) { shiftsMdc[sec][mod][parn] += shift; if(coef >= 0.) shiftSysMdc[sec][mod] = sys; // Otherwise = -1 if(shift != 0.) isGeomChanged = kTRUE; mShFlag[sec][mod][parn] = kTRUE; sigmaMdc[sec][mod][parn] = coef; } else { shiftsLay[sec][mod][lay][parn] += shift; if(coef >= 0.) shiftSysLay[sec][mod][lay] = sys; if(shift != 0.) isGeomChanged = kTRUE; lShFlag[sec][mod][lay][parn] = kTRUE; sigmaLay[sec][mod][lay][parn] = coef; } } } arr->Delete(); delete arr; return exitFlag; } void HMdcMille::printErrorAndExit(const char* func,const char* form,const char* str) { #warning the following warning because of string literal can be ignored! if(str==0) Error(func,form); else Error(func,form,str); exit(1); } void HMdcMille::creatSumOfShiftsFile(void) { FILE* file = fopen(sumOfShiftsFName.Data(),"w"); if(file == 0) printErrorAndExit("creatSumOfShiftsFile", "Can't open file %s",sumOfShiftsFName.Data()); TString sys; if(shiftType == 0) sys = "all shifts done in sector sys."; if(shiftType == 1) sys = "all shifts done in mdc sys."; Double_t mNWr = Double_t(nWiresTot)/nTracks; Double_t mNLy = Double_t(nLayersTot)/nTracks; fprintf(file,"SumOfShifts: %i iteration. System %i - %s\n",iteration,shiftType,sys.Data()); fprintf(file,"Statistics: %i tracks; %8.4f ; %8.4f \n",nTracks,mNWr,mNLy); for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) for(Int_t l=-1;l<6;l++) { if(l>=0 && useMdcShParOnly) continue; for(Int_t p=0;p=0 && !lShFlag[s][m][l][p]) continue; fprintf(file,"%11i",packLabel(s,m,l,p)); // Print label if(l < 0) fprintf(file," %18.10e %9.4e ! %is. %im.: mdc", shiftsMdc[s][m][p],sigmaMdc[s][m][p],s+1,m+1); else { fprintf(file," %18.10e %9.4e ! %is. %im. %il.: layer", shiftsLay[s][m][l][p],sigmaLay[s][m][l][p],s+1,m+1,l+1); if(p==7 || p==8) fprintf(file," part II"); } fprintf(file," %s\n",shitsInfo[p].Data()); // Print description } } fclose(file); if(doCopySumShFile) copyFile("cp",sumOfShiftsFName); if(doCopyResFile) copyFile("cp",pedeResFName,".res"); if(doCopyHisFile) copyFile("mv",pedeResFName,".his"); if(doCopyLogFile) copyFile("mv",pedeResFName,".log"); } void HMdcMille::copyFile(const char* op,TString& file,const char* ext) { TString cpFile(op); if(ext) cpFile += " "+file+ext+" "+file+ext+stepIter; else cpFile += " "+file+" "+file+stepIter; gSystem->Exec(cpFile); } void HMdcMille::creatPedeTaskFile(const char* fileName) { pedeTaskFileName = fileName; if(pedeTaskFileName.Length()<=0) printErrorAndExit("creatPedeTaskFile", "File name of pede task list is not specified!"); crPedeTaskFile = kTRUE; // Set default method: method = "inversion "; numOfIter = 3; accuracy = 0.001; mthDescr = "Gauss matrix inversion"; bandwidth = -1; } void HMdcMille::useInversionMethod(Int_t nit,Float_t acc) { method = "inversion "; numOfIter = nit; accuracy = acc; mthDescr = "Gauss matrix inversion"; bandwidth = -1; } void HMdcMille::useBandCholeskyMethod(Int_t nit,Float_t acc,Int_t bwidth) { method = "bandcholesky "; numOfIter = nit; accuracy = acc; mthDescr = "Cholesky"; bandwidth = bwidth; } void HMdcMille::useCholeskyMethod(Int_t nit,Float_t acc,Int_t bwidth) { method = "cholesky "; numOfIter = nit; accuracy = acc; mthDescr = "Cholesky"; bandwidth = bwidth; } void HMdcMille::useSparseGMRESMethod(Int_t nit,Float_t acc) { method = "sparseGMRES "; numOfIter = nit; accuracy = acc; mthDescr = "minimal residual"; bandwidth = -1; } void HMdcMille::useFullGMRESMethod(Int_t nit,Float_t acc) { method = "fullGMRES "; numOfIter = nit; accuracy = acc; mthDescr = "minimal residual"; bandwidth = -1; } void HMdcMille::useDiagonalizationMethod(Int_t nit,Float_t acc) { method = "diagonalization"; numOfIter = nit; accuracy = acc; mthDescr = "diagonalization"; bandwidth = -1; } void HMdcMille::setParConstrainFile(const char* file) { parConstrainFile = file; } void HMdcMille::writePedeTaskFile(void) { if(!crPedeTaskFile) return; if(pedeInParFName.Length() <= 0) printErrorAndExit("writePedeTaskFile", "File name of parameters list not specified!"); FILE* file = fopen(pedeTaskFileName.Data(),"w"); if(file==0) printErrorAndExit("writePedeTaskFile","Can't open file %s",pedeTaskFileName.Data()); fprintf(file,"Cfiles !\n"); if(parConstrainFile.Length() > 0) fprintf(file,"%s ! constraints text file\n",parConstrainFile.Data()); else fprintf(file,"!parConstrain.txt ! constraints text file\n"); fprintf(file,"%s ! parameter text file\n",pedeInParFName.Data()); for(Int_t nf=0;nf 0) fprintf(file,"bandwidth %i ! width of band matrix\n",bandwidth); fprintf(file,"method %s %i %g ! %s\n",method.Data(),numOfIter,accuracy,mthDescr.Data()); fclose(file); } void HMdcMille::setMaxNumWiresCut(Int_t* wc) { for(Int_t s=0;s<6;s++) nWiresCut[s] = wc[s]; } void HMdcMille::setMaxNumWiresCut(Int_t cs1,Int_t cs2,Int_t cs3,Int_t cs4,Int_t cs5,Int_t cs6) { nWiresCut[0] = cs1; nWiresCut[1] = cs2; nWiresCut[2] = cs3; nWiresCut[3] = cs4; nWiresCut[4] = cs5; nWiresCut[5] = cs6; }