Changed coordinate system
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderAlgoOmni.cxx
index 9170a8ef406b8e8a8a1f4f8fc96c9c8b0f9ae294..763d5759ea0c7cc3ac93ea67f334b341e09f966f 100644 (file)
 /*
  
 $Log$
+Revision 1.4  2003/09/19 13:16:20  mhorner
+Added additional jet energy info
 
 
+Revision 1.3  2003/09/04 12:49:56  mhorner
+Changed hadron correction and added saving EMCAL and track contributions
 
 */
 
@@ -37,7 +41,6 @@ $Log$
 #include "AliEMCALJetFinderAlgoUA1Unit.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCAL.h"
-#include "AliEMCALGetter.h"
 #include "AliEMCALDigit.h"
 #include "TParticle.h"
 #include "AliRun.h"
@@ -57,7 +60,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
   fESeed             = 5.0;       //Default value
   fConeRad           = 0.3;       //Default value
   fJetEMin           = 10.0;      //Default value
-  fEtMin             = 0.28;      //Default value
+  fEtMin             = 0.0;      //Default value
   fMinMove           = 0.05;      //From original UA1 JetFinder
   fMaxMove           = 0.15;      //From original UA1 JetFinder
   fBGMaxMove         = 0.035;     //From original UA1 JetFinder
@@ -142,15 +145,14 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
  void AliEMCALJetFinderAlgoOmni::FillUnitArray(AliEMCALJetFinderAlgoUA1FillUnitFlagType_t flag)
    {
      if (fDebug>1) Info("FillUnitArray","Beginning FillUnitArray");
-     //     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 
          //   if (pEMCAL){ 
          //         AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
          //     }else
          //    {
-     //AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
-     AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry() ; 
-       //    }
+     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
+        //    }
          
      AliEMCALJetFinderAlgoUA1FillUnitFlagType_t option = flag;
      Int_t         numTracks, numDigits;
@@ -177,17 +179,20 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
             Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy(); 
             Float_t unitEnergyNoCuts = fUnitNoCuts[towerID-1].GetUnitEnergy();
 
-            /*
-            //OLD WAY:              //Do Hadron Correction
+            
+            //OLD WAY:   //Do Hadron Correction
               if(fHadCorr != 0)
               {
                 Double_t   fullP = myPart->P();
                 Double_t   hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
                 unitEnergy -= hCEnergy*TMath::Sin(myPart->Theta());
+                unitEnergyNoCuts -= hCEnergy*TMath::Sin(myPart->Theta());
                 fUnit[towerID-1].SetUnitEnergy(unitEnergy);
+                fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts);
               } //end Hadron Correction loop 
-            */      
-            
+             
+    
+            /*
              //Do Hadron Correction with propagate phi for the track
              if(fHadCorr != 0)
                {
@@ -204,7 +209,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
                  phi += deltaPhi;
                  //Get new tower id for cell that track would curve into
                   Int_t towerID2;
-                  if(phi<(1.0/3.0)*TMath::Pi() || phi>TMath::Pi())
+                  if(phi>(TMath::Pi()/180.0)*geom->GetArm1PhiMax() || phi<(TMath::Pi()/180.0)*geom->GetArm1PhiMin())
                     {
                       towerID2 = -1;
                     }
@@ -225,6 +230,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
                      fUnitNoCuts[towerID2-1].SetUnitEnergy(unitEnergy2NoCuts);
                    }//end if for towerID2
                }//end Hadron Correction loop
+            */
 
              fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts + pT);
             //Do Pt cut on tracks
