PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadCorQA.cxx

   1 #include "AliAnalysisTaskEmcalJetHadCorQA.h"
   2
   3 #include <TCanvas.h>
   4 #include <TChain.h>
   5 #include <TClonesArray.h>
   6 #include <TH1F.h>
   7 #include <TH2F.h>
   8 #include <TH3F.h>
   9 #include <THnSparse.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12 #include <TParameter.h>
  13 #include <TParticle.h>
  14 #include <TTree.h>
  15 #include <TVector3.h>
  16
  17 #include "AliAODEvent.h"
  18 #include "AliAnalysisManager.h"
  19 #include "AliAnalysisTask.h"
  20 #include "AliCentrality.h"
  21 #include "AliESDEvent.h"
  22 #include "AliESDInputHandler.h"
  23 #include "AliEmcalJet.h"
  24 #include "AliVCluster.h"
  25 #include "AliVTrack.h"
  26 #include "AliRhoParameter.h"
  27 #include "AliEmcalParticle.h"
  28 #include "AliPicoTrack.h"
  29 #include "AliEMCALGeometry.h"
  30
  31 ClassImp(AliAnalysisTaskEmcalJetHadCorQA)
  32
  33 //________________________________________________________________________
  34 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() :
  35   AliAnalysisTaskEmcalJet("spectra",kFALSE),
  36   fCalo2Name(),
  37   fCaloClusters2(),
  38   fHistRhovsCent(0),
  39   fHistNjetvsCent(0)
  40 {
  41   for (int i = 0;i<3;i++){
  42   fHistNEFvsPt[i] = 0;
  43   fHistNTMatchvsPt[i] = 0;
  44   fHistNCMatchvsPt[i] = 0;
  45   fHistHadCorvsPt[i] = 0;
  46   fHistNEFvsPtBias[i] = 0;
  47   fHistNTMatchvsPtBias[i] = 0;
  48   fHistNCMatchvsPtBias[i] = 0;
  49   fHistHadCorvsPtBias[i] = 0;
  50   }
  51   // Default constructor.
  52
  53   SetMakeGeneralHistograms(kTRUE);
  54 }
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
  58   AliAnalysisTaskEmcalJet(name,kTRUE),
  59   fCalo2Name(),
  60   fCaloClusters2(),
  61   fHistRhovsCent(0),
  62   fHistNjetvsCent(0)
  63  {
  64   for (int i = 0;i<3;i++){
  65   fHistNEFvsPt[i] = 0;
  66   fHistNTMatchvsPt[i] = 0;
  67   fHistNCMatchvsPt[i] = 0;
  68   fHistHadCorvsPt[i] = 0;
  69   fHistNEFvsPtBias[i] = 0;
  70   fHistNTMatchvsPtBias[i] = 0;
  71   fHistNCMatchvsPtBias[i] = 0;
  72   fHistHadCorvsPtBias[i] = 0;
  73   }
  74    SetMakeGeneralHistograms(kTRUE);
  75  }
  76
  77 //________________________________________________________________________
  78 void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
  79 {
  80   if (! fCreateHisto)
  81     return;
  82   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  83   fHistRhovsCent             = new TH2F("RhovsCent",              "RhovsCent",             100, 0.0, 100.0, 500, 0, 500);
  84   fHistNjetvsCent            = new TH2F("NjetvsCent",             "NjetvsCent",            100, 0.0, 100.0, 100, 0, 100);
  85
  86   for (int i = 0;i<3;i++){
  87     char name[200];
  88     TString title;
  89     sprintf(name,"NEFvsPt%i",i);
  90     fHistNEFvsPt[i]         = new TH2F(name, name, 100,0,1,500,-250,250);
  91     fOutput->Add(fHistNEFvsPt[i]);
  92     sprintf(name,"NTMatchvsPt%i",i);
  93     fHistNTMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);
  94     fOutput->Add(fHistNTMatchvsPt[i]);
  95     sprintf(name,"NCMatchvsPt%i",i);
  96     fHistNCMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);
  97     fOutput->Add(fHistNCMatchvsPt[i]);
  98     sprintf(name,"HadCorvsPt%i",i);
  99     fHistHadCorvsPt[i]      = new TH2F(name, name, 1000,0,500,500,-250,250);
 100     fOutput->Add(fHistHadCorvsPt[i]);
 101     sprintf(name,"NEFvsPtBias%i",i);
 102     fHistNEFvsPtBias[i]     = new TH2F(name, name, 100,0,1,500,-250,250);
 103     fOutput->Add(fHistNEFvsPtBias[i]);
 104     sprintf(name,"NTMatchvsPtBias%i",i);
 105     fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
 106     fOutput->Add(fHistNTMatchvsPtBias[i]);
 107     sprintf(name,"NCMatchvsPtBias%i",i);
 108     fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
 109     fOutput->Add(fHistNCMatchvsPtBias[i]);
 110     sprintf(name,"HadCorvsPtBias%i",i);
 111     fHistHadCorvsPtBias[i]  = new TH2F(name, name, 1000,0,500,500,-250,250);
 112     fOutput->Add(fHistHadCorvsPtBias[i]);
 113   }
 114   fHistNTMatchvsPtvsNtack0   = new TH3F("NTMmatchvsPtvsNtrack0",  "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
 115
 116   // for (Int_t i = 0;i<6;++i){
 117   //   name = TString(Form("JetPtvsTrackPt_%i",i));
 118   //   title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
 119   //   fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
 120   //   fOutput->Add(fHistJetPtvsTrackPt[i]);
 121
 122   // }
 123   fOutput->Add(fHistRhovsCent);
 124   fOutput->Add(fHistNjetvsCent);
 125   fOutput->Add(fHistNTMatchvsPtvsNtack0);
 126    PostData(1, fOutput);
 127 }
 128
 129 //________________________________________________________________________
 130
 131 Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const
 132 {
 133   // Get centrality bin.
 134
 135   Int_t centbin = -1;
 136   if (cent>=0 && cent<10)
 137     centbin = 0;
 138   else if (cent>=10 && cent<30)
 139     centbin = 1;
 140   else if (cent>=30 && cent<50)
 141     centbin = 2;
 142   else if (cent>50)
 143     centbin =3;
 144   return centbin;
 145 }
 146
 147 //________________________________________________________________________
 148
 149 Float_t AliAnalysisTaskEmcalJetHadCorQA:: RelativePhi(Double_t mphi,Double_t vphi) const
 150 {
 151   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 152   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 153   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 154   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 155   double dphi = mphi-vphi;
 156   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 157   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 158
 159   return dphi;//dphi in [-Pi, Pi]
 160 }
 161
 162 //________________________________________________________________________
 163 void AliAnalysisTaskEmcalJetHadCorQA::ExecOnce(){
 164   AliAnalysisTaskEmcalJet::ExecOnce();
 165
 166 //   AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
 167 //   if (fCaloName == "CaloClusters")
 168 //     am->LoadBranch("CaloClusters");
 169
 170   if (!fCalo2Name.IsNull() && !fCaloClusters2){
 171     fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 172     if (!fCaloClusters2){
 173       AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
 174       fInitialized = kFALSE;
 175       return;
 176     }
 177     //  fCaloClusters2(),
 178  //    else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
 179 //       AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
 180 //       fCaloClusters2 = 0;
 181 //       fInitialized = kFALSE;
 182 //       return;
 183     //  }
 184   }
 185 }
 186
 187
 188 //________________________________________________________________________
 189 Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
 190 {
 191   Int_t centbin = GetCentBin(fCent);
 192   //for pp analyses we will just use the first centrality bin
 193   if (centbin == -1)
 194     centbin = 0;
 195   if (centbin>2)
 196     return kTRUE;
 197   if (!fTracks)
 198     return kTRUE;
 199
 200   const Int_t nCluster = fCaloClusters->GetEntriesFast();
 201   const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
 202   const Int_t nTrack = fTracks->GetEntriesFast();
 203
 204   TString fRhoScaledName = fRhoName;
 205   fRho = GetRhoFromEvent(fRhoScaledName);
 206   fRhoVal = fRho->GetVal();
 207   fHistRhovsCent->Fill(fCent,fRhoVal);
 208   const Int_t Njets = fJets->GetEntriesFast();
 209
 210   Int_t NjetAcc = 0;
 211
 212
 213   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 214     Int_t TrackMatch = 0;
 215      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 216      if (!jet)
 217        continue;
 218      if (jet->Area()==0)
 219        continue;
 220      if (jet->Pt()<0.1)
 221        continue;
 222      if (jet->MaxTrackPt()>100)
 223        continue;
 224      if (! AcceptJet(jet))
 225        continue;
 226      NjetAcc++;
 227      vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
 228      vector<Int_t> cluster_id2;
 229      Double_t Esub = 0; //total E subtracted from the jet
 230      Double_t jetPt = -500;
 231      jetPt = jet->Pt()-jet->Area()*fRhoVal;
 232      for (int i = 0;i<nTrack;i++){
 233        AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
 234        if (!emctrack)
 235          continue;
 236        if (!emctrack->IsEMCAL())
 237          continue;
 238        AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
 239        if (! track)
 240          continue;
 241        if (! AcceptTrack(track))
 242          continue;
 243        if (! IsJetTrack(jet,i,false))
 244          continue;
 245        Int_t iClus = track->GetEMCALcluster();
 246        if (iClus<0)
 247          continue;
 248        bool ischecked = false;
 249        for (Int_t icid = 0;i<cluster_id.size();icid++)
 250          if (cluster_id[icid] == iClus)
 251            ischecked = true; // we've already looked at this uncorrected cluster
 252        if (ischecked)
 253          continue; //no need to go further
 254        AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 255        AliVCluster *cluster = emcluster->GetCluster();
 256        if (! cluster)
 257          continue;
 258        if (! AcceptCluster(cluster))
 259          continue;
 260        Double_t etadiff = 999;
 261        Double_t phidiff = 999;
 262        AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
 263        if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
 264          continue;
 265        TrackMatch++; //this cluster has been matched, let's add it to the list
 266        cluster_id.push_back(iClus);
 267        Int_t ismatch = -1;
 268        //now we need to find its matched corrected cluster if any
 269        for (Int_t ic = 0;ic < nCluster2;ic++){
 270          //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
 271          AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
 272          if (! emcluster2)
 273            continue;
 274          if (! AcceptCluster(emcluster2))
 275            continue;
 276          if (! IsJetCluster(jet,ic,false))
 277            continue;
 278          TLorentzVector nPart;
 279          cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
 280          TLorentzVector nPart2;
 281          emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
 282          float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
 283          if (R < 0.001){//this cluster stayed in the jet!
 284            cluster_id2.push_back(ic);
 285            ismatch = ic;
 286          }
 287        } // end of cluster2 loop
 288        //we only get here if the track was matched to a cluster that hadn't been looked at before
 289        if (ismatch < 0) //this cluster was entirely deleted
 290          Esub+=cluster->E();
 291        else{ //get the corrected cluster
 292          AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
 293          Esub+=(cluster->E() - emcluster2temp->E());
 294        }
 295
 296      } // end of track loop
 297      fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
 298      fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
 299      //        fHistNCMatchvsPtvsCent(0),
 300      fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
 301      if (jet->MaxTrackPt()<5.0)
 302        continue;
 303      fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
 304      fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
 305      fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
 306      if (centbin == 0)
 307        fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,nTrack);
 308   }
 309   fHistNjetvsCent->Fill(fCent,NjetAcc);
 310   return kTRUE;
 311 }
 312
 313
 314
 315
 316