updates from salvatore

diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.cxx b/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.cxx
index cc7c7dd79d185d1cc0c9fc9fec922c25243b611b..7f82f96339a382f6d45a75220d75d540fbc1dac1 100644 (file)
--- a/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.cxx
+++ b/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.cxx
@@ -39,11 +39,13 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
   fKtJetsName("KtJets"),
+  fEmbJetsName("EmbJets"),
   fTrgClusName("ClustersL1GAMMAFEE"),
   fTracks(0),
   fCaloClusters(0),
   fJets(0),
   fKtJets(0),
+  fEmbJets(0),
   fTrgClusters(0),
   fCent(0),
   fCentBin(-1),
@@ -78,6 +80,9 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
     fHistDeltaPtRCExLJ[i] = 0;
     fHistRCPt[i] = 0;
     fHistRCPtExLJ[i] = 0;
+    fHistEmbJets[i] = 0;
+    fHistEmbPart[i] = 0;
+    fHistDeltaPtMedKtEmb[i] = 0;
   }
 
   // Output slot #1 writes into a TH1 container
@@ -98,11 +103,13 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
   fKtJetsName("KtJets"),
+  fEmbJetsName("EmbJets"),
   fTrgClusName("ClustersL1GAMMAFEE"),
   fTracks(0),
   fCaloClusters(0),
   fJets(0),
   fKtJets(0),
+  fEmbJets(0),
   fTrgClusters(0),
   fCent(0),
   fCentBin(-1),
@@ -137,6 +144,9 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
     fHistDeltaPtRCExLJ[i] = 0;
     fHistRCPt[i] = 0;
     fHistRCPtExLJ[i] = 0;
+    fHistEmbJets[i] = 0;
+    fHistEmbPart[i] = 0;
+    fHistDeltaPtMedKtEmb[i] = 0;
   }
 
   // Output slot #1 writes into a TH1 container
@@ -310,6 +320,27 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistRCPtExLJ[i]);
+
+    histname = "fHistEmbJets_";
+    histname += i;
+    fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
+    fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet P_{T} [GeV/c]");
+    fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistEmbJets[i]);
+
+    histname = "fHistEmbPart_";
+    histname += i;
+    fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
+    fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle P_{T} [GeV/c]");
+    fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistEmbPart[i]);
+
+    histname = "fHistDeltaPtMedKtEmb_";
+    histname += i;
+    fHistDeltaPtMedKtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
+    fHistDeltaPtMedKtEmb[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
+    fHistDeltaPtMedKtEmb[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtMedKtEmb[i]);
   }
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
@@ -339,6 +370,13 @@ void AliAnalysisTaskSAJF::RetrieveEventObjects()
       AliWarning(Form("Could not retrieve Kt jets!")); 
     }
   }
+
+  if (strcmp(fEmbJetsName,"")) {
+    fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
+    if (!fEmbJets) {
+      AliWarning(Form("Could not retrieve emb jets!")); 
+    }
+  }
     
   if (strcmp(fTrgClusName,"")) {
     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
@@ -437,6 +475,24 @@ Int_t AliAnalysisTaskSAJF::GetNumberOfKtJets() const
     return 0;
 }
 
