Adding more bins in QA (Alis)
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
index e8cb828..792da47 100644 (file)
 #include <TH1D.h>            //HoughResponse()
 #include <TClonesArray.h>    //CkovAngle()
 #include <AliESDtrack.h>     //CkovAngle()
+#include <AliESDfriendTrack.h>     //CkovAngle()
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDRecon::AliHMPIDRecon():
-  TTask("RichRec","RichPat"),
+  TNamed("RichRec","RichPat"),
   fPhotCnt(-1),
   fPhotFlag(0x0),
+  fPhotClusIndex(0x0),    
   fPhotCkov(0x0),
   fPhotPhi(0x0),
   fPhotWei(0x0),
@@ -60,8 +62,9 @@ void AliHMPIDRecon::InitVars(Int_t n)
 //..
 //Init some variables
 //..
-  if(n<0) return;
+  if(n<=0) return;
   fPhotFlag = new Int_t[n];
+  fPhotClusIndex  = new Int_t[n];
   fPhotCkov = new Double_t[n];
   fPhotPhi  = new Double_t[n];
   fPhotWei  = new Double_t[n];
@@ -73,30 +76,28 @@ void AliHMPIDRecon::DeleteVars()const
 //..
 //Delete variables
 //..
-  delete fPhotFlag;
-  delete fPhotCkov;
-  delete fPhotPhi;
-  delete fPhotWei;
+  delete [] fPhotFlag;
+  delete [] fPhotClusIndex;
+  delete [] fPhotCkov;
+  delete [] fPhotPhi;
+  delete [] fPhotWei;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean)
+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
 //              pCluLst  - list of clusters for this chamber   
 //   Returns:            - track ckov angle, [rad], 
-    
-
+       
   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);
@@ -106,8 +107,11 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
   
   fPhotCnt=0;
   
+  Int_t nPads = 0;
+  
   for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
-    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
+    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster
+    nPads+=pClu->Size();    
     if(iClu == index) {                                                                       // this is the MIP! not a photon candidate: just store mip info
       mipX = pClu->X();
       mipY = pClu->Y();
@@ -116,10 +120,12 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
       continue;                                                             
     }
     chId=pClu->Ch();
+    if(pClu->Q()>2*fParam->QCut()) continue;
     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)
+      fPhotPhi [fPhotCnt]=phiCer;
+      fPhotClusIndex[fPhotCnt]=iClu;                                                             //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
       fPhotCnt++;                                                                             //increment counter of photon candidates
     }
   }//clusters loop
@@ -127,29 +133,33 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
   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
     return;
   }
-  
-  
+    
   fMipPos.Set(mipX,mipY);
-  
-  
+    
 //PATTERN RECOGNITION STARTED: 
+  if(fPhotCnt>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
+  else                           fIsWEIGHT = kFALSE;
+    
+  Int_t iNrec=FlagPhot(HoughResponse(),pCluLst,pTrk);                                                      //flag photons according to individual theta ckov with respect to most probable
   
-  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 
+  pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                  //store mip info 
 
-  if(iNrec<1){
-    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
+  if(iNrec<nMinPhotAcc){
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                    //no photon candidates are accepted
     return;
   }
+  
+  Int_t occupancy = (Int_t)(1000*(nPads/(6.*80.*48.))); 
+  
   Double_t thetaC = FindRingCkov(pCluLst->GetEntries());                                    //find the best reconstructed theta Cherenkov
 //    FindRingGeom(thetaC,2);
-  pTrk->SetHMPIDsignal(thetaC);                                                             //store theta Cherenkov
-  pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //store errors squared
-
+  pTrk->SetHMPIDsignal(thetaC+occupancy);                                                   //store theta Cherenkov and chmaber occupancy
+  pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //store experimental ring angular resolution squared
+  
   DeleteVars();
 }//CkovAngle()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -353,7 +363,7 @@ Double_t AliHMPIDRecon::FindRingCkov(Int_t)
   return weightThetaCerenkov;
 }//FindCkovRing()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
+Int_t AliHMPIDRecon::FlagPhot(Double_t ckov,TClonesArray *pCluLst, AliESDtrack *pTrk)
 {
 // Flag photon candidates if their individual ckov angle is inside the window around ckov angle returned by  HoughResponse()
 // Arguments: ckov- value of most probable ckov angle for track as returned by HoughResponse()
@@ -374,14 +384,27 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
   Int_t iInsideCnt = 0; //count photons which Theta ckov inside the window
   for(Int_t i=0;i<fPhotCnt;i++){//photon candidates loop
     fPhotFlag[i] = 0;
-    if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax)   { 
-      fPhotFlag[i]=2;    
+    if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) { 
+      fPhotFlag[i]=2;
+      AddObjectToFriends(pCluLst,i,pTrk);
       iInsideCnt++;
     }
   }
+      
   return iInsideCnt;
+  
 }//FlagPhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void  AliHMPIDRecon::AddObjectToFriends(TClonesArray *pCluLst, Int_t photonIndex, AliESDtrack *pTrk)
+{
+// Add AliHMPIDcluster object to ESD friends
+    
+  AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(fPhotClusIndex[photonIndex]);     
+  AliHMPIDCluster *pClus = new AliHMPIDCluster(*pClu);
+  pClus->SetChi2(fPhotCkov[photonIndex]);  
+  pTrk->AddCalibObject(pClus);   
+}    
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const
 {
 // Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses
@@ -464,6 +487,7 @@ Double_t AliHMPIDRecon::HoughResponse()
     Double_t sumPhots=phots->Integral(bin1,bin2);
     if(sumPhots<3) continue;                            // if less then 3 photons don't trust to this ring
     Double_t sumPhotsw=photsw->Integral(bin1,bin2);
+    if((Double_t)((i+0.5)*fDTheta)>0.7) continue;
     resultw->Fill((Double_t)((i+0.5)*fDTheta),sumPhotsw);
   } 
 // evaluate the "BEST" theta ckov as the maximum value of histogramm
@@ -474,3 +498,32 @@ 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(ipadx<0 || ipady>160 || ipady<0 || ipady>144 || ch<0 || ch>6) continue;
+      if(fParam->IsDeadPad(ipadx,ipady,ch)) continue;
+      nPhi++;
+    }//point loop
+  return ((Double_t)nPhi/(Double_t)nStep); 
+  }//if
+  return -1;
+}