From 5637918c0f5b51643ee1f5bf00cb8ebe92e73d43 Mon Sep 17 00:00:00 2001
From: hristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Sun, 24 Oct 2010 15:34:58 +0000
Subject: [PATCH] HMPID friends (Giacomo)

---
 HMPID/AliHMPIDCluster.h     |  1 +
 HMPID/AliHMPIDRecon.cxx     | 46 +++++++++++++++++++++++++++++--------
 HMPID/AliHMPIDRecon.h       |  5 ++--
 STEER/AliESDfriendTrack.cxx | 21 +++++++++++++++++
 STEER/AliESDfriendTrack.h   |  8 ++++++-
 test/gun/rec.C              |  2 ++
 6 files changed, 70 insertions(+), 13 deletions(-)

diff --git a/HMPID/AliHMPIDCluster.h b/HMPID/AliHMPIDCluster.h
index c7a7914500a..337a89ed131 100644
--- a/HMPID/AliHMPIDCluster.h
+++ b/HMPID/AliHMPIDCluster.h
@@ -63,6 +63,7 @@ public:
          void           SetY     (Double_t y                               ){fYY=y;}                                               // Setter
          void           SetSize  (Int_t size                               ){fSi=size;}                                            // Setter
          void           FindClusterSize(Int_t i,Int_t *pSigmaCut);                                                                 //Find the clusterSize of deconvoluted clusters 
+ virtual void	        Clear(const Option_t*) { delete [] fDigs; fDigs=0; delete [] fParam; fParam=0; }
          
 protected:
   Int_t         fCh;          //chamber number
diff --git a/HMPID/AliHMPIDRecon.cxx b/HMPID/AliHMPIDRecon.cxx
index 27069b2a598..9195dd67c43 100644
--- a/HMPID/AliHMPIDRecon.cxx
+++ b/HMPID/AliHMPIDRecon.cxx
@@ -27,12 +27,14 @@
 #include <TH1D.h>            //HoughResponse()
 #include <TClonesArray.h>    //CkovAngle()
 #include <AliESDtrack.h>     //CkovAngle()
+#include <AliESDfriendTrack.h>     //CkovAngle()
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDRecon::AliHMPIDRecon():
   TTask("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;
@@ -86,7 +90,9 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
 //              pCluLst  - list of clusters for this chamber   
 //   Returns:            - track ckov angle, [rad], 
     
-
+   
+  AliESDfriendTrack *pFriendTrk = (AliESDfriendTrack*)pTrk->GetFriendTrack();
+  
   const Int_t nMinPhotAcc = 3;                      // Minimum number of photons required to perform the pattern recognition
   
   Int_t nClusTot = pCluLst->GetEntries();
@@ -119,7 +125,8 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
     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 +134,31 @@ 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);
-  
+    
+  TClonesArray *pPhotCluLst = pFriendTrk->GetHmpPhotClus();
   
 //PATTERN RECOGNITION STARTED: 
   
-  Int_t iNrec=FlagPhot(HoughResponse());                                                      //flag photons according to individual theta ckov with respect to most probable
+  Int_t iNrec=FlagPhot(HoughResponse(),pCluLst,pPhotCluLst);                                                      //flag photons according to individual theta ckov with respect to most probable
+  
   pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                    //store mip info 
 
-  if(iNrec<nMinPhotAcc){
+  if(iNrec<1){
     pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
     return;
   }
+  
   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
-
+  
   DeleteVars();
 }//CkovAngle()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -353,7 +362,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, TClonesArray *pPhotCluLst)
 {
 // 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()
@@ -362,6 +371,7 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
 // Photon Flag:  Flag = 0 initial set; 
 //               Flag = 1 good candidate (charge compatible with photon); 
 //               Flag = 2 photon used for the ring;
+  Int_t *PhotIndex = new Int_t[fPhotCnt];
   
   Int_t steps = (Int_t)((ckov )/ fDTheta); //how many times we need to have fDTheta to fill the distance between 0  and thetaCkovHough
 
@@ -376,10 +386,26 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
     fPhotFlag[i] = 0;
     if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax)	{ 
       fPhotFlag[i]=2;	  
+      PhotIndex[iInsideCnt]=fPhotClusIndex[i];
       iInsideCnt++;
     }
   }