+//________________________________________________________________________
+AliEmcalJet* AliAnalysisTaskSAJF::GetEmbJet(const Int_t i) const
+{
+  if (fEmbJets)
+    return dynamic_cast<AliEmcalJet*>(fEmbJets->At(i));
+  else
+    return 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
+{
+  if (fEmbJets)
+    return fEmbJets->GetEntriesFast();
+  else
+    return 0;
+}
+
 //________________________________________________________________________
 AliVCluster* AliAnalysisTaskSAJF::GetTrgCluster(const Int_t i) const
 {
@@ -487,7 +543,7 @@ void AliAnalysisTaskSAJF::FillHistograms()
     jet2 = GetJet(max2JetIndex);
 
   if (jet2) {
-    fHistLeadingJetPt[fCentBin]->Fill(jet2->Pt());
+    //fHistLeadingJetPt[fCentBin]->Fill(jet2->Pt());
     jet2->SortConstituents();
   }
 
@@ -520,6 +576,21 @@ void AliAnalysisTaskSAJF::FillHistograms()
     fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
     fHistRCPtVSRhoPart->Fill(RCptExLJ / A, rhoClus + rhoTracks);
   }
+
+  Float_t maxEmbJetPt = 0;
+  Float_t maxEmbPartPt = 0;
+
+  Bool_t embOK = DoEmbJetLoop(maxEmbJetPt, maxEmbPartPt);
+
+  if (embOK) {
+    fHistEmbJets[fCentBin]->Fill(maxEmbJetPt);
+    fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
+    fHistDeltaPtMedKtEmb[fCentBin]->Fill(maxEmbJetPt - A * rhoKt - maxEmbPartPt);
+  }
+  else {
+    if (maxEmbPartPt != 0)
+      AliWarning("Embedded particle is not the leading particle of the leading jet!");
+  }
 }
 
 //________________________________________________________________________
@@ -529,7 +600,6 @@ void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 
   Int_t njets = GetNumberOfJets();
-  //cout << "num of jets = " <<  njets << endl;
 
   Float_t maxJetPt = 0;
   Float_t max2JetPt = 0;
@@ -587,8 +657,6 @@ Float_t AliAnalysisTaskSAJF::DoKtJetLoop()
   Float_t ktJetsMedian = 0;
   Int_t nktjets =  GetNumberOfKtJets();
 
-  //cout << "num of ktjets = " <<  nktjets << endl;
-
   Int_t NoOfZeroJets = 0;
   if (nktjets > 0) {
     
@@ -636,8 +704,86 @@ Float_t AliAnalysisTaskSAJF::DoKtJetLoop()
 }
 
 //________________________________________________________________________
-Float_t AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex)
+Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(Float_t &maxJetPt, Float_t &maxPartPt)
 {
+  Double_t vertex[3] = {0, 0, 0};
+  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+
+  Int_t nembjets = GetNumberOfEmbJets();
+
+  Int_t maxJetIdx = -1;
+  Int_t maxTrackIdx = -1;
+  Int_t maxClusIdx = -1;
+
+  Float_t maxTrackPt = 0;
+  Float_t maxClusEt = 0;
+
+  for (Int_t ij = 0; ij < nembjets; ij++) {
+      
+    AliEmcalJet* jet = GetEmbJet(ij);
+      
+    if (!jet) {
+      AliError(Form("Could not receive jet %d", ij));
+      continue;
+    } 
+      
+    if (jet->Pt() <= 0)
+       continue;
+ 
+    if (!AcceptJet(jet))
+      continue;
+    
+    if (jet->Pt() > maxJetPt) {
+      maxJetPt = jet->Pt();
+      maxJetIdx = ij;
+    }
+  }
+
+  if (maxJetPt <= 0)
+    return kFALSE;
+
+  AliEmcalJet* jet = GetEmbJet(maxJetIdx);
+
+  for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
+    AliVTrack *track = jet->TrackAt(it, fTracks);
+
+    if (!track) continue;
+     
+    if (track->Pt() > maxTrackPt) {
+      maxTrackPt = track->Pt();
+      maxTrackIdx = it;
+    }
+  }
+
+  for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
+    AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
+    
+    if (!cluster) continue;
+    
+    TLorentzVector nPart;
+    cluster->GetMomentum(nPart, vertex);
+
+    if (nPart.Et() > maxClusEt) {
+      maxClusEt = nPart.Et();
+      maxClusIdx = ic;
+    } 
+  }
+
+  if (maxClusEt > maxTrackPt) {
+    AliVCluster *cluster = jet->ClusterAt(maxClusIdx, fCaloClusters);
+    maxPartPt = maxClusEt;
+    return (Bool_t)(cluster->Chi2() == 100);
+  }
+  else {
+    AliVTrack *track = jet->TrackAt(maxTrackIdx, fTracks);
+    maxPartPt = maxTrackPt;
+    return (Bool_t)(track->GetLabel() == 100);
+  }
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex)
+{ 
   AliEmcalJet* jet = GetJet(maxJetIndex);
   AliEmcalJet* jet2 = 0;
   if (max2JetIndex >= 0)
@@ -645,6 +791,7 @@ Float_t AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex)
 
   Float_t rhoTracks = 0;
   Int_t ntracks = GetNumberOfTracks();
