]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITSClusterFinderV2SSD.cxx
commented define _ClusterTopology_ - to be used only for the special productions
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
index 7f7427a27c38ec870af121196edab19898af7bef..2d8c9716bb51d8bbde86b8b411ca1a6211d08099 100644 (file)
 //                                                                        //
 ///////////////////////////////////////////////////////////////////////////
 
+#include "AliITSClusterFinderV2SSD.h"
+
 #include <Riostream.h>
-#include "AliLog.h"
+#include <TGeoGlobalMagField.h>
 
-#include "AliITSClusterFinderV2SSD.h"
+#include "AliLog.h"
+#include "AliMagF.h"
 #include "AliITSRecPoint.h"
 #include "AliITSRecPointContainer.h"
 #include "AliITSgeomTGeo.h"
-#include "AliITSDetTypeRec.h"
 #include "AliRawReader.h"
 #include "AliITSRawStreamSSD.h"
 #include <TClonesArray.h>
@@ -66,15 +68,37 @@ const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
 ClassImp(AliITSClusterFinderV2SSD)
 
 
-AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
-                                                                            fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1)
+  AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), fLorentzShiftP(0), fLorentzShiftN(0)
 {
 //Default constructor
+  static AliITSRecoParam *repa = NULL;  
+  if(!repa){
+    repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
+    if(!repa){
+      repa = AliITSRecoParam::GetHighFluxParam();
+      AliWarning("Using default AliITSRecoParam class");
+    }
+  }
 
+  if (repa->GetCorrectLorentzAngleSSD()) {
+    AliMagF* field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField()); 
+    if (field == 0) {
+      AliError("Cannot get magnetic field from TGeoGlobalMagField");
+    }
+    else {
+      Float_t bField = field->SolenoidField();
+      // NB: spatial shift has opposite sign for lay 5 and 6, but strip numbering also changes direction, so no sign-change 
+      // Shift due to ExB on drift N-side, units: strip width 
+      fLorentzShiftP = -repa->GetTanLorentzAngleElectronsSSD() * 150.e-4/95.e-4 * bField / 5.0;
+      // Shift due to ExB on drift P-side, units: strip width 
+      fLorentzShiftN = -repa->GetTanLorentzAngleHolesSSD() * 150.e-4/95.e-4 * bField / 5.0;
+      AliDebug(1,Form("bField %f Lorentz Shift P-side %f N-side %f",bField,fLorentzShiftN,fLorentzShiftP));
+    }
+  }
 }
  
 //______________________________________________________________________
-AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf),                                              fLastSSD1(cf.fLastSSD1)
+AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1), fLorentzShiftP(cf.fLorentzShiftP), fLorentzShiftN(cf.fLorentzShiftN)
 {
   // Copy constructor
 }
@@ -161,8 +185,11 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
   cout<<"this is module "<<fModule;
   cout<<endl;
   cout<<endl;
+  Int_t layer = 4;
+  if (fModule>fLastSSD1) 
+    layer = 5;
   */
-
   //--------------------------------------------------------
   // start 1D-clustering from the first digit in the digits array
   //
