]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HMPID/AliHMPIDRecon.cxx
Completely reengineered version of CMake build system (Johny)
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
index 7af3e5b7c2fbb36749354d747f5614bc7da57ce8..27069b2a598108cf7b39e40a3f5717767f5c8eab 100644 (file)
@@ -60,7 +60,7 @@ void AliHMPIDRecon::InitVars(Int_t n)
 //..
 //Init some variables
 //..
-  if(n<0) return;
+  if(n<=0) return;
   fPhotFlag = new Int_t[n];
   fPhotCkov = new Double_t[n];
   fPhotPhi  = new Double_t[n];
@@ -73,13 +73,13 @@ void AliHMPIDRecon::DeleteVars()const
 //..
 //Delete variables
 //..
-  delete fPhotFlag;
-  delete fPhotCkov;
-  delete fPhotPhi;
-  delete fPhotWei;
+  delete [] fPhotFlag;
+  delete [] fPhotCkov;
+  delete [] fPhotPhi;
+  delete [] fPhotWei;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean,Float_t xRa,Float_t yRa)
 {
 // Pattern recognition method based on Hough transform
 // Arguments:   pTrk     - track for which Ckov angle is to be found
@@ -89,52 +89,49 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t n
 
   const Int_t nMinPhotAcc = 3;                      // Minimum number of photons required to perform the pattern recognition
   
-  
   Int_t nClusTot = pCluLst->GetEntries();
   if(nClusTot>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
   else                           fIsWEIGHT = kFALSE;
 
   InitVars(nClusTot);
   
-  Float_t xRa,yRa,th,ph;
-  pTrk->GetHMPIDtrk(xRa,yRa,th,ph);        //initialize this track: th and ph angles at middle of RAD 
+  Float_t xPc,yPc,th,ph;
+  pTrk->GetHMPIDtrk(xPc,yPc,th,ph);        //initialize this track: th and ph angles at middle of RAD 
   SetTrack(xRa,yRa,th,ph);
 
   fParam->SetRefIdx(nmean);
 
-  Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;                                                                           
-  fPhotCnt=0;                                                      
+  Float_t mipX=-1,mipY=-1;
+  Int_t chId=-1,mipQ=-1,sizeClu = -1;
+  
+  fPhotCnt=0;
+  
   for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
+    if(iClu == index) {                                                                       // this is the MIP! not a photon candidate: just store mip info
+      mipX = pClu->X();
+      mipY = pClu->Y();
+      mipQ=(Int_t)pClu->Q();
+      sizeClu=pClu->Size();
+      continue;                                                             
+    }
     chId=pClu->Ch();
-    if(pClu->Q()>qthre){                                                                      //charge compartible with MIP clusters      
-      Float_t dX=fPc.X()-pClu->X(),dY=fPc.Y()-pClu->Y(),d =TMath::Sqrt(dX*dX+dY*dY);          //distance between current cluster and intersection point
-      if( d < dMin) {mipId=iClu; dMin=d;mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();} //current cluster is closer, overwrite data for min cluster
-    }else{                                                                                    //charge compatible with photon cluster
-      Double_t thetaCer,phiCer;
-      if(FindPhotCkov(pClu->X(),pClu->Y(),thetaCer,phiCer)){                                  //find ckov angle for this  photon candidate
-        fPhotCkov[fPhotCnt]=thetaCer;                                                         //actual theta Cerenkov (in TRS)
-        fPhotPhi [fPhotCnt]=phiCer;                                                           //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
-       //PH        Printf("photon n. %i reconstructed theta = %f",fPhotCnt,fPhotCkov[fPhotCnt]);
-        fPhotCnt++;                                                                           //increment counter of photon candidates
-      }
+    Double_t thetaCer,phiCer;
+    if(FindPhotCkov(pClu->X(),pClu->Y(),thetaCer,phiCer)){                                    //find ckov angle for this  photon candidate
+      fPhotCkov[fPhotCnt]=thetaCer;                                                           //actual theta Cerenkov (in TRS)
+      fPhotPhi [fPhotCnt]=phiCer;                                                             //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
+      fPhotCnt++;                                                                             //increment counter of photon candidates
     }
   }//clusters loop
+
+  pTrk->SetHMPIDmip(mipX,mipY,mipQ,fPhotCnt);                                                 //store mip info in any case 
+  pTrk->SetHMPIDcluIdx(chId,index+1000*sizeClu);                                              //set index of cluster
   
-  if(fPhotCnt<=nMinPhotAcc) {                                                                 //no reconstruction with <=3 photon candidates
+  if(fPhotCnt<nMinPhotAcc) {                                                                  //no reconstruction with <=3 photon candidates
     pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //set the appropriate flag
-    pTrk->SetHMPIDmip(mipX,mipY,mipQ,fPhotCnt);                                               //store mip info 
-    pTrk->SetHMPIDcluIdx(-1,-1);                                                              //set index of cluster
     return;
   }
   
-  if(mipId==-1) {
-    pTrk->SetHMPIDcluIdx(chId,9999);                                                          //set index of cluster
-    pTrk->SetHMPIDsignal(kMipQdcCut);
-    return;
-  }                                                                                           //no clusters with QDC more the threshold at all
-    pTrk->SetHMPIDcluIdx(chId,mipId);                                                         //set index of cluster
-    if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;}                   //closest cluster with enough charge is still too far from intersection
   
   fMipPos.Set(mipX,mipY);
   
@@ -144,7 +141,7 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t n
   Int_t iNrec=FlagPhot(HoughResponse());                                                      //flag photons according to individual theta ckov with respect to most probable
   pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                    //store mip info 
 
-  if(iNrec<1){
+  if(iNrec<nMinPhotAcc){
     pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
     return;
   }
@@ -477,3 +474,31 @@ Double_t AliHMPIDRecon::HoughResponse()
   return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov   
 }//HoughResponse()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  Double_t AliHMPIDRecon::FindRingExt(Double_t ckov,Int_t ch,Double_t xPc,Double_t yPc,Double_t thRa,Double_t phRa)
+{
+// To find the acceptance of the ring even from external inputs. 
+//    
+//       
+  Double_t xRa = xPc - (fParam->RadThick()+fParam->WinThick()+fParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
+  Double_t yRa = yPc - (fParam->RadThick()+fParam->WinThick()+fParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);
+  
+  Int_t nStep = 500;
+  Int_t nPhi = 0;  
+
+  Int_t ipc,ipadx,ipady;
+    
+  if(ckov>0){
+    SetTrack(xRa,yRa,thRa,phRa);
+    for(Int_t j=0;j<nStep;j++){
+      TVector2 pos; pos=TracePhot(ckov,j*TMath::TwoPi()/(Double_t)(nStep-1));
+      if(fParam->IsInDead(pos.X(),pos.Y())) continue;
+      fParam->Lors2Pad(pos.X(),pos.Y(),ipc,ipadx,ipady);
+      ipadx+=(ipc%2)*fParam->kPadPcX;
+      ipady+=(ipc/2)*fParam->kPadPcY;
+      if(fParam->IsDeadPad(ipadx,ipady,ch)) continue;
+      nPhi++;
+    }//point loop
+  return ((Double_t)nPhi/(Double_t)nStep); 
+  }//if
+  return -1;
+}