+ //
+ 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; isector<nROCs*4; isector++){ //set all ellemts of crosstalk matrix to 0
+ TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(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;row<nrows;row++){ // loop over rows
+ //
+ //
+ Int_t wireSegmentID = fkParam->GetWireSegment(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;i<ncl;i++) { // loop over clusters
+ AliTPCclusterMI *clXtalk= (iside>0)?(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