@@ -221,7 +248,14 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
 
        if(flag5) {
          //cout<<"here1"<<endl;
-         c[*n].SetY(y/q);
+         Float_t dLorentz = 0;
+         if (!flag) { // P-side is neg clust
+           dLorentz = fLorentzShiftN;
+         }
+         else { // N-side is p clust
+           dLorentz = fLorentzShiftP;
+         }
+         c[*n].SetY(y/q+dLorentz);
          c[*n].SetQ(q);
          c[*n].SetNd(nd);
         CheckLabels2(milab);
@@ -232,14 +266,14 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
 
           //Split suspiciously big cluster
           if (nd>4&&nd<25) {
-            c[*n].SetY(y/q-0.25*nd);
+            c[*n].SetY(y/q-0.25*nd+dLorentz);
             c[*n].SetQ(0.5*q);
             (*n)++;
             if (*n==kMax) {
               Error("FindClustersSSD","Too many 1D clusters !");
               return;
             }
-            c[*n].SetY(y/q+0.25*nd);
+            c[*n].SetY(y/q+0.25*nd+dLorentz);
             c[*n].SetQ(0.5*q);
             c[*n].SetNd(nd);
             c[*n].SetLabels(milab);
@@ -311,8 +345,15 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
   if(flag5) {
 
     // cout<<"here2"<<endl;
-
-    c[*n].SetY(y/q);
+    Float_t dLorentz = 0;
+    if (!flag) { // P-side is neg clust
+      dLorentz = fLorentzShiftN;
+    }
+    else { // N-side is p clust
+      dLorentz = fLorentzShiftP;
+    }
+    
+    c[*n].SetY(y/q + dLorentz);
     c[*n].SetQ(q);
     c[*n].SetNd(nd);
     c[*n].SetLabels(lab);
@@ -321,14 +362,14 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
       
       //Split suspiciously big cluster
       if (nd>4 && nd<25) {
-       c[*n].SetY(y/q-0.25*nd);
+       c[*n].SetY(y/q-0.25*nd + dLorentz);
        c[*n].SetQ(0.5*q);
        (*n)++;
        if (*n==kMax) {
          Error("FindClustersSSD","Too many 1D clusters !");
          return;
        }
-       c[*n].SetY(y/q+0.25*nd);
+       c[*n].SetY(y/q+0.25*nd + dLorentz);
        c[*n].SetQ(0.5*q);
        c[*n].SetNd(nd);
        c[*n].SetLabels(lab);
@@ -367,9 +408,10 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
 
 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader){
 
-    //------------------------------------------------------------
+  //------------------------------------------------------------
   // This function creates ITS clusters from raw data
   //------------------------------------------------------------
+  fNClusters = 0;
   rawReader->Reset();
   AliITSRawStreamSSD inputSSD(rawReader);
   FindClustersSSD(&inputSSD);
@@ -392,11 +434,15 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
       AliWarning("Using default AliITSRecoParam class");
     }
   }
+  if (fRawID2ClusID) { // RS: reset references from 1D clusters to rawID's
+    fRawIDRef[0].Reset();
+    fRawIDRef[1].Reset();
+  }
   Int_t nClustersSSD = 0;
   const Int_t kNADC = 12;
   const Int_t kMaxADCClusters = 1000;
 
-  Int_t strips[kNADC][2][kMaxADCClusters][2]; // [ADC],[side],[istrip], [0]=istrip [1]=signal
+  Int_t strips[kNADC][2][kMaxADCClusters][3]; // [ADC],[side],[istrip], [0]=istrip [1]=signal [2]=rawID (for embedding, RS)
   Int_t nStrips[kNADC][2];
 
   for( int i=0; i<kNADC; i++ ){
@@ -410,7 +456,9 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
   //*
   //* Loop over modules DDL+AD
   //*
-  
+  int countRW = 0; //RS
+  if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
+
   while (kTRUE) {
 
     bool next = input->Next();
@@ -461,6 +509,11 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
 //CM channels are always present even everything is suppressed 
          continue;
        }
+       
+       /*      Int_t layer = 4;
+       if (fModule>fLastSSD1) 
+         layer = 5;
+       */
 
        AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
        if( !cal ){
@@ -486,18 +539,26 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
          Int_t nDigits = 0;
          Int_t ostrip = -2;
          Bool_t snFlag = 0;
+
+         Float_t dLorentz = 0;
+         if (side==0) { // P-side is neg clust
+           dLorentz = fLorentzShiftN;
+         }
+         else { // N-side is pos clust
+           dLorentz = fLorentzShiftP;
+         }
          
          Int_t n = nStrips[adc][side];
          for( int istr = 0; istr<n+1; istr++ ){
            
            bool stripOK = 1;
-           Int_t strip=0;
+           Int_t strip=0, rwID = 0;
            Float_t signal=0.0, noise=0.0, gain=0.0;
            
            if( istr<n ){
              strip = strips[adc][side][istr][0];
              signal = strips[adc][side][istr][1];
-             
+             rwID   = strips[adc][side][istr][2]; // RS
              //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
 
              if( cal ){
@@ -522,7 +583,7 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
                }else {
                  
                  Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
-                 cluster.SetY( y / q + dStrip);
+                 cluster.SetY( y / q + dStrip + dLorentz);
                  cluster.SetQ(q);
                  cluster.SetNd(nDigits);
                  cluster.SetLabels(lab);
@@ -532,10 +593,10 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
                  if( repa->GetUseUnfoldingInClusterFinderSSD()
                      && nDigits > 4 && nDigits < 25 
                      ){
-                   cluster.SetY(y/q + dStrip - 0.25*nDigits);      
+                   cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);           
                    cluster.SetQ(0.5*q);          
                    Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
-                   cluster2.SetY(y/q + dStrip + 0.25*nDigits);     
+                   cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);          
                    cluster2.SetQ(0.5*q);
                    cluster2.SetNd(nDigits);
                    cluster2.SetLabels(lab);      
@@ -559,6 +620,7 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
              nDigits++;
              //nstat[side]++;
              ostrip = strip;
+             if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
 
            }
          } //* end loop over strips
@@ -639,11 +701,13 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
     }
     strips[adc][side][n][0] = strip;
     strips[adc][side][n][1] = signal;    
+    strips[adc][side][n][2] = countRW;    
     n++;
 
     //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
     //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
-
+    //
+    countRW++; //RS
   } //* End main loop over the input
   
   AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
@@ -651,8 +715,8 @@ void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
 
 
 void AliITSClusterFinderV2SSD::
-FindClustersSSD(Ali1Dcluster* neg, Int_t nn, 
-               Ali1Dcluster* pos, Int_t np,
+FindClustersSSD(const Ali1Dcluster* neg, Int_t nn, 
+               const Ali1Dcluster* pos, Int_t np,
                TClonesArray *clusters) {
   //------------------------------------------------------------
   // Actual SSD cluster finder
@@ -674,7 +738,7 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
 
 //  TClonesArray &cl=*clusters;
   
-  AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
+  AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
   if (fModule>fLastSSD1) 
     seg->SetLayer(6);
   else 
@@ -700,8 +764,8 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
 
   if ((np*nn) > fgPairsSize) {
 
-    if (fgPairs) delete [] fgPairs;
-    fgPairsSize = 4*np*nn;
+    delete [] fgPairs;
+    fgPairsSize = 2*np*nn;
     fgPairs = new Short_t[fgPairsSize];
   }
   memset(fgPairs,0,sizeof(Short_t)*np*nn);
@@ -872,7 +936,7 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
        }
 
        AliITSRecPoint * cl2;
-    cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
+       cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
          
     cl2->SetChargeRatio(ratio);        
     cl2->SetType(1);
@@ -1235,7 +1299,7 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
       //
       if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){ 
        Float_t minchargediff =4.;
-       Float_t minchargeratio =0.2;
+       //      Float_t minchargeratio =0.2;
 
        Int_t j=-1;
        for (Int_t di=0;di<cnegative[ip];di++){
@@ -1246,7 +1310,7 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
          if (TMath::Abs(ratio)<0.2){
            j =jc;
            minchargediff = TMath::Abs(chargedif);
-           minchargeratio = TMath::Abs(ratio);
+           //      minchargeratio = TMath::Abs(ratio);
          }
        }
        if (j<0) continue;  // not proper cluster      
@@ -1453,6 +1517,26 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
            lp[5]=-4.32e-05;
          }
        }
+       // 
+       if (fRawID2ClusID) { // set rawID <-> clusterID correspondence for embedding
+         const int kMaxRefRW = 200;
+         UInt_t nrefsRW,refsRW[kMaxRefRW];
+         nrefsRW = fRawIDRef[0].GetReferences(j,refsRW,kMaxRefRW); // n-side
+         for (int ir=nrefsRW;ir--;) {
+           int rwid = (int)refsRW[ir];
+           if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
+           (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
+         }
+         //
+         nrefsRW = fRawIDRef[1].GetReferences(ip,refsRW,kMaxRefRW); // p-side
+         for (int ir=nrefsRW;ir--;) {
+           int rwid = (int)refsRW[ir];
+           if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
+           (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
+         }
+         //
+         milab[0] = fNClusters+1;  // RS: assign id as cluster label
+       }
 
        AliITSRecPoint * cl2;
          cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
@@ -1466,7 +1550,8 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
          }
          cused1[ip]++;
          cused2[j]++;      
-       ncl++;
+         ncl++;
+         fNClusters++;
        
       } // manyXmany
     } // loop over Pside 1Dclusters