@@ -258,15 +264,16 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
         for(Int_t i=0; i<fNumUnits; i++)
           {
              if (fDebug>10) Info("FillUnitArray","Setting all units outside jets");
-            fUnit[i].SetUnitFlag(kOutJet);           //Set all units to be outside a jet initially
+            //Set all units to be outside a jet initially
+            fUnit[i].SetUnitFlag(kOutJet);           
             fUnit[i].SetUnitID(i+1);
             Float_t eta;
             Float_t phi;
             geom->EtaPhiFromIndex(fUnit[i].GetUnitID(), eta, phi);
             fUnit[i].SetUnitEta(eta);
             fUnit[i].SetUnitPhi(phi*TMath::Pi()/180.0);
-
-            fUnitNoCuts[i].SetUnitFlag(kOutJet);           //Set all units to be outside a jet initially
+            //Set all units to be outside a jet initially
+            fUnitNoCuts[i].SetUnitFlag(kOutJet);          
             fUnitNoCuts[i].SetUnitID(i+1);
             eta = 0.0;
             phi = 0.0;
@@ -464,7 +471,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
    {
      //Stores the resulting jet information in appropriate storage structure (TO BE DECIDED!!!!)
      if (fDebug>1) Info("StoreJetInfo","Storing Jet Information");
-    
+     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance("EMCAL_5655_21", "");
      //Store:
      //fJetESum is the final jet energy (background has been subtracted)
      //fJetEta is the final jet Eta
@@ -482,8 +489,14 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
       cout<<"The jet phi is ---->" <<fJetPhi <<endl;
 
      Int_t             numberTracks = fInputPointer->GetNTracks();
+     Int_t            numberDigits = fInputPointer->GetNDigits();
+     AliEMCALDigit     *myD;
      TParticle         *myP;
      Int_t             numTracksInCone = 0;
+     Float_t           trackEnergy = 0.0;
+     Float_t           trackEnergyPtCut =0.0; 
+     Float_t           emcalEnergy = 0.0;
+     Float_t           emcalEnergyBGSub = 0.0;
 
      for(Int_t counter=0; counter<numberTracks; counter++)
        {
@@ -513,14 +526,52 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor");
         if(rad<=fConeRad)
           {
             pTArray[index] = myP->Pt();
+            //Calculate track contribution within jetcone
+            trackEnergy += myP->Pt();
+            if(myP->Pt() >= fPtCut) trackEnergyPtCut += myP->Pt();
             etaArray[index] = eta;
             phiArray[index] = phi;
             pdgArray[index] = myP->GetPdgCode();
             index++;
+            if(fHadCorr != 0)
+            {
+                    Double_t   fullP = myP->P();
+                    Double_t   hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta);
+                    emcalEnergy -= hCEnergy*TMath::Sin(myP->Theta());
+                    emcalEnergyBGSub -= hCEnergy*TMath::Sin(myP->Theta());
+            } //end Hadron Correction loop 
+                          
           }//end if
        }//end for
 
+     //Loop over digits to find EMCal contribution within jetcone
+     for(Int_t counter3=0; counter3<numberDigits; counter3++)
+       {
+        myD = fInputPointer->GetDigit(counter3);
+        //GET DIGIT ETA, PHI so that can check if inside R!
+        Float_t eta = 0.0;
+        Float_t phi = 0.0;
+        Int_t ID = myD->GetId();
+        geom->EtaPhiFromIndex(ID, eta, phi);
+        Float_t deta = fJetEta-eta;
+        Float_t dphi = fJetPhi -(TMath::Pi()/180.0)*phi;
+        Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi));
+        if(rad<=fConeRad)
+          {
+        Int_t amplitude = myD->GetAmp();     //Gets the integer valued amplitude of the digit
+        Float_t amp = (Float_t)amplitude;        //Need to typecast to Float_t before doing real energy conversion
+        Float_t digitEnergy = amp/10000000.0;    //Factor of 10 million needed to convert!
+        emcalEnergy += digitEnergy;
+        emcalEnergyBGSub += (digitEnergy - fEBGAve);
+          }//end if
+       }//end count3 for
+
+     //Save in JET object
      fJet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray);
+     fJet.SetEMCALEnergy(emcalEnergy);
+     fJet.SetEMCALEnergyBGSub(emcalEnergyBGSub);
+     fJet.SetTrackEnergy(trackEnergy);
+     fJet.SetTrackEnergyPtCut(trackEnergyPtCut);
      fOutputObject.AddJet(&fJet);
      delete[] pTArray;
      delete[] etaArray;