+
   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
     AliVTrack* track = GetTrack(iTracks);         
     if(!track) {
@@ -712,6 +859,8 @@ Float_t AliAnalysisTaskSAJF::DoClusterLoop(Int_t maxJetIndex, Int_t max2JetIndex
     
     if (!(cluster->IsEMCAL())) continue;
     
+    if (!AcceptCluster(cluster)) continue;
+
     TLorentzVector nPart;
     cluster->GetMomentum(nPart, vertex);
 
@@ -739,6 +888,8 @@ Float_t AliAnalysisTaskSAJF::DoClusterLoop(Int_t maxJetIndex, Int_t max2JetIndex
        }  
        
        if (!(cluster2->IsEMCAL())) continue;
+
+       if (!AcceptCluster(cluster)) continue;
        
        if (IsJetCluster(jet, ic2)) continue;
        
@@ -820,6 +971,8 @@ Float_t AliAnalysisTaskSAJF::GetRigidConePt(AliEmcalJet *jet, Float_t minD)
     }  
     
     if (!(cluster->IsEMCAL())) continue;
+
+    if (!AcceptCluster(cluster)) continue;
     
     TLorentzVector nPart;
     cluster->GetMomentum(nPart, vertex);
@@ -918,8 +1071,21 @@ Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet* jet) const
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track) const
+Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC) const
+{
+  if (!acceptMC && clus->Chi2() == 100)
+    return kFALSE;
+
+  return kTRUE;
+}
+
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track, Bool_t acceptMC) const
 {
+  if (!acceptMC && track->GetLabel() == 100)
+    return kFALSE;
+  
   if (fAnaType == kFullAcceptance) {
     return kTRUE;
   }

diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.h b/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.h
index 4ce135fc6c93eef1373a1c29e72e6ad8c653c26b..dfcfee4822451d859fa1a7bc2264208e76d6a2cf 100644 (file)
--- a/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.h
+++ b/PWGGA/EMCALTasks/AliAnalysisTaskSAJF.h
@@ -32,10 +32,11 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   void                        Terminate(Option_t *option);
   void                        Init();
 
+  void                        SetTracksName(const char *n)                         { fTracksName    = n          ; }
   void                        SetClusName(const char *n)                           { fCaloName      = n          ; }
   void                        SetJetsName(const char *n)                           { fJetsName      = n          ; }
   void                        SetKtJetsName(const char *n)                         { fKtJetsName    = n          ; }
-  void                        SetTracksName(const char *n)                         { fTracksName    = n          ; }
+  void                        SetEmbJetsName(const char *n)                        { fEmbJetsName   = n          ; }
   void                        SetTrgClusName(const char *n)                        { fTrgClusName   = n          ; }
   void                        SetAnaType(SAJFAnaType type)                         { fAnaType       = type       ; }
   void                        SetJetRadius(Float_t r)                              { fJetRadius     = r          ; } 
@@ -43,47 +44,53 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
 
  protected:
 
-  AliVTrack                  *GetTrack(const Int_t i)               const;
-  Int_t                       GetNumberOfTracks()                   const;
-  AliVCluster                *GetCaloCluster(const Int_t i)         const;
-  Int_t                       GetNumberOfCaloClusters()             const;
-  AliEmcalJet                *GetJet(const Int_t i)                 const;
-  Int_t                       GetNumberOfJets()                     const;
-  AliEmcalJet                *GetKtJet(const Int_t i)               const;
-  Int_t                       GetNumberOfKtJets()                   const;
-  AliVCluster                *GetTrgCluster(const Int_t i)          const;
-  Int_t                       GetNumberOfTrgClusters()              const;
-  void                        FillHistograms()                           ;
-  void                        RetrieveEventObjects()                     ;
-  Bool_t                      AcceptTrack(AliVTrack* track)         const;
-  Bool_t                      AcceptJet(AliEmcalJet* jet)           const;
-  Bool_t                      AcceptJet(Float_t eta, Float_t phi)   const;
+  AliVTrack                  *GetTrack(const Int_t i)                                              const;
+  Int_t                       GetNumberOfTracks()                                                  const;
+  AliVCluster                *GetCaloCluster(const Int_t i)                                        const;
+  Int_t                       GetNumberOfCaloClusters()                                            const;
+  AliEmcalJet                *GetJet(const Int_t i)                                                const;
+  Int_t                       GetNumberOfJets()                                                    const;
+  AliEmcalJet                *GetKtJet(const Int_t i)                                              const;
+  Int_t                       GetNumberOfKtJets()                                                  const;
+  AliEmcalJet                *GetEmbJet(const Int_t i)                                             const;
+  Int_t                       GetNumberOfEmbJets()                                                 const;
+  AliVCluster                *GetTrgCluster(const Int_t i)                                         const;
+  Int_t                       GetNumberOfTrgClusters()                                             const;
+  Bool_t                      AcceptTrack(AliVTrack* track, Bool_t acceptMC = kFALSE)              const;
+  Bool_t                      AcceptCluster(AliVCluster* clus, Bool_t acceptMC = kFALSE)           const;
+  Bool_t                      AcceptJet(AliEmcalJet* jet)                                          const;
+  Bool_t                      AcceptJet(Float_t eta, Float_t phi)                                  const;
   Bool_t                      IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted = kTRUE)    const;
   Bool_t                      IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted = kTRUE)   const;
