]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HMPID/AliHMPIDRecon.cxx
add calculation and histograms for MC cross section
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
index 85495c40d6f1247ebb352c3348425e47fe8ef583..792da47cbc64c909f162de2bed36735a08692561 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),
@@ -62,6 +64,7 @@ void AliHMPIDRecon::InitVars(Int_t n)
 //..
   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];
@@ -74,6 +77,7 @@ void AliHMPIDRecon::DeleteVars()const
 //Delete variables
 //..
   delete [] fPhotFlag;
+  delete [] fPhotClusIndex;
   delete [] fPhotCkov;
   delete [] fPhotPhi;
   delete [] fPhotWei;
@@ -85,13 +89,10 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
 // 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);
   
@@ -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
@@ -495,6 +519,7 @@ Double_t AliHMPIDRecon::HoughResponse()
       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