]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALJetFinder.cxx
iostream.h replaced by Riostream.h
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinder.cxx
index 85b815400fa5c59507512fff158aa569b97b12cf..2d8c89de1e758b5105891f4b6d6db46801fb8c5e 100644 (file)
 
 /*
 $Log$
+Revision 1.34  2003/01/14 10:50:18  alibrary
+Cleanup of STEER coding conventions
+
+Revision 1.33  2003/01/10 10:48:19  morsch
+SetSamplingFraction() removed from constructor.
+
+Revision 1.32  2003/01/10 10:26:40  morsch
+Sampling fraction initialized from geometry class.
+
+Revision 1.31  2003/01/08 17:13:41  schutz
+added the HCAL section
+
+Revision 1.30  2002/12/09 16:26:28  morsch
+- Nummber of particles per jet increased to 1000
+- Warning removed.
+
+Revision 1.29  2002/11/21 17:01:40  alibrary
+Removing AliMCProcess and AliMC
+
+Revision 1.28  2002/11/20 14:13:16  morsch
+- FindChargedJets() added.
+- Destructor corrected.
+- Geometry cuts taken from AliEMCALGeometry.
+
 Revision 1.27  2002/11/15 17:39:10  morsch
 GetPythiaParticleName removed.
 
@@ -107,40 +131,39 @@ Revision 1.3  2002/01/18 05:07:56  morsch
 
 #include <stdio.h>
 // From root ...
-#include <TROOT.h>
-#include <TClonesArray.h>
-#include <TTree.h>
+#include <TArrayF.h>
+#include <TAxis.h>
 #include <TBranchElement.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
 #include <TFile.h>
 #include <TH1.h>
 #include <TH2.h>
-#include <TArrayF.h>
-#include <TCanvas.h>
 #include <TList.h>
+#include <TPDGCode.h>
 #include <TPad.h>
-#include <TPaveText.h>
-#include <TAxis.h>
-#include <TStyle.h>
 #include <TParticle.h>
 #include <TParticlePDG.h>
+#include <TPaveText.h>
 #include <TPythia6Calls.h>
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TTree.h>
 
 // From AliRoot ...
-#include "AliEMCALJetFinder.h"
-#include "AliEMCALFast.h"
-#include "AliEMCALGeometry.h"
-#include "AliEMCALHit.h"
+#include "AliEMCAL.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALDigitizer.h"
+#include "AliEMCALFast.h"
+#include "AliEMCALGeometry.h"
 #include "AliEMCALHadronCorrection.h"
+#include "AliEMCALHit.h"
+#include "AliEMCALJetFinder.h"
 #include "AliEMCALJetMicroDst.h"
-#include "AliRun.h"
+#include "AliHeader.h"
 #include "AliMagF.h"
 #include "AliMagFCM.h"
-#include "AliEMCAL.h"
-#include "AliHeader.h"
-#include "AliPDG.h"
-#include "AliMC.h"
+#include "AliRun.h"
 
 // Interface to FORTRAN
 #include "Ecommon.h"
@@ -224,7 +247,6 @@ AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
     SetEfficiencySim();
     SetDebug();
     SetHadronCorrection();
-    SetSamplingFraction();
     SetIncludeK0andN();
 
     SetParametersForBgSubtraction();
@@ -310,7 +332,10 @@ void AliEMCALJetFinder::Init()
 //  Geometry 
     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
     AliEMCALGeometry* geom = 
-       AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+    AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+
+//    SetSamplingFraction(geom->GetSampling());
+
     fNbinEta = geom->GetNZ();
     fNbinPhi = geom->GetNPhi();
     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
@@ -641,7 +666,8 @@ void AliEMCALJetFinder::DumpLego()
             }
        } // phi
     } // eta
-    fNcell--;
+//  today
+//    fNcell--;
 }
 
 void AliEMCALJetFinder::ResetMap()
@@ -736,7 +762,7 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
        fPhiT[part]      = phi;
        fPdgT[part]      = mpart;
        
-
+// today
        if (part < 2) continue;
 
        // move to fLego->Fill because hadron correction must apply 
@@ -752,9 +778,12 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
                if (!fK0N) { 
                    continue;
                } else {
+// today
+/*
                    if (mpart != kNeutron    &&
                        mpart != kNeutronBar &&
                        mpart != kK0Long) continue;
+*/
                }
            }
        }
@@ -792,8 +821,7 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 // Tracking Efficiency and TPC acceptance goes here ...
        Float_t eff;
        if (fEffic && TMath::Abs(chTmp)) {
-         //        eff =  AliEMCALFast::Efficiency(1,p);
-            eff = 0.95; // for testing 25-feb-2002
+           eff =  AliEMCALFast::Efficiency(2,p);
             if(fhEff) fhEff->Fill(p, eff);
            if (AliEMCALFast::RandomReject(eff)) {
               if(fDebug >= 5) printf(" reject due to unefficiency ");
@@ -828,7 +856,7 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 
              if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n", 
              phi, phiHC, -eTdpH); // correction is negative
-             fLego->Fill(eta, phiHC, -eTdpH);
+             fLego->Fill(eta, phiHC, -eTdpH );
               fhLegoHadrCorr->Fill(eta, phiHC, eTdpH);
          }
         }