-  Float_t                     GetArea()                        const;
-  void                        DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex);
-  Float_t                     DoKtJetLoop();
-  Float_t                     DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex);
-  Float_t                     DoClusterLoop(Int_t maxJetIndex, Int_t max2JetIndex);
-  Float_t                     GetRigidConePt(AliEmcalJet *jet = 0,  Float_t minD = 0.8);
+  Float_t                     GetArea()                                                            const;
+  void                        RetrieveEventObjects()                                                    ;
+  void                        FillHistograms()                                                          ;
+  void                        DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)                        ;
+  Float_t                     DoKtJetLoop()                                                             ;
+  Bool_t                      DoEmbJetLoop(Float_t &maxJetPt, Float_t &maxPartPt)                       ;
+  Float_t                     DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex)                        ;
+  Float_t                     DoClusterLoop(Int_t maxJetIndex, Int_t max2JetIndex)                      ;
+  Float_t                     GetRigidConePt(AliEmcalJet *jet = 0,  Float_t minD = 0.8)                 ;
 
 
-  SAJFAnaType                 fAnaType;                    // analysis type
-  Float_t                     fMinEta;                     // minimum eta accepatance
-  Float_t                     fMaxEta;                     // maximum eta accepatance
-  Float_t                     fMinPhi;                     // minimum phi accepatance
-  Float_t                     fMaxPhi;                     // maximum phi accepatance  
+  SAJFAnaType                 fAnaType;                    // Analysis type
+  Float_t                     fMinEta;                     // Minimum eta accepatance
+  Float_t                     fMaxEta;                     // Maximum eta accepatance
+  Float_t                     fMinPhi;                     // Minimum phi accepatance
+  Float_t                     fMaxPhi;                     // Maximum phi accepatance  
   Float_t                     fJetRadius;                  // jet radius
   TList                      *fOutput;                     // Output list
