]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
hadronic correction and scale task from Rosi
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 21:10:14 +0000 (21:10 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 21:10:14 +0000 (21:10 +0000)
PWGGA/CMakelibPWGGAEMCALJetTasks.pkg
PWGGA/EMCALJetTasks/AliAnalysisTaskScale.cxx [new file with mode: 0644]
PWGGA/EMCALJetTasks/AliAnalysisTaskScale.h [new file with mode: 0644]
PWGGA/EMCALJetTasks/AliHadCorrTask.cxx [new file with mode: 0644]
PWGGA/EMCALJetTasks/AliHadCorrTask.h [new file with mode: 0644]
PWGGA/PWGGAEMCALJetTasksLinkDef.h

index a2b63c96dec4926d2269d98c9a27380f01ae8626..f050c5243cf483b38d5bafc682ea4247851c8eeb 100644 (file)
@@ -25,7 +25,9 @@
 #--------------------------------------------------------------------------------#
 
 set ( SRCS 
+ EMCALJetTasks/AliAnalysisTaskScale.cxx
  EMCALJetTasks/AliEmcalJetTask.cxx
+ EMCALJetTasks/AliHadCorrTask.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGGA/EMCALJetTasks/AliAnalysisTaskScale.cxx b/PWGGA/EMCALJetTasks/AliAnalysisTaskScale.cxx
new file mode 100644 (file)
index 0000000..1cbde06
--- /dev/null
@@ -0,0 +1,159 @@
+// $Id$
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TList.h"
+#include "TH1.h"
+#include "TH2.h"
+#include <TClonesArray.h>
+#include <TParticle.h>
+#include <TLorentzVector.h>
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliVCluster.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliCentrality.h"
+
+#include "AliAnalysisTaskScale.h"
+
+ClassImp(AliAnalysisTaskScale)
+
+//________________________________________________________________________
+AliAnalysisTaskScale::AliAnalysisTaskScale(const char *name) 
+  : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), fHistCentrality(0), fHistPtTPCvsCent(0), fHistPtEMCALvsCent(0), fHistEtvsCent(0),  fHistScalevsCent(0),  fHistDeltaScalevsCent(0), fHistPtTPCvsNtrack(0), fHistPtEMCALvsNtrack(0), fHistEtvsNtrack(0),  fHistScalevsNtrack(0),  fHistDeltaScalevsNtrack(0),
+    fTracksName("tracks"),
+    fClustersName("clusters")
+
+{
+  // Constructor
+
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScale::UserCreateOutputObjects()
+{
+  OpenFile(1);
+  fOutputList = new TList();
+  fOutputList->SetOwner();
+
+  fHistCentrality = new TH1F("fHistCentrality","Centrality",101,-1,100);
+
+  fHistPtTPCvsCent = new TH2F("fHistPtTPCvsCent","rho vs cent",101,-1,100,500,0,1000);
+  fHistPtEMCALvsCent = new TH2F("fHistPtEMCALvsCent","rho vs cent",101,-1,100,500,0,1000);
+  fHistEtvsCent = new TH2F("fHistEtvsCent","rho vs cent",101,-1,100,500,0,1000);
+  fHistScalevsCent = new TH2F("fHistScalevsCent","rho vs cent",101,-1,100,400,0,4.0);
+  fHistDeltaScalevsCent = new TH2F("fHistDeltaScalevsCent","rho vs cent",101,-1,100,400,-2.0,2.0);
+  
+  
+  fHistPtTPCvsNtrack = new TH2F("fHistPtTPCvsNtrack","rho vs cent",500,0,2500,500,0,1000);
+  fHistPtEMCALvsNtrack = new TH2F("fHistPtEMCALvsNtrack","rho vs cent",500,0,2500,500,0,1000);
+  fHistEtvsNtrack = new TH2F("fHistEtvsNtrack","rho vs cent",500,0,2500,500,0,1000);
+  fHistScalevsNtrack = new TH2F("fHistScalevsNtrack","rho vs cent",500,0,2500,400,0,4.0);
+  fHistDeltaScalevsNtrack = new TH2F("fHistDeltaScalevsNtrack","rho vs cent",500,0,2500,400,-2.0,2.0);
+    
+  fOutputList->Add(fHistCentrality);
+  fOutputList->Add(fHistPtTPCvsCent);
+  fOutputList->Add(fHistPtEMCALvsCent);
+  fOutputList->Add(fHistEtvsCent);
+  fOutputList->Add(fHistScalevsCent);
+  fOutputList->Add(fHistDeltaScalevsCent);
+  fOutputList->Add(fHistPtTPCvsNtrack);
+  fOutputList->Add(fHistPtEMCALvsNtrack);
+  fOutputList->Add(fHistEtvsNtrack);
+  fOutputList->Add(fHistScalevsNtrack);
+  fOutputList->Add(fHistDeltaScalevsNtrack);
+  
+
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskScale::UserExec(Option_t *) 
+{
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!fESD) {
+    printf("ERROR: fESD not available\n");
+    return;
+  }
+
+  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+  TClonesArray *tracks = 0;
+  TClonesArray *clusters = 0;
+  AliCentrality *centrality = 0;
+  TList *l = InputEvent()->GetList();
+
+  tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
+  if (!tracks) {
+    AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
+    return;
+  }
+    
+  clusters = dynamic_cast<TClonesArray*>(l->FindObject(fClustersName));
+  if (!clusters){
+    AliError(Form("Pointer to clusters %s == 0", fClustersName.Data() ));
+    return;
+  }
+  
+
+  centrality = dynamic_cast<AliCentrality*>(l->FindObject("Centrality"));
+  float fCent = centrality->GetCentralityPercentile("V0M");
+  fHistCentrality->Fill(fCent);
+
+  float scale = 0.000066*fCent*fCent-0.0015*fCent+1.5; //march 25th fit
+  float ptTPC = 0;float ptEMCAL=0;
+  const Int_t Ntracks = tracks->GetEntries();
+  for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
+    AliVTrack *track = static_cast<AliVTrack*>(tracks->At(iTracks));
+    if (!track)
+      continue;
+    if (fabs(track->Eta())>0.7)
+      continue;
+    ptTPC+=track->Pt();
+    if ((track->Phi()>3.14)||(track->Phi()<1.4))
+      continue;
+    ptEMCAL+=track->Pt();
+  } //track loop 
+  
+  Double_t vertex[3] = {0, 0, 0};
+  
+  float Et = 0;
+  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+  const Int_t Nclus = clusters->GetEntries();
+  for (Int_t iClus = 0, iN = 0; iClus < Nclus; ++iClus) {
+    
+    AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(iClus));
+    if (!c->IsEMCAL())
+      continue;
+    TLorentzVector nPart;
+    c->GetMomentum(nPart, vertex);
+    Et += nPart.P();
+  }
+  
+  float scalecalc = ((Et+ptEMCAL)/2.44)*(8.8/ptTPC); //2.44 = EMCAL Acc, 8.8 = TPC Acc
+  fHistPtTPCvsCent->Fill(fCent,ptTPC);
+  fHistPtEMCALvsCent->Fill(fCent,ptEMCAL);
+  fHistEtvsCent->Fill(fCent,Et);
+  fHistScalevsCent->Fill(fCent,scalecalc);
+  fHistDeltaScalevsCent->Fill(fCent,scalecalc-scale);
+
+  fHistPtTPCvsNtrack->Fill(Ntracks,ptTPC);
+  fHistPtEMCALvsNtrack->Fill(Ntracks,ptEMCAL);
+  fHistEtvsNtrack->Fill(Ntracks,Et);
+  fHistScalevsNtrack->Fill(Ntracks,scalecalc);
+  fHistDeltaScalevsNtrack->Fill(Ntracks,scalecalc-scale);
+
+  PostData(1, fOutputList);
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskScale::Terminate(Option_t *) 
+{}
+
+
+
+
diff --git a/PWGGA/EMCALJetTasks/AliAnalysisTaskScale.h b/PWGGA/EMCALJetTasks/AliAnalysisTaskScale.h
new file mode 100644 (file)
index 0000000..86f8796
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef AliAnalysisTaskScale_cxx
+#define AliAnalysisTaskScale_cxx
+
+// $Id$
+
+class TList;
+class TH1F;
+class TH1I;
+class TH2F;
+//class TH3F;
+class AliESDEvent;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskScale : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskScale() : AliAnalysisTaskSE(), fESD(0), fOutputList(0), fHistCentrality(0), fHistPtTPCvsCent(0), fHistPtEMCALvsCent(0), fHistEtvsCent(0),  fHistScalevsCent(0),  fHistDeltaScalevsCent(0), fHistPtTPCvsNtrack(0), fHistPtEMCALvsNtrack(0), fHistEtvsNtrack(0),  fHistScalevsNtrack(0),  fHistDeltaScalevsNtrack(0){}
+  AliAnalysisTaskScale(const char *name);
+  virtual ~AliAnalysisTaskScale() {}
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  virtual void   SetTracksName(const char *n)        { fTracksName = n; }
+  virtual void   SetClustersName(const char *n)          { fClustersName = n; }
+  
+ private:
+  TString                fTracksName;             // name of track collection
+  TString                fJetsName;             // name of jet collection
+  TString                fClustersName;             // name of clusters collection
+  AliESDEvent *fESD;    //! ESD object
+  TList       *fOutputList; //! Output list
+  TH1F        *fHistCentrality;
+
+  TH2F        *fHistPtTPCvsCent; 
+  TH2F        *fHistPtEMCALvsCent; 
+  TH2F        *fHistEtvsCent;  
+  TH2F        *fHistScalevsCent; 
+  TH2F        *fHistDeltaScalevsCent;  
+
+  TH2F        *fHistPtTPCvsNtrack; 
+  TH2F        *fHistPtEMCALvsNtrack; 
+  TH2F        *fHistEtvsNtrack;  
+  TH2F        *fHistScalevsNtrack; 
+  TH2F        *fHistDeltaScalevsNtrack;  
+   
+
+
+  AliAnalysisTaskScale(const AliAnalysisTaskScale&); // not implemented
+  AliAnalysisTaskScale& operator=(const AliAnalysisTaskScale&); // not implemented
+  
+  ClassDef(AliAnalysisTaskScale, 1); // example of analysis
+};
+
+#endif
diff --git a/PWGGA/EMCALJetTasks/AliHadCorrTask.cxx b/PWGGA/EMCALJetTasks/AliHadCorrTask.cxx
new file mode 100644 (file)
index 0000000..c0fab35
--- /dev/null
@@ -0,0 +1,246 @@
+// $Id$
+
+#include "AliHadCorrTask.h"
+#include <TClonesArray.h>
+#include <TParticle.h>
+#include "AliEmcalJet.h"
+#include "AliAnalysisManager.h"
+#include "AliESDtrack.h"
+#include "AliFJWrapper.h"
+#include "AliESDCaloCluster.h"
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TChain.h"
+#include <TLorentzVector.h>
+#include <AliCentrality.h>
+#include "AliPicoTrack.h"
+
+ClassImp(AliHadCorrTask)
+
+//________________________________________________________________________
+AliHadCorrTask::AliHadCorrTask(const char *name) : 
+  AliAnalysisTaskSE("AliHadCorrTask"),
+  fTracksName("Tracks"),
+  fCaloName("CaloClusters"),
+  fMinPt(0.15),
+  fOutputList(0),
+  fHistEbefore(0),
+  fHistEafter(0),
+  fHistNclusvsCent(0),
+  fHistNclusMatchvsCent(0),
+  fOutCaloName("CaloClustersOut"),
+  fOutClusters(0)
+{
+
+  for(Int_t i=0; i<4; i++){
+    fHistMatchEtaPhi[i]=0x0;
+    fHistMatchEvsP[i]=0x0;
+    fHistMatchdRvsEP[i]=0x0;
+  }
+
+  // Standard constructor.
+  cout << "Constructor for HadCorrTask " << name <<endl;
+  if (!name)
+    return;
+
+  SetName(name);
+  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
+
+  DefineInput(0,TChain::Class());
+  DefineOutput(1,TList::Class());
+}
+
+AliHadCorrTask::AliHadCorrTask() : 
+  AliAnalysisTaskSE("AliHadCorrTask"),
+  fTracksName("Tracks"),
+  fCaloName("CaloClusters"),
+  fOutClusters(0)
+{
+  // Standard constructor.
+
+  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
+}
+
+
+//________________________________________________________________________
+AliHadCorrTask::~AliHadCorrTask()
+{
+  // Destructor
+}
+
+//________________________________________________________________________
+void AliHadCorrTask::UserCreateOutputObjects()
+{
+
+  fOutClusters = new TClonesArray("AliESDCaloCluster");
+  fOutClusters->SetName(fOutCaloName);
+  OpenFile(1);
+  fOutputList = new TList();
+  fOutputList->SetOwner();
+
+  char name[200];
+
+  for(int icent=0; icent<4; icent++){
+    sprintf(name,"fHistMatchEtaPhi_%i",icent);
+    fHistMatchEtaPhi[icent] = new TH2F(name,name,400,-0.2,0.2,1600,-0.8,0.8);
+    sprintf(name,"fHistMatchEvsP_%i",icent);
+    fHistMatchEvsP[icent]=new TH2F(name,name,400,0.,200.,1000,0.,10.);
+    sprintf(name,"fHistMatchdRvsEP_%i",icent);
+    fHistMatchdRvsEP[icent]=new TH2F(name,name,1000,0.,1.,1000,0.,10.);
+
+    
+  fOutputList->Add(fHistMatchEtaPhi[icent]);
+  fOutputList->Add(fHistMatchEvsP[icent]);
+  fOutputList->Add(fHistMatchdRvsEP[icent]);
+
+  }
+  fHistNclusvsCent=new TH1F("Nclusvscent","NclusVsCent",100,0,100);
+  fHistNclusMatchvsCent=new TH1F("NclusMatchvscent","NclusMatchVsCent",100,0,100);
+  fHistEbefore=new TH1F("Ebefore","Ebefore",100,0,100);
+  fHistEafter=new TH1F("Eafter","Eafter",100,0,100);
+
+  fOutputList->Add(fHistNclusMatchvsCent);
+  fOutputList->Add(fHistNclusvsCent);
+  fOutputList->Add(fHistEbefore);
+  fOutputList->Add(fHistEafter);
+
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliHadCorrTask::UserExec(Option_t *) 
+{
+
+  if (!(InputEvent()->FindListObject(fOutCaloName)))
+    InputEvent()->AddObject(fOutClusters);
+  else fOutClusters->Delete();
+
+
+  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+  TClonesArray *tracks = 0;
+  TClonesArray *clus   = 0;
+  TList *l = InputEvent()->GetList();
+
+  AliCentrality *centrality = 0;
+  centrality = dynamic_cast<AliCentrality*>(l->FindObject("Centrality"));
+  Float_t fCent=-1; 
+  if(centrality)
+    fCent = centrality->GetCentralityPercentile("V0M");
+  else
+    fCent=99;//probably pp data
+
+  if(fCent<0) return;
+
+  if (fTracksName == "Tracks")
+    am->LoadBranch("Tracks");
+  tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
+  if (!tracks) {
+    AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
+    return;
+  }
+  const Int_t Ntrks = tracks->GetEntries();
+  
+  if (fCaloName == "CaloClusters")
+    am->LoadBranch("CaloClusters");
+  clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
+  if (!clus) {
+    AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
+    return;
+  }
+
+  Int_t centbin=-1;
+  if(fCent>=0 && fCent<10) centbin=0;
+  else if(fCent>=10 && fCent<30) centbin=1;
+  else if(fCent>=30 && fCent<50) centbin=2;
+  else if(fCent>=50 && fCent<100) centbin=3;
+  if (clus) {
+    Double_t vertex[3] = {0, 0, 0};
+    InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+    const Int_t Nclus = clus->GetEntries();
+    for (Int_t iClus = 0, iN = 0, clusCount=0; iClus < Nclus; ++iClus) {
+      AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
+      if (!c)
+        continue;
+      if (!c->IsEMCAL())
+        continue;
+      TLorentzVector nPart;
+
+      c->SetEmcCpvDistance(-1);
+      c->SetTrackDistance(999,999);
+      Double_t dEtaMin  = 1e9;
+      Double_t dPhiMin  = 1e9;
+      Int_t    imin     = -1;
+      for(Int_t t = 0; t<Ntrks; ++t) {
+        AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(t));
+        Double_t etadiff=999;
+        Double_t phidiff=999;
+        AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
+        Double_t dR = TMath::Sqrt(etadiff*etadiff+phidiff*phidiff);
+        Double_t dRmin = TMath::Sqrt(dEtaMin*dEtaMin+dPhiMin*dPhiMin);
+        if(dR > 25) 
+          continue;
+       if(dR<dRmin){
+          dEtaMin = etadiff;
+          dPhiMin = phidiff;
+          imin = t;
+        }
+      }
+      c->SetEmcCpvDistance(imin);
+      c->SetTrackDistance(dPhiMin, dEtaMin);
+      c->GetMomentum(nPart, vertex);
+      Double_t energy = nPart.P();
+      if(energy<fMinPt) continue;
+      if (imin>=0) {
+       Double_t dPhiMin = c->GetTrackDx();
+       Double_t dEtaMin = c->GetTrackDz();
+       fHistMatchEtaPhi[centbin]->Fill(dEtaMin,dPhiMin);
+      }
+
+      fHistNclusvsCent->Fill(fCent);
+    
+      if (fHadCorr>0) {
+       if (imin>=0) {
+         Double_t dPhiMin = c->GetTrackDx();
+         Double_t dEtaMin = c->GetTrackDz();
+         Double_t dR=TMath::Sqrt(dEtaMin*dEtaMin+dPhiMin*dPhiMin);
+             
+         AliVTrack *t = dynamic_cast<AliVTrack*>(tracks->At(imin));
+         if (t) {
+           if (t->Pt()<fMinPt)
+             continue;
+           if (t->P()>0) fHistMatchEvsP[centbin]->Fill(energy,energy/t->P());
+           if (t->P()>0) fHistMatchdRvsEP[centbin]->Fill(dR,energy/t->P());
+           fHistEbefore->Fill(fCent,energy);
+           if (dPhiMin<0.05 && dEtaMin<0.025) { // pp cuts!!!
+              energy -= fHadCorr*t->P();
+             fHistNclusMatchvsCent->Fill(fCent);
+           }
+           if (energy<0)
+             continue;
+         }
+       }
+       fHistEafter->Fill(fCent,energy);
+
+      }//end had correction if
+      if (energy>0){//Output new corrected clusters
+       AliESDCaloCluster *oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*>(c)));
+       oc->SetE(energy);
+       clusCount++;
+      }
+    }
+
+
+  }
+}
+
+//________________________________________________________________________
+void AliHadCorrTask::Terminate(Option_t *) 
+{
+
+}
+
+//________________________________________________________________________
diff --git a/PWGGA/EMCALJetTasks/AliHadCorrTask.h b/PWGGA/EMCALJetTasks/AliHadCorrTask.h
new file mode 100644 (file)
index 0000000..3e5b600
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef ALIHADCORRTASK_H
+#define ALIHADCORRTASK_H
+
+// $Id$
+
+class TClonesArray;
+class AliESDtrackCuts;
+class TList;
+class TH1;
+class TH2;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliHadCorrTask : public AliAnalysisTaskSE {
+ public:
+  AliHadCorrTask(const char *name); 
+  AliHadCorrTask();
+  virtual ~AliHadCorrTask();
+
+  void         UserCreateOutputObjects();
+  void         UserExec(Option_t *option);
+  void         Terminate(Option_t *option);
+
+  void         SetClusName(const char *n)            { fCaloName   = n; }
+  void         SetHadCorr(Double_t c)                { fHadCorr    = c; }
+  void         SetTracksName(const char *n)          { fTracksName = n; }
+  void         SetMinPt(Double_t min)                { fMinPt = min;  }
+  void         SetOutClusName(const char *n)         { fOutCaloName = n;}
+
+ protected:
+  void FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t fCent);
+
+  TString                fTracksName;             // name of track collection
+  TString                fCaloName;               // name of calo cluster collection
+
+  Double_t               fHadCorr;                // hadronic correction
+  TClonesArray          *fJets;                   //!jet collection
+  Double_t               fMinPt; 
+  TString                fOutCaloName;            //name of output clusters
+  TClonesArray          *fOutClusters;           //output cluster collection
+
+  TH2                   *fHistMatchEtaPhi[4];
+  TList                 *fOutputList;
+  TH2                   *fHistMatchEvsP[4];
+  TH2                   *fHistMatchdRvsEP[4];
+  TH1                   *fHistNclusvsCent;
+  TH1                   *fHistNclusMatchvsCent;
+  TH1                   *fHistEbefore;
+  TH1                   *fHistEafter;
+
+
+
+ private:
+  AliHadCorrTask(const AliHadCorrTask&);            // not implemented
+  AliHadCorrTask &operator=(const AliHadCorrTask&); // not implemented
+
+  ClassDef(AliHadCorrTask, 1) 
+};
+#endif
index e7849a446fd315fbb37a6966aa68346be7d43447..c206165f6d3838930538808ad4d8c1698df242d2 100644 (file)
@@ -4,5 +4,8 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
+#pragma link C++ class AliAnalysisTaskScale+;
 #pragma link C++ class AliEmcalJetTask+;
+#pragma link C++ class AliHadCorrTask+;
+
 #endif