Calculate crosstalk correction in separate function. Do calculation in 2 itterations...

diff --git a/TPC/Rec/AliTPCtracker.cxx b/TPC/Rec/AliTPCtracker.cxx
index 7314ded998640ff138a15ce027e02ad081b30f36..d61293395d9ca136ee570142ee0cf356a2218a72 100644 (file)
--- 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; jmatrix<nTimeBinsAll; jmatrix++){
@@ -1345,19 +1345,12 @@ Int_t  AliTPCtracker::LoadClusters(const TClonesArray *arr)
 
 Int_t  AliTPCtracker::LoadClusters()
 {
-  //
-  // load clusters to the memory  
-  Int_t nROCs   = 72;
+ //
+  // load clusters to the memory
   static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
   //
   //  TTree * tree = fClustersArray.GetTree();
   AliInfo("LoadClusters()\n");
-  // reset crosstalk matrix
-  //
-  for (Int_t isector=0; isector<nROCs*2; isector++){  //set all ellemts of crosstalk matrix to 0
-    TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(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; i<j; i++) {
     br->GetEntry(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; icl<clrow->GetArray()->GetEntriesFast(); icl++){
-      AliTPCclusterMI *clXtalk= static_cast<AliTPCclusterMI*>(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; 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 &sector= (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
     //
-    // 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; isector<nROCs*4; isector++){  //set all ellemts of crosstalk matrix to 0 
+      TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(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; isector<nROCs; isector++){  //set all ellemts of crosstalk matrix to 0
       TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(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="<<eventNr<<
       "counterAll="<<counterAll<<
       "counterOut="<<counterOut<<
+      "matSectotCluster.="<<&matSectorCluster<<                   // 
       //
       "filteredSector="<<filteredSector<<                        //  counter filtered sectors                   
       "filteredSectorTime="<<filteredSectorTime<<                //  counter filtered time bins

diff --git a/TPC/Rec/AliTPCtracker.h b/TPC/Rec/AliTPCtracker.h
index 1330733122d9184a385e1470ca00cc5a8d110062..dce6b1351b53e19230b8982a9d20d6b4f892959a 100644 (file)
--- a/TPC/Rec/AliTPCtracker.h
+++ b/TPC/Rec/AliTPCtracker.h
@@ -84,6 +84,7 @@ public:
   void   Transform(AliTPCclusterMI * cluster);
   void ApplyTailCancellation();
   void ApplyXtalkCorrection();
+  void CalculateXtalkCorrection();
   void GetTailValue(Float_t ampfactor,Double_t &ionTailMax,Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1);
   //
   void FillESD(const TObjArray* arr);

diff --git a/test/testdEdx/rec.C b/test/testdEdx/rec.C
index 7e0900edcb2e7f34c798fbebc2baf0605eb310cb..de2bd5fa75ce0a07721a6dbd9217c079f8f56620 100644 (file)
--- a/test/testdEdx/rec.C
+++ b/test/testdEdx/rec.C
@@ -55,6 +55,7 @@ void rec(Bool_t useIonTail, Double_t crossTalkCorrection) {
   //
   // Run reconstruction
   //
+  AliTPCReconstructor::SetStreamLevel(AliTPCtracker::kStreamCrosstalkMatrix | AliTPCtracker::kStreamXtalk  | AliTPCtracker::kStreamIonTail);  // enable steaming of information for ion tail cancelation and crosstalk
   TStopwatch timer;
   timer.Start();
   reco.Run();