1) Use fiducial cuts in eta to select jets
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2005 07:56:48 +0000 (07:56 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jul 2005 07:56:48 +0000 (07:56 +0000)
   Modification to AliUA1JetFinder::FindJets() in  AliUA1JetFinder.cxx
   to select only jets within the acceptable eta range

2) Put a flag if accepted track belongs to pythia

   Modification to:

   AliJetReader.h to add an array to keep the flags and a method
   to get the array.
   AliJetReader.cxx to initialize the array in the constructur
   AliJetESDReader.cxx to fill it in AliJetESDReader::FillMomentumArray
   Flag set to 1 for particles whose label is less
   than 10000 (ie from Pythia)

3) Compute percentage of energy in jet coming from signal

   Modification to

   AliUA1JetFinder::FindJets() in  AliUA1JetFinder.cxx
   to compute the percentage

   AliJet.h to set an array to store the percentages and functions to
   set in the percentages and to get at them

   AliJet.cxx to implement the functions and initialize the array

4) Modify destructors

   AliUA1JetFinder.cxx: now it is an empty destructor
   (fLego is not deleted anymore)

   AliJetFinder.cxx: fLeading is not deleted anymore

   AliJetReaderHeader.cxx:  now it is an empty destructor
   (TStrings not deleted anymore)

(G. Contreras)

JETAN/AliJetESDReader.cxx
JETAN/AliJetReader.cxx
JETAN/AliJetReader.h
JETAN/AliJetReaderHeader.cxx
JETAN/AliUA1JetFinder.cxx

index 956bb06..4a2ce34 100755 (executable)
@@ -45,7 +45,7 @@ AliJetESDReader::AliJetESDReader()
 AliJetESDReader::~AliJetESDReader()
 {
   // Destructor
-  delete fReaderHeader;
+  //  delete fReaderHeader;
 }
 
 //____________________________________________________________________________
@@ -108,6 +108,9 @@ void AliJetESDReader::FillMomentumArray(Int_t event)
   fChainMC->GetEntry(event);
   // get number of tracks in event (for the loop)
   nt = fESD->GetNumberOfTracks();
+  // tmporary storage of signal flags
+  Int_t* flag  = new Int_t[nt];
+
   // get cuts set by user
   Float_t ptMin = ((AliJetESDReaderHeader*) fReaderHeader)->GetPtCut();
 
@@ -117,7 +120,8 @@ void AliJetESDReader::FillMomentumArray(Int_t event)
       UInt_t status = track->GetStatus();
       if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
 
-      if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() && TMath::Abs(track->GetLabel()) > 10000)  continue;
+      if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
+         && TMath::Abs(track->GetLabel()) > 10000)  continue;
     
       track->GetPxPyPz(p);
       pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]); // pt of the track
@@ -127,8 +131,12 @@ void AliJetESDReader::FillMomentumArray(Int_t event)
       new ((*fMomentumArray)[goodTrack]) 
          TLorentzVector(p[0], p[1], p[2],
                         TMath::Sqrt(pt * pt +p[2] * p[2]));
+      flag[goodTrack]=0;
+      if (TMath::Abs(track->GetLabel()) < 10000) flag[goodTrack]=1;
       goodTrack++;
   }
+  // set the signal flags
+  fSignalFlag.Set(goodTrack,flag);
 
   printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
 }
index d60fc23..3835a8a 100755 (executable)
@@ -39,6 +39,7 @@ AliJetReader::AliJetReader()
   fMomentumArray = new TClonesArray("TLorentzVector",2000);
   fArrayMC = 0;
   fAliHeader = 0;
+  fSignalFlag = TArrayI();
 }
 
 ////////////////////////////////////////////////////////////////////////
index 1dc85a0..cb8979d 100755 (executable)
@@ -12,6 +12,7 @@
   
 #include <TObject.h>
 #include <TChain.h>
+#include <TArrayI.h>
 
 class TClonesArray;
 class AliJetReaderHeader;
@@ -30,10 +31,12 @@ class AliJetReader : public TObject
   virtual Int_t GetChainEntries() {return fChain->GetEntries();} 
   virtual AliJetReaderHeader* GetReaderHeader() { return fReaderHeader;}
   virtual AliHeader* GetAliHeader() { return fAliHeader;}
+  virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];}
 
   // Setters
   virtual void FillMomentumArray(Int_t) {}
-  virtual void SetReaderHeader(AliJetReaderHeader* header) {fReaderHeader = header;}
+  virtual void SetReaderHeader(AliJetReaderHeader* header) 
+    {fReaderHeader = header;}
          
   // others
   virtual void OpenInputFiles() {}
