Improved hadron correction (Magali, Mark)

[u/mrichter/AliRoot.git] / JETAN / AliDAJetFinder.cxx

diff --git a/JETAN/AliDAJetFinder.cxx b/JETAN/AliDAJetFinder.cxx
index c3d2cf7cdd999678ac6d3aac169d3ab6d43a6e89..130dc076180f43077e0b290190df8e195ab5c90d 100644 (file)
--- a/JETAN/AliDAJetFinder.cxx
+++ b/JETAN/AliDAJetFinder.cxx
@@ -96,23 +96,28 @@ void AliDAJetFinder::InitDetAnn(Double_t &dEtSum,Double_t **xData,TVectorD *vPx,
 {
 //Initialise the variables used by the algorithm
        fBeta=0.1;
-       TClonesArray *lvArray = fReader->GetMomentumArray();
-       fNin = lvArray->GetEntries();
        fNclustMax= ((AliDAJetHeader*)fHeader)->GetFixedCl() ? 
            ((AliDAJetHeader*)fHeader)->GetNclustMax() : 
            TMath::Max((Int_t)TMath::Sqrt(fNin),5);
+       TClonesArray *lvArray = fReader->GetMomentumArray();
+       Int_t nEntr = lvArray->GetEntries();
+       fNin=0;
+       for (Int_t iEn=0; iEn<nEntr; iEn++) if (fReader->GetCutFlag(iEn)==1) fNin++;
        Double_t *xEta = new Double_t[fNin];
-    Double_t *xPhi = new Double_t[fNin];
+       Double_t *xPhi = new Double_t[fNin];
        xData[0]=xEta; xData[1]=xPhi;
        vPx->ResizeTo(fNin);
-       for (Int_t iIn=0; iIn<fNin; iIn++){
-               TLorentzVector *lv=(TLorentzVector*)lvArray->At(iIn);
+       Int_t iIn=0;
+       for (Int_t iEn=0; iEn<nEntr; iEn++){
+               if (fReader->GetCutFlag(iEn)==0) continue;
+               TLorentzVector *lv=(TLorentzVector*)lvArray->At(iEn);
                xEta[iIn] = lv->Eta();
                xPhi[iIn] = lv->Phi()<0 ? lv->Phi() + 2*TMath::Pi() : lv->Phi();
                (*vPx)(iIn)=lv->Pt();
                dEtSum+=(*vPx)(iIn);
+               iIn++;
        }
-       for (Int_t iIn=0; iIn<fNin; iIn++) (*vPx)(iIn)=(*vPx)(iIn)/dEtSum;
+       for (iIn=0; iIn<fNin; iIn++) (*vPx)(iIn)=(*vPx)(iIn)/dEtSum;
 
        Int_t njdim=2*fNclustMax+1;
        mPyx->ResizeTo(fNin,njdim);
@@ -122,17 +127,16 @@ void AliDAJetFinder::InitDetAnn(Double_t &dEtSum,Double_t **xData,TVectorD *vPx,
        (*vPy)(0)=1;
        TMatrixDColumn(*mPyx,0)=1;
        Double_t ypos=0,xpos=0;
-       for (Int_t iIn=0; iIn<fNin; iIn++){
+       for (iIn=0; iIn<fNin; iIn++){
                (*mY)(0,0)+=(*vPx)(iIn)*xEta[iIn];
                ypos+=(*vPx)(iIn)*TMath::Sin(xPhi[iIn]);
                xpos+=(*vPx)(iIn)*TMath::Cos(xPhi[iIn]);
        }
        (*mY)(1,0)=(atan2(ypos,xpos)>0) ? atan2(ypos,xpos) : atan2(ypos,xpos)+2*TMath::Pi();
-       lvArray->Delete();
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk,TVectorD *vPy,TMatrixD *mY)
+void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk,  TVectorD *vPy,  TMatrixD *mY) const
 {
        for(Int_t iClust=0; iClust<nc; iClust++){
                (*vPy)(iClust)=(*vPy)(iClust)/2;
@@ -143,7 +147,7 @@ void AliDAJetFinder::DoubleClusters(Int_t nc,Int_t &nk,TVectorD *vPy,TMatrixD *m
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::Annealing(Int_t nk,Double_t **xData,TVectorD *vPx,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
+void AliDAJetFinder::Annealing(Int_t nk,Double_t **xData,  TVectorD *vPx,  TVectorD *vPy,  TMatrixD *mPyx,  TMatrixD *mY)
 {
 // Main part of the algorithm
        const Double_t pi=TMath::Pi();
@@ -234,11 +238,11 @@ void AliDAJetFinder::Annealing(Int_t nk,Double_t **xData,TVectorD *vPx,TVectorD
     delete y;
     delete y1;
     delete ry;
-       delete [] m;
+    delete [] m;
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy,TMatrixD *mPyx,TMatrixD *mY)
+void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy,  TMatrixD *mPyx,TMatrixD *mY)
 {
        static Bool_t growcl=true;
        
@@ -295,7 +299,7 @@ void AliDAJetFinder::NumCl(Int_t &nc,Int_t &nk,TVectorD *vPy,TMatrixD *mPyx,TMat
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut)
+void AliDAJetFinder::ReduceClusters(Int_t **iSame,Int_t nc,Int_t &ncout,Int_t **cont,Int_t *nSameOut) const
 {
        Int_t *nSame = new Int_t[nc];
        Int_t *iperm = new Int_t[nc];
@@ -374,7 +378,7 @@ void AliDAJetFinder::EndDetAnn(Int_t &nk,Double_t **xData,Int_t *xx,Double_t etx
 }
 
 //-----------------------------------------------------------------------------------
-void AliDAJetFinder::StoreJets(Int_t nk,Double_t **xData,Int_t *xx,TMatrixD *mY)
+void AliDAJetFinder::StoreJets(Int_t nk,Double_t **xData,  Int_t *xx,  TMatrixD *mY)
 {
 //evaluate significant clusters properties
        const Double_t pi=TMath::Pi();
@@ -463,8 +467,15 @@ void AliDAJetFinder::StoreJets(Int_t nk,Double_t **xData,Int_t *xx,TMatrixD *mY)
                        pz = (*mY)(3,iClust)/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-(*mY)(0,iClust))));
                        en = TMath::Sqrt(px * px + py * py + pz * pz);
                        AliAODJet jet(px, py, pz, en);
-                       if (fromAod) 
-                           for (Int_t iIn=0; iIn<fNin; iIn++) if (xx[iIn]==iClust) jet.AddTrack(refs->At(iIn));
+                       if (fromAod){
+                               Int_t iIn=0;
+                               Int_t nEntr = fReader->GetMomentumArray()->GetEntries();
+                               for (Int_t iEn=0; iEn<nEntr; iEn++){
+                                       if (fReader->GetCutFlag(iEn)==0) continue;
+                                       if (xx[iIn]==iClust) jet.AddTrack(refs->At(iEn));
+                                       iIn++;
+                               }
+                       }
                        AddJet(jet);
                        printf("jet %d, Eta: %f, Phi: %f, Et: %f\n",iClust,jet.Eta(),jet.Phi(),jet.Pt());
                }