From 3c8c25a3ca426bf268f989f393980e3c5e71bada Mon Sep 17 00:00:00 2001 From: mivanov Date: Tue, 4 Nov 2014 17:05:23 +0100 Subject: [PATCH] Calculate crosstalk correction in separate function. Do calculation in 2 itterations.SetStreamLevel to enable streaming of additional information --- TPC/Rec/AliTPCtracker.cxx | 201 ++++++++++++++++++++++++-------------- TPC/Rec/AliTPCtracker.h | 1 + test/testdEdx/rec.C | 1 + 3 files changed, 127 insertions(+), 76 deletions(-) diff --git a/TPC/Rec/AliTPCtracker.cxx b/TPC/Rec/AliTPCtracker.cxx index 7314ded9986..d61293395d9 100644 --- a/TPC/Rec/AliTPCtracker.cxx +++ b/TPC/Rec/AliTPCtracker.cxx @@ -472,9 +472,9 @@ AliTracker(), Int_t nROCs = 72; Int_t nTimeBinsAll = par->GetMaxTBin(); Int_t nWireSegments = 11; - fCrossTalkSignalArray = new TObjArray(nROCs*2); // + fCrossTalkSignalArray = new TObjArray(nROCs*4); // fCrossTalkSignalArray->SetOwner(kTRUE); - for (Int_t isector=0; isector<2*nROCs; isector++){ + for (Int_t isector=0; isector<4*nROCs; isector++){ TMatrixD * crossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll); for (Int_t imatrix = 0; imatrix<11; imatrix++) for (Int_t jmatrix = 0; jmatrixAt(isector); - if (crossTalkMatrix)(*crossTalkMatrix)*=0; - } TTree * tree = fInput; TBranch * br = tree->GetBranch("Segment"); @@ -1366,61 +1359,15 @@ Int_t AliTPCtracker::LoadClusters() // Conversion of pad, row coordinates in local tracking coords. // Could be skipped here; is already done in clusterfinder - Int_t j=Int_t(tree->GetEntries()); for (Int_t i=0; iGetEntry(i); // Int_t sec,row; fkParam->AdjustSectorRow(clrow->GetID(),sec,row); - - // wire segmentID and nPadsPerSegment to be used for Xtalk calculation - Int_t wireSegmentID = fkParam->GetWireSegment(sec,row); - Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID)); - TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec)); - TMatrixD &crossTalkSignalBelow = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs)); - Int_t nCols=crossTalkSignal.GetNcols(); - Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge(); for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ - AliTPCclusterMI *clXtalk= static_cast(clrow->GetArray()->At(icl)); - Transform((AliTPCclusterMI*)(clXtalk)); - Int_t timeBinXtalk = clXtalk->GetTimeBin(); - Double_t rmsPadMin2=0.5*0.5+(fkParam->GetDiffT()*fkParam->GetDiffT())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); // minimal PRF width - 0.5 is the PRF in the pad units - we should et it from fkparam getters - Double_t rmsTimeMin2=1+(fkParam->GetDiffL()*fkParam->GetDiffL())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetZWidth()*fkParam->GetZWidth()); // minimal PRF width - 1 is the TRF in the time bin units - we should et it from fkParam getters - Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth()); - Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); - if (rmsPadMin2>rmsPad2){ - rmsPad2=rmsPadMin2; - } - if (rmsTimeMin2>rmsTime2){ - rmsTime2=rmsTimeMin2; - } - - Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.; - Double_t qTotXtalk = 0.; - Double_t qTotXtalkMissing = 0.; - for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) { - if (itb<0 || itb>=nCols) continue; - Double_t missingCharge=0; - Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2)); - if (missingChargeFactor>0) { - for (Int_t dpad=-2; dpad<=2; dpad++){ - Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf; - qPad-=2.*crossTalkSignal[wireSegmentID][itb]; // missing part due crosttalk - feedback part - factor 2 used - if (TMath::Nint(qPad)<=fkParam->GetZeroSup()){ - missingCharge+=qPad+2.*crossTalkSignal[wireSegmentID][itb]; - }else{ - missingCharge+=2.*crossTalkSignal[wireSegmentID][itb]; - } - } - } - qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor; - qTotXtalkMissing = missingCharge; - crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment; - crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment; - } + Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl))); } - // AliTPCtrackerRow * tpcrow=0; Int_t left=0; @@ -1443,12 +1390,123 @@ Int_t AliTPCtracker::LoadClusters() tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k))); } } - if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) { + // + clrow->Clear("C"); + LoadOuterSectors(); + LoadInnerSectors(); + + cout << " =================================================================================================================================== " << endl; + cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl; + cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl; + cout << " =================================================================================================================================== " << endl; + + if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation(); + if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) CalculateXtalkCorrection(); + if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection(); + //if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters(); + return 0; +} + + +void AliTPCtracker::CalculateXtalkCorrection(){ + // + // + // + const Int_t nROCs = 72; + + // 0.) reset crosstalk matrix + // + for (Int_t isector=0; isectorAt(isector); + if (crossTalkMatrix)(*crossTalkMatrix)*=0; + } + + // + // 1.) Filling part -- loop over clusters + // + Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge(); + for (Int_t iter=0; iter<2;iter++){ + for (Int_t isector=0; isector<36; isector++){ // loop over sectors + for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C + AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18]; + Int_t nrows = sector.GetNRows(); + Int_t sec=0; + if (isector<18) sec=isector+18*iside; + if (isector>=18) sec=36+isector+18*iside; + for (Int_t row = 0;rowGetWireSegment(sec,row); + Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID)); + TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec)); + TMatrixD &crossTalkSignalCache = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs*2)); // this is the cache value of the crosstalk from previous iteration + TMatrixD &crossTalkSignalBelow = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs)); + Int_t nCols=crossTalkSignal.GetNcols(); + // + AliTPCtrackerRow& tpcrow = sector[row]; + Int_t ncl = tpcrow.GetN1(); // number of clusters in the row + if (iside>0) ncl=tpcrow.GetN2(); + for (Int_t i=0;i0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i)); + + Int_t timeBinXtalk = clXtalk->GetTimeBin(); + Double_t rmsPadMin2=0.5*0.5+(fkParam->GetDiffT()*fkParam->GetDiffT())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); // minimal PRF width - 0.5 is the PRF in the pad units - we should et it from fkparam getters + Double_t rmsTimeMin2=1+(fkParam->GetDiffL()*fkParam->GetDiffL())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetZWidth()*fkParam->GetZWidth()); // minimal PRF width - 1 is the TRF in the time bin units - we should et it from fkParam getters + Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth()); + Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); + if (rmsPadMin2>rmsPad2){ + rmsPad2=rmsPadMin2; + } + if (rmsTimeMin2>rmsTime2){ + rmsTime2=rmsTimeMin2; + } + + Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.; + Double_t qTotXtalk = 0.; + Double_t qTotXtalkMissing = 0.; + for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) { + if (itb<0 || itb>=nCols) continue; + Double_t missingCharge=0; + Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2)); + if (missingChargeFactor>0) { + for (Int_t dpad=-2; dpad<=2; dpad++){ + Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf; + if (TMath::Nint(qPad-crossTalkSignalCache[wireSegmentID][itb])<=fkParam->GetZeroSup()){ + missingCharge+=qPad+crossTalkSignalCache[wireSegmentID][itb]; + }else{ + missingCharge+=crossTalkSignalCache[wireSegmentID][itb]; + } + } + } + qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor; + qTotXtalkMissing = missingCharge; + crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment; + crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment; + } // end of time bin loop + } // end of cluster loop + } // end of rows loop + } // end of side loop + } // end of sector loop + // + // copy crosstalk matrix to cached used for next itteration // - // dump the crosstalk matrices to tree for further investigation - // a.) to estimate fluctuation of pedestal in indiviula wire segments - // b.) to check correlation between regions - // c.) to check relative conribution of signal below threshold to crosstalk + for (Int_t isector=0; isectorAt(isector); + TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2); + if (crossTalkMatrix){ + (*crossTalkMatrixCache)*=0; + (*crossTalkMatrixCache)+=(*crossTalkMatrix); + } + } + } // end of 2 iterations + + // + // 2.) dump the crosstalk matrices to tree for further investigation + // a.) to estimate fluctuation of pedestal in indiviula wire segments + // b.) to check correlation between regions + // c.) to check relative conribution of signal below threshold to crosstalk + + if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) { for (Int_t isector=0; isectorAt(isector); TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs); @@ -1471,22 +1529,12 @@ Int_t AliTPCtracker::LoadClusters() } - // - clrow->Clear("C"); - LoadOuterSectors(); - LoadInnerSectors(); - cout << " =================================================================================================================================== " << endl; - cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl; - cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl; - cout << " =================================================================================================================================== " << endl; - - if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation(); - if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection(); - //if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters(); - return 0; } + + + void AliTPCtracker::FilterOutlierClusters(){ // // filter outlier clusters @@ -1716,6 +1764,7 @@ void AliTPCtracker::FilterOutlierClusters(){ "eventNr="<