+      
+  Int_t nPhot = 0;
+     
+  for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
+    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster
+    for(Int_t j=0; j<iInsideCnt; j++){
+      if(iClu==PhotIndex[j]) {
+      new ((*pPhotCluLst)[nPhot++]) AliHMPIDCluster(*pClu);  //add this raw cluster 
+     } 
+    }
+  } 
+                                                                                      
+  delete [] PhotIndex;
+  
   return iInsideCnt;
+  
 }//FlagPhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const
diff --git a/HMPID/AliHMPIDRecon.h b/HMPID/AliHMPIDRecon.h
index 8b1f8dc1268..b0dcdc921ba 100644
--- a/HMPID/AliHMPIDRecon.h
+++ b/HMPID/AliHMPIDRecon.h
@@ -33,7 +33,7 @@ public :
   Double_t FindRingCkov (Int_t iNclus                                                       );     //best ckov for ring formed by found photon candidates
   void     FindRingGeom (Double_t ckovAng,Int_t level=1                                     );     //estimated area of ring in cm^2 and portion accepted by geometry
   TVector2 IntWithEdge  (TVector2 p1,TVector2 p2                                            )const;//find intercection between plane and lines of 2 thetaC
-  Int_t    FlagPhot     (Double_t ckov                                                      );     //is photon ckov near most probable track ckov
+  Int_t    FlagPhot     (Double_t ckov,TClonesArray *pCluLst,TClonesArray *pPhotCluLst      );     //is photon ckov near most probable track ckov
   Double_t HoughResponse(                                                                   );     //most probable track ckov angle
   void     Propagate    (const TVector3  dir,      TVector3 &pos,Double_t z                 )const;//propagate photon alogn the line  
   void     Refract      (      TVector3 &dir,                    Double_t n1,    Double_t n2)const;//refract photon on the boundary
@@ -59,6 +59,7 @@ public :
 protected:
   Int_t     fPhotCnt;                           // counter of photons candidate
   Int_t    *fPhotFlag;                          // flags of photon candidates
+  Int_t    *fPhotClusIndex;                     // cluster index of photon candidates
   Double_t *fPhotCkov;                          // Ckov angles of photon candidates, [rad]
   Double_t *fPhotPhi;                           // phis of photons candidates, [rad]
   Double_t *fPhotWei;                           // weigths of photon candidates
@@ -81,7 +82,7 @@ private:
   AliHMPIDRecon(const AliHMPIDRecon& r);              //dummy copy constructor
   AliHMPIDRecon &operator=(const AliHMPIDRecon& r);   //dummy assignment operator
 //
-  ClassDef(AliHMPIDRecon,1)
+  ClassDef(AliHMPIDRecon,2)
 };
 
 #endif // #ifdef AliHMPIDRecon_cxx
diff --git a/STEER/AliESDfriendTrack.cxx b/STEER/AliESDfriendTrack.cxx
index 508ffe2d5a4..f5b0ad6781d 100644
--- a/STEER/AliESDfriendTrack.cxx
+++ b/STEER/AliESDfriendTrack.cxx
@@ -21,6 +21,7 @@
 #include "AliTrackPointArray.h"
 #include "AliESDfriendTrack.h"
 #include "TObjArray.h"
+#include "TClonesArray.h"
 #include "AliKalmanTrack.h"
 
 ClassImp(AliESDfriendTrack)
@@ -28,6 +29,7 @@ ClassImp(AliESDfriendTrack)
 AliESDfriendTrack::AliESDfriendTrack(): 
 TObject(), 
 f1P(0), 
+fHmpPhotClus(new TClonesArray("AliHMPIDCluster",1000)),
 fnMaxITScluster(0),
 fnMaxTPCcluster(0),
 fnMaxTRDcluster(0),
@@ -52,11 +54,15 @@ fTRDIn(0)
   //for (i=0; i<kMaxITScluster; i++) fITSindex[i]=-2;
   //for (i=0; i<kMaxTPCcluster; i++) fTPCindex[i]=-2;
   //for (i=0; i<kMaxTRDcluster; i++) fTRDindex[i]=-2;
+  
+  //fHmpPhotClus->SetOwner(kTRUE); 
+  
 }
 
 AliESDfriendTrack::AliESDfriendTrack(const AliESDfriendTrack &t): 
 TObject(t),
 f1P(t.f1P),
+fHmpPhotClus(t.fHmpPhotClus),        
 fnMaxITScluster(t.fnMaxITScluster),
 fnMaxTPCcluster(t.fnMaxTPCcluster),
 fnMaxTRDcluster(t.fnMaxTRDcluster),