@@ -47,6 +50,9 @@ class AliJetReader : public TObject
   AliESD                  *fESD;           // pointer to esd
   AliJetReaderHeader      *fReaderHeader;  // pointer to header
   AliHeader               *fAliHeader;     // pointer to event header
+  TArrayI fSignalFlag;   // to flag if a particle comes from pythia or 
+                        // from the underlying event
+
   ClassDef(AliJetReader,1)
 };
  
index 193857c..cf39a0a 100755 (executable)
@@ -58,7 +58,7 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name):
 AliJetReaderHeader::~AliJetReaderHeader()
 {
   // destructor
-  delete fComment;
-  delete fDir;
-  delete fPattern;
+  //  delete fComment;
+  //delete fDir;
+  //delete fPattern;
 }
index 946bd4d..2e291de 100755 (executable)
@@ -28,6 +28,7 @@
 #include "AliUA1JetFinder.h"
 #include "AliUA1JetHeader.h"
 #include "UA1Common.h"
+#include "AliJetReaderHeader.h"
 #include "AliJetReader.h"
 #include "AliJet.h"
 
@@ -55,7 +56,7 @@ AliUA1JetFinder::~AliUA1JetFinder()
   // destructor
   //
 
-  delete fLego;            
+  // delete fLego;            
   
   // reset and delete header
 }
@@ -148,38 +149,52 @@ void AliUA1JetFinder::FindJets()
   ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, 
                 minmove, maxmove, mode, precbg, ierror);
 
-  // download jet info. 
-  Int_t* mult  = new Int_t[UA1JETS.njet];
+
   // sort jets
   Int_t * idx  = new Int_t[UA1JETS.njet];
   TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
+
+  // download jet info.   
+  AliJetReaderHeader* rheader = fReader->GetReaderHeader();
+  for(Int_t i = 0; i < UA1JETS.njet; i++) {
+    // reject events outside acceptable eta range
+    if (((UA1JETS.etaj[1][idx[i]])> (rheader->GetFiducialEtaMax()))
+       || ((UA1JETS.etaj[1][idx[i]]) < (rheader->GetFiducialEtaMin())))
+      continue;
+    Float_t px, py,pz,en; // convert to 4-vector
+    px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
+    py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
+    pz = UA1JETS.etj[idx[i]] /
+      TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
+    en = TMath::Sqrt(px * px + py * py + pz * pz);
+    fJets->AddJet(px, py, pz, en);
+  }
   
+  // find multiplicities and relationship jet-particle
+  // find also percentage of pt from pythia
   Int_t* injet = new Int_t[nIn];
   for (Int_t i = 0; i < nIn; i++) injet[i]= -1;
-
-  for(Int_t i = 0; i < UA1JETS.njet; i++) {
-      Float_t px, py,pz,en; // convert to 4-vector
-      px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
-      py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
-      pz = UA1JETS.etj[idx[i]] /
-          TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
-      en = TMath::Sqrt(px * px + py * py + pz * pz);
-      
-      fJets->AddJet(px, py, pz, en);
-      // find multiplicities and relationship jet-particle
-      mult[i] = 0;
-      for (Int_t j = 0; j < nIn; j++) {
-         Float_t deta = etaT[j] - UA1JETS.etaj[1][idx[i]];
-         Float_t dphi = phiT[j] - UA1JETS.phij[1][idx[i]];
-         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
-         if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
-         Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
-         if (dr < fHeader->GetRadius() && injet[j] == -1) {
-             injet[j] = i;
-             mult[i]++;
-         }
+  Int_t* mult  = new Int_t[fJets->GetNJets()];
+  Float_t* percentage  = new Float_t[fJets->GetNJets()];
+
+  for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
+    Float_t pt_sig = 0.0;
+    mult[i] = 0;
+    for (Int_t j = 0; j < nIn; j++) {
+      Float_t deta = etaT[j] - fJets->GetEta(i);
+      Float_t dphi = phiT[j] - fJets->GetPhi(i);
+      if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+      if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+      Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+      if (dr < fHeader->GetRadius() && injet[j] == -1) {
+       if (fReader->GetSignalFlag(j) == 1) pt_sig+=ptT[j];
+       injet[j] = i;
+       mult[i]++;
       }
+    }
+    percentage[i] = pt_sig/((Double_t) fJets->GetPt(i));    
   }
+  fJets->SetPtFromSignal(percentage);
   fJets->SetMultiplicities(mult);
   fJets->SetInJet(injet);
 }