@@ -856,6 +884,14 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
         }
 
     } // primary loop    
+    Float_t etsum = 0.;
+    for(Int_t i=0; i<fLego->GetSize(); i++) {
+       Float_t etc =  (*fLego)[i];
+       if (etc > fMinCellEt) etsum += etc;
+    }
+    
+    printf("\nFillFromTracks: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
+
     DumpLego();
 }
 
@@ -910,7 +946,7 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag)
            Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
            Float_t phi    =    TMath::ATan2(y,x);
 
-           if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
+           if (fDebug >= 11) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
 //         printf("\n Hit %f %f %f %f", x, y, z, eloss);
            
             etH = fSamplingF*eloss*TMath::Sin(theta);
@@ -919,8 +955,18 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag)
        } // Hit Loop
     } // Track Loop
     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
-    for(Int_t i=0; i<fLego->GetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i]; 
+    Float_t etsum = 0;
+    
+    for(Int_t i=0; i<fLego->GetSize(); i++) {
+       (*fhLegoEMCAL)[i] = (*fLego)[i];
+       Float_t etc =  (*fLego)[i];
+       if (etc > fMinCellEt) etsum += etc;
+    }
+    
+    printf("\nFillFromHits: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
+    
     //    DumpLego(); ??
+    
 }
 
 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
@@ -962,10 +1008,12 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag)
     nbytes += branchDr->GetEntry(0);
 //
 //  Get digitizer parameters
-    Float_t towerADCped  = digr->GetTowerpedestal();
-    Float_t towerADCcha  = digr->GetTowerchannel();
-    Float_t preshoADCped = digr->GetPreShopedestal();
-    Float_t preshoADCcha = digr->GetPreShochannel();
+    Float_t preADCped = digr->GetPREpedestal();
+    Float_t preADCcha = digr->GetPREchannel();
+    Float_t ecADCped  = digr->GetECpedestal();
+    Float_t ecADCcha  = digr->GetECchannel();
+    Float_t hcADCped  = digr->GetHCpedestal();
+    Float_t hcADCcha  = digr->GetHCchannel();
 
     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
     AliEMCALGeometry* geom = 
@@ -973,10 +1021,8 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag)
     
     if (fDebug) {
        Int_t   ndig = digs->GetEntries();
-       printf("\n Number of Digits: %d %d\n", ndig, nent);
-       printf("\n Parameters: %f %f %f %f\n", 
-              towerADCped, towerADCcha, preshoADCped, preshoADCcha );
-       printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
+       Info("FillFromDigits","Number of Digits: %d %d\n Parameters: PRE : %f %f EC: %f %f HC: %f %f\n Geometry: %d %d", 
+            ndig, nent, preADCped, preADCcha, ecADCped, ecADCcha, hcADCped, hcADCcha, geom->GetNEta(), geom->GetNPhi());
     }
     
 //
@@ -985,22 +1031,29 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag)
     TIter next(digs);
     while ((sdg = (AliEMCALDigit*)(next())))
     {
-       Double_t pedestal;
-       Double_t channel;
-       if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi())) 
-       {
-           pedestal = preshoADCped;
-           channel  = preshoADCcha; 
-       } else {
-           pedestal = towerADCped;
-           channel  = towerADCcha; 
+       Double_t pedestal = 0.;
+       Double_t channel  = 0.;
+       if (geom->IsInPRE(sdg->GetId())) {
+         pedestal = preADCped;
+         channel  = preADCcha; 
+       } 
+       else if (geom->IsInECAL(sdg->GetId()))  {
+         pedestal = ecADCped;
+         channel  = ecADCcha; 
+       }
+       else if (geom->IsInHCAL(sdg->GetId()))  {
+         pedestal = hcADCped;
+         channel  = hcADCcha; 
        }
+       else 
+         Fatal("FillFromDigits", "unexpected digit is number!") ; 
        
        Float_t eta = sdg->GetEta();
        Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
        Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
        
-       if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
+       if (fDebug > 1) 
+         Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
                           eta, phi, amp, sdg->GetAmp());
        
        fLego->Fill(eta, phi, fSamplingF*amp);
@@ -1430,7 +1483,7 @@ void AliEMCALJetFinder::FindTracksInJetCone()
        Int_t nT0 = nT + nTB;
        printf("Total number of tracks %d\n", nT0);
        
-       if (nT0 > 50) nT0 = 50;
+       if (nT0 > 1000) nT0 = 1000;
        
        Float_t* ptT  = new Float_t[nT0];
        Float_t* etaT = new Float_t[nT0];
@@ -1505,14 +1558,11 @@ Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curl
 // Propagates phi angle to EMCAL radius
 //
   static Float_t b = 0.0, rEMCAL = -1.0;
-  if(rEMCAL<0) {
 // Get field in kGS
-    b =  gAlice->Field()->SolenoidField();
+  b =  gAlice->Field()->SolenoidField();
 // Get EMCAL radius in cm 
-    rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
-    printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
-  }
-    Float_t dPhi = 0.;
+  rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
+  Float_t dPhi = 0.;
 //
 //
 // bending radies
@@ -1791,7 +1841,7 @@ void AliEMCALJetFinder::FindChargedJet()
        Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
        Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
        Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
-       Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
+//     Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
        
 //
 //      Sum up all energy