-  TString                     fTracksName;                 // name of track collection
-  TString                     fCaloName;                   // name of calo cluster collection
-  TString                     fJetsName;                   // name of jet collection
-  TString                     fKtJetsName;                 // name of kt jet collection
-  TString                     fTrgClusName;                // name of trg clus name
+  TString                     fTracksName;                 // Name of track collection
+  TString                     fCaloName;                   // Name of calo cluster collection
+  TString                     fJetsName;                   // Name of jet collection
+  TString                     fKtJetsName;                 // Name of kt jet collection
+  TString                     fEmbJetsName;                // Name of embedded jets collection
+  TString                     fTrgClusName;                // Name of trg clus name
   TClonesArray               *fTracks;                     //!Tracks
   TClonesArray               *fCaloClusters;               //!Clusters
   TClonesArray               *fJets;                       //!Jets
   TClonesArray               *fKtJets;                     //!Kt Jets
+  TClonesArray               *fEmbJets;                    //!Embedded Jets
   TClonesArray               *fTrgClusters;                //!Trg Clusters
   AliCentrality              *fCent;                       // Event centrality
   Int_t                       fCentBin;                    // Event centrality bin
@@ -111,6 +118,9 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   TH1F                       *fHistDeltaPtRCExLJ[4];       // deltaPt = Pt(RC) - A * rhoKt, imposing min distance from leading jet
   TH1F                       *fHistRCPt[4];                // Random cone pt
   TH1F                       *fHistRCPtExLJ[4];            // Random cone pt, imposing min distance from leading jet
+  TH1F                       *fHistEmbJets[4];             // Pt distribution of embedded jets
+  TH1F                       *fHistEmbPart[4];             // Pt distribution of embedded particle
+  TH1F                       *fHistDeltaPtMedKtEmb[4];     // deltaPt = Pt(embjet) - A * rhoKt - Pt(embtrack)
   Int_t                       fNbins;                      // No. of pt bins
   Float_t                     fMinPt;                      // Min pt
   Float_t                     fMaxPt;                      // Max pt
@@ -119,6 +129,6 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   AliAnalysisTaskSAJF(const AliAnalysisTaskSAJF&);            // not implemented
   AliAnalysisTaskSAJF &operator=(const AliAnalysisTaskSAJF&); // not implemented
 
-  ClassDef(AliAnalysisTaskSAJF, 1) // Isolated photons task
+  ClassDef(AliAnalysisTaskSAJF, 1) // jet finder task
 };
 #endif

diff --git a/PWGGA/EMCALTasks/macros/AddTaskSAJF.C b/PWGGA/EMCALTasks/macros/AddTaskSAJF.C
index 73553c65c8e1148e71cfe746bc20c8d51587fd0d..d08f8aebebb3b886f346bca620092c5a9df56a58 100644 (file)
--- a/PWGGA/EMCALTasks/macros/AddTaskSAJF.C
+++ b/PWGGA/EMCALTasks/macros/AddTaskSAJF.C
@@ -6,9 +6,10 @@ AliAnalysisTaskSAJF* AddTaskSAJF(
   const char *nclusters          = "CaloClusters",
   const char *njets              = "Jets",
   const char *nktjets            = "KtJets",
-  Float_t jetradius              = 0.4,
+  const char *nembjets           = "EmbJets",
+  Double_t    jetradius          = 0.4,
   const char *ntrgclusters       = "ClustersL1GAMMAFEE",
-  UInt_t type                    = AliAnalysisTaskSAJF::kEMCAL
+  UInt_t      type               = AliAnalysisTaskSAJF::kEMCAL
 )
 {  
   // Get the pointer to the existing analysis manager via the static access method.
@@ -38,6 +39,7 @@ AliAnalysisTaskSAJF* AddTaskSAJF(
   phTask->SetClusName(nclusters);
   phTask->SetJetsName(njets);
   phTask->SetKtJetsName(nktjets);
+  phTask->SetEmbJetsName(nembjets);
   phTask->SetTrgClusName(ntrgclusters);
   phTask->SetJetRadius(jetradius);