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 using std::vector;
  31
  32 ClassImp(AliAnalysisTaskEmcalJetHadCorQA)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() :
  36   AliAnalysisTaskEmcalJet("jethadcor",kFALSE),
  37   fCalo2Name(),
  38   fCaloClusters2(),
  39   fMCParticlesName(),
  40   fMCParticles(),
  41   fHistRhovsCent(0),
  42   fHistNjetvsCent(0),
  43   fHistNTMatchvsPtvsNtack0(0)
  44 {
  45   for (int i = 0;i<3;i++){
  46   fHistNEFvsPt[i] = 0;
  47   fHistNTMatchvsPt[i] = 0;
  48   fHistNCMatchvsPt[i] = 0;
  49   fHistHadCorvsPt[i] = 0;
  50   fHistNEFvsPtBias[i] = 0;
  51   fHistNconvsPt[i] = 0;
  52   fHistNtvsPt[i] = 0;
  53   fHistNcvsPt[i] = 0;
  54   fHistNTMatchvsPtBias[i] = 0;
  55   fHistNCMatchvsPtBias[i] = 0;
  56   fHistHadCorvsPtBias[i] = 0;
  57   fHistNconvsPtBias[i] = 0;
  58   fHistNtvsPtBias[i] = 0;
  59   fHistNcvsPtBias[i] = 0;
  60   }
  61   // Default constructor.
  62
  63   SetMakeGeneralHistograms(kTRUE);
  64 }
  65
  66 //________________________________________________________________________
  67 AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
  68   AliAnalysisTaskEmcalJet(name,kTRUE),
  69   fCalo2Name(),
  70   fCaloClusters2(),
  71   fMCParticlesName(),
  72   fMCParticles(),
  73   fHistRhovsCent(0),
  74   fHistNjetvsCent(0),
  75   fHistNTMatchvsPtvsNtack0(0)
  76  {
  77   for (int i = 0;i<3;i++){
  78   fHistNEFvsPt[i] = 0;
  79   fHistNTMatchvsPt[i] = 0;
  80   fHistNCMatchvsPt[i] = 0;
  81   fHistHadCorvsPt[i] = 0;
  82   fHistNEFvsPtBias[i] = 0;
  83   fHistNconvsPt[i] = 0;
  84   fHistNtvsPt[i] = 0;
  85   fHistNcvsPt[i] = 0;
  86   fHistNTMatchvsPtBias[i] = 0;
  87   fHistNCMatchvsPtBias[i] = 0;
  88   fHistHadCorvsPtBias[i] = 0;
  89   fHistNconvsPtBias[i] = 0;
  90   fHistNtvsPtBias[i] = 0;
  91   fHistNcvsPtBias[i] = 0;
  92   }
  93    SetMakeGeneralHistograms(kTRUE);
  94  }
  95
  96 //________________________________________________________________________
  97 void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
  98 {
  99   if (! fCreateHisto)
 100     return;
 101   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 102   fHistRhovsCent             = new TH2F("RhovsCent",              "RhovsCent",             100, 0.0, 100.0, 500, 0, 500);
 103   fHistNjetvsCent            = new TH2F("NjetvsCent",             "NjetvsCent",            100, 0.0, 100.0, 100, 0, 100);
 104
 105   for (int i = 0;i<3;i++){
 106     char name[200];
 107     TString title;
 108     sprintf(name,"NEFvsPt%i",i);
 109     fHistNEFvsPt[i]         = new TH2F(name, name, 100,0,1,500,-250,250);
 110     fOutput->Add(fHistNEFvsPt[i]);
 111     sprintf(name,"NTMatchvsPt%i",i);
 112     fHistNTMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);
 113     fOutput->Add(fHistNTMatchvsPt[i]);
 114     sprintf(name,"NCMatchvsPt%i",i);
 115     fHistNCMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);
 116     fOutput->Add(fHistNCMatchvsPt[i]);
 117     sprintf(name,"HadCorvsPt%i",i);
 118     fHistHadCorvsPt[i]      = new TH2F(name, name, 1000,0,500,500,-250,250);
 119     fOutput->Add(fHistHadCorvsPt[i]);
 120     sprintf(name,"NconvsPt%i",i);
 121     fHistNconvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 122     fOutput->Add(fHistNconvsPt[i]);
 123     sprintf(name,"NtvsPt%i",i);
 124     fHistNtvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 125     fOutput->Add(fHistNtvsPt[i]);
 126     sprintf(name,"NcvsPt%i",i);
 127     fHistNcvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 128     fOutput->Add(fHistNcvsPt[i]);
 129     sprintf(name,"NEFvsPtBias%i",i);
 130     fHistNEFvsPtBias[i]     = new TH2F(name, name, 100,0,1,500,-250,250);
 131     fOutput->Add(fHistNEFvsPtBias[i]);
 132     sprintf(name,"NTMatchvsPtBias%i",i);
 133     fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
 134     fOutput->Add(fHistNTMatchvsPtBias[i]);
 135     sprintf(name,"NCMatchvsPtBias%i",i);
 136     fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
 137     fOutput->Add(fHistNCMatchvsPtBias[i]);
 138     sprintf(name,"HadCorvsPtBias%i",i);
 139     fHistHadCorvsPtBias[i]  = new TH2F(name, name, 1000,0,500,500,-250,250);
 140     fOutput->Add(fHistHadCorvsPtBias[i]);
 141    sprintf(name,"NconvsPtBias%i",i);
 142     fHistNconvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 143     fOutput->Add(fHistNconvsPtBias[i]);
 144     sprintf(name,"NtvsPtBias%i",i);
 145     fHistNtvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 146     fOutput->Add(fHistNtvsPtBias[i]);
 147     sprintf(name,"NcvsPtBias%i",i);
 148     fHistNcvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
 149     fOutput->Add(fHistNcvsPtBias[i]);
 150   }
 151   fHistNTMatchvsPtvsNtack0   = new TH3F("NTMmatchvsPtvsNtrack0",  "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
 152
 153   // for (Int_t i = 0;i<6;++i){
 154   //   name = TString(Form("JetPtvsTrackPt_%i",i));
 155   //   title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
 156   //   fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
 157   //   fOutput->Add(fHistJetPtvsTrackPt[i]);
 158
 159   // }
 160   fOutput->Add(fHistRhovsCent);
 161   fOutput->Add(fHistNjetvsCent);
 162   fOutput->Add(fHistNTMatchvsPtvsNtack0);
 163    PostData(1, fOutput);
 164 }
 165
 166 //________________________________________________________________________
 167
 168 Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const
 169 {
 170   // Get centrality bin.
 171
 172   Int_t centbin = -1;
 173   if (cent>=0 && cent<10)
 174     centbin = 0;
 175   else if (cent>=10 && cent<30)
 176     centbin = 1;
 177   else if (cent>=30 && cent<50)
 178     centbin = 2;
 179   else if (cent>50)
 180     centbin =3;
 181   return centbin;
 182 }
 183
 184 //________________________________________________________________________
 185
 186 Float_t AliAnalysisTaskEmcalJetHadCorQA:: RelativePhi(Double_t mphi,Double_t vphi) const
 187 {
 188   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 189   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 190   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 191   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 192   double dphi = mphi-vphi;
 193   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 194   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 195
 196   return dphi;//dphi in [-Pi, Pi]
 197 }
 198
 199 //________________________________________________________________________
 200 void AliAnalysisTaskEmcalJetHadCorQA::ExecOnce(){
 201   AliAnalysisTaskEmcalJet::ExecOnce();
 202
 203 //   AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
 204 //   if (fCaloName == "CaloClusters")
 205 //     am->LoadBranch("CaloClusters");
 206
 207   if (!fCalo2Name.IsNull() && !fCaloClusters2){
 208     fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 209     if (!fCaloClusters2){
 210       AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
 211       fInitialized = kFALSE;
 212       return;
 213     }
 214   }
 215
 216   if (!fMCParticlesName.IsNull() && !fMCParticles){
 217     fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
 218     if (!fMCParticles){
 219       AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
 220       fInitialized = kFALSE;
 221       return;
 222     }
 223   }
 224     //TString             fMCParticlesName;
 225     // TClonesArray       *fMCParticles;
 226     //  fCaloClusters2(),
 227  //    else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
 228 //       AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
 229 //       fCaloClusters2 = 0;
 230 //       fInitialized = kFALSE;
 231 //       return;
 232     //  }
 233 }
 234
 235
 236 //________________________________________________________________________
 237 Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
 238 {
 239   Int_t centbin = GetCentBin(fCent);
 240   //for pp analyses we will just use the first centrality bin
 241   if (GetBeamType()==0)
 242     centbin = 0;
 243   if (centbin>2)
 244     return kTRUE;
 245   if (!fTracks)
 246     return kTRUE;
 247   if (!fCaloClusters)
 248     return kTRUE;
 249   if (!fCaloClusters2)
 250     return kTRUE;
 251   const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
 252   const Int_t nTrack = fTracks->GetEntriesFast();
 253   //  const Int_t nMC = fMCParticles->GetEntriesFast();
 254
 255   TString fRhoScaledName = fRhoName;
 256   fRho = GetRhoFromEvent(fRhoScaledName);
 257   fRhoVal = fRho->GetVal();
 258   if (GetBeamType()==0)
 259     fRhoVal = 0;
 260   fHistRhovsCent->Fill(fCent,fRhoVal);
 261   const Int_t Njets = fJets->GetEntriesFast();
 262
 263   Int_t NjetAcc = 0;
 264
 265
 266   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 267     Int_t TrackMatch = 0;
 268      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 269      if (!jet)
 270        continue;
 271      if (jet->Area()==0)
 272        continue;
 273      if (jet->Pt()<0.1)
 274        continue;
 275      if (jet->MaxTrackPt()>100)
 276        continue;
 277      if (! AcceptJet(jet))
 278        continue;
 279      NjetAcc++;
 280      vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
 281      vector<Int_t> cluster_id2;
 282      Double_t Esub = 0; //total E subtracted from the jet
 283      Double_t jetPt = -500;
 284      jetPt = jet->Pt()-jet->Area()*fRhoVal;
 285      for (int i = 0;i<nTrack;i++){
 286        AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
 287        if (!emctrack)
 288          continue;
 289        if (!emctrack->IsEMCAL())
 290          continue;
 291        AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
 292        if (! track)
 293          continue;
 294        if (! AcceptTrack(track))
 295          continue;
 296        if (! IsJetTrack(jet,i,false))
 297          continue;
 298        Int_t iClus = track->GetEMCALcluster();
 299        if (iClus<0)
 300          continue;
 301        //we have the id of the matched cluster of a track constituent
 302        //       cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
 303        bool ischecked = false;
 304        for (UInt_t icid = 0;icid<cluster_id.size();icid++)
 305          if (cluster_id[icid] == iClus)
 306            ischecked = true; // we've already looked at this uncorrected cluster
 307        if (ischecked)
 308          continue; //no need to go further
 309        AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 310        AliVCluster *cluster = emcluster->GetCluster();
 311        if (! cluster)
 312          continue;
 313        if (! AcceptCluster(cluster))
 314          continue;
 315        Double_t etadiff = 999;
 316        Double_t phidiff = 999;
 317        AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
 318        if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
 319          continue;
 320        TrackMatch++; //this cluster has been matched, let's add it to the list
 321        cluster_id.push_back(iClus);
 322        Int_t ismatch = -1;
 323        //now we need to find its matched corrected cluster if any
 324        for (Int_t ic = 0;ic < nCluster2;ic++){
 325          //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
 326          AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
 327          if (! emcluster2)
 328            continue;
 329          if (! AcceptCluster(emcluster2))
 330            continue;
 331          if (! IsJetCluster(jet,ic,false))
 332            continue;
 333          TLorentzVector nPart;
 334          cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
 335          TLorentzVector nPart2;
 336          emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
 337          float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
 338          if (R < 0.001){//this cluster stayed in the jet!
 339            cluster_id2.push_back(ic);
 340            ismatch = ic;
 341          }
 342        } // end of cluster2 loop
 343        //we only get here if the track was matched to a cluster that hadn't been looked at before
 344        if (ismatch < 0) //this cluster was entirely deleted
 345          Esub+=cluster->E();
 346        else{ //get the corrected cluster
 347          AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
 348          Esub+=(cluster->E() - emcluster2temp->E());
 349        }
 350
 351      } // end of track loop
 352      fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
 353      fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
 354      fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
 355      fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
 356      fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
 357      fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
 358      fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
 359      if (jet->MaxTrackPt()<5.0)
 360        continue;
 361      fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
 362      fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
 363      fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
 364      fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
 365      fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
 366      fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
 367      fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
 368      if (centbin == 0)
 369        fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
 370   }
 371   fHistNjetvsCent->Fill(fCent,NjetAcc);
 372   return kTRUE;
 373 }
 374
 375
 376
 377
 378