@@ -106,6 +112,8 @@ fTRDIn(0)
   if (t.fITSOut) fITSOut = new AliExternalTrackParam(*(t.fITSOut));
   if (t.fTRDIn)  fTRDIn = new AliExternalTrackParam(*(t.fTRDIn));
 
+  fHmpPhotClus = new TClonesArray(*t.fHmpPhotClus);  
+  
 }
 
 AliESDfriendTrack::~AliESDfriendTrack() {
@@ -115,6 +123,7 @@ AliESDfriendTrack::~AliESDfriendTrack() {
    delete fPoints;
    if (fCalibContainer) fCalibContainer->Delete();
    delete fCalibContainer;
+   delete fHmpPhotClus;   
    delete fITStrack;
    delete fTRDtrack;
    delete fTPCOut;
@@ -167,6 +176,18 @@ void AliESDfriendTrack::SetTRDIn(const AliExternalTrackParam  &param)  {
   fTRDIn=new AliExternalTrackParam(param);
 } 
 
+
+void AliESDfriendTrack::SetHmpPhotClus(TClonesArray  *array)  {
+  //
+  // 
+  //
+  fHmpPhotClus->Clear("C");
+  for(Int_t iClu=0; iClu<array->GetEntriesFast(); iClu++){
+     TObject *pClu = (TObject*)array->UncheckedAt(iClu);
+     new ((*fHmpPhotClus)[iClu]) TObject(*pClu);
+  }
+}
+
 void AliESDfriendTrack::SetITSIndices(Int_t* indices, Int_t n){
 
 	//
diff --git a/STEER/AliESDfriendTrack.h b/STEER/AliESDfriendTrack.h
index 7e2acc6d9d5..3edfa97706a 100644
--- a/STEER/AliESDfriendTrack.h
+++ b/STEER/AliESDfriendTrack.h
@@ -8,6 +8,7 @@
 //-------------------------------------------------------------------------
 
 #include <TObject.h>
+#include <TClonesArray.h>
 #include <AliExternalTrackParam.h>
 
 class AliTrackPointArray;
@@ -47,6 +48,8 @@ public:
   void SetITSOut(const AliExternalTrackParam &param);
   void SetTRDIn(const AliExternalTrackParam  &param);
   //
+  void SetHmpPhotClus(TClonesArray *array);
+  
   const AliExternalTrackParam * GetTPCOut() const {return  fTPCOut;} 
   const AliExternalTrackParam * GetITSOut() const {return fITSOut;} 
   const AliExternalTrackParam * GetTRDIn()  const {return fTRDIn;} 
@@ -59,12 +62,15 @@ public:
   Int_t GetMaxTPCcluster() {return fnMaxTPCcluster;}
   Int_t GetMaxTRDcluster() {return fnMaxTRDcluster;}
 
+  TClonesArray *GetHmpPhotClus() const {return fHmpPhotClus;}   
+  
   // bit manipulation for filtering
   void SetSkipBit(Bool_t skip){SetBit(23,skip);}
   Bool_t TestSkipBit() {return TestBit(23);}
 
 protected:
   Float_t f1P;                     // 1/P (1/(GeV/c))
+  TClonesArray *fHmpPhotClus; // TClonesArray of reconstructed photon clusters  
   Int_t fnMaxITScluster; // Max number of ITS clusters
   Int_t fnMaxTPCcluster; // Max number of TPC clusters
   Int_t fnMaxTRDcluster; // Max number of TRD clusters
@@ -85,7 +91,7 @@ protected:
 private:
   AliESDfriendTrack &operator=(const AliESDfriendTrack & /* t */) {return *this;}
 
-  ClassDef(AliESDfriendTrack,4) //ESD friend track
+  ClassDef(AliESDfriendTrack,5) //ESD friend track
 };
 
 #endif
diff --git a/test/gun/rec.C b/test/gun/rec.C
index e5b3ba99a2d..574df50f7eb 100644
--- a/test/gun/rec.C
+++ b/test/gun/rec.C
@@ -9,6 +9,8 @@ void rec() {
 			  Form("local://%s",gSystem->pwd()));
   reco.SetRunPlaneEff(kTRUE);
 
+  reco.SetFractionFriends(1.);
+
   TStopwatch timer;
   timer.Start();
   reco.Run();
-- 
2.43.5