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 {
 202   AliAnalysisTaskEmcalJet::ExecOnce();
 203
 204 //   AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
 205 //   if (fCaloName == "CaloClusters")
 206 //     am->LoadBranch("CaloClusters");
 207
 208   if (!fCalo2Name.IsNull() && !fCaloClusters2){
 209     fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 210     if (!fCaloClusters2){
 211       AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
 212       fInitialized = kFALSE;
 213       return;
 214     }
 215   }
 216
 217   if (!fMCParticlesName.IsNull() && !fMCParticles){
 218     fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
 219     if (!fMCParticles){
 220       AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
 221       fInitialized = kFALSE;
 222       return;
 223     }
 224   }
 225     //TString             fMCParticlesName;
 226     // TClonesArray       *fMCParticles;
 227     //  fCaloClusters2(),
 228  //    else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
 229 //       AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
 230 //       fCaloClusters2 = 0;
 231 //       fInitialized = kFALSE;
 232 //       return;
 233     //  }
 234 }
 235
 236
 237 //________________________________________________________________________
 238 Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
 239 {
 240   Int_t centbin = GetCentBin(fCent);
 241   //for pp analyses we will just use the first centrality bin
 242   if (GetBeamType()==0)
 243     centbin = 0;
 244   if (centbin>2)
 245     return kTRUE;
 246   if (!fTracks)
 247     return kTRUE;
 248   if (!fCaloClusters)
 249     return kTRUE;
 250   if (!fCaloClusters2)
 251     return kTRUE;
 252   const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
 253   const Int_t nTrack = fTracks->GetEntriesFast();
 254   //  const Int_t nMC = fMCParticles->GetEntriesFast();
 255
 256   TString fRhoScaledName = fRhoName;
 257   fRho = GetRhoFromEvent(fRhoScaledName);
 258   fRhoVal = fRho->GetVal();
 259   if (GetBeamType()==0)
 260     fRhoVal = 0;
 261   fHistRhovsCent->Fill(fCent,fRhoVal);
 262   const Int_t Njets = fJets->GetEntriesFast();
 263
 264   Int_t NjetAcc = 0;
 265
 266
 267   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 268     Int_t TrackMatch = 0;
 269      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 270      if (!jet)
 271        continue;
 272      if (jet->Area()==0)
 273        continue;
 274      if (jet->Pt()<0.1)
 275        continue;
 276      if (jet->MaxTrackPt()>100)
 277        continue;
 278      if (! AcceptJet(jet))
 279        continue;
 280      NjetAcc++;
 281      vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
 282      vector<Int_t> cluster_id2;
 283      Double_t Esub = 0; //total E subtracted from the jet
 284      Double_t jetPt = -500;
 285      jetPt = jet->Pt()-jet->Area()*fRhoVal;
 286      for (int i = 0;i<nTrack;i++){
 287        AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
 288        if (!emctrack)
 289          continue;
 290        if (!emctrack->IsEMCAL())
 291          continue;
 292        AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
 293        if (! track)
 294          continue;
 295        if (! AcceptTrack(track))
 296          continue;
 297        if (! IsJetTrack(jet,i,false))
 298          continue;
 299        Int_t iClus = track->GetEMCALcluster();
 300        if (iClus<0)
 301          continue;
 302        //we have the id of the matched cluster of a track constituent
 303        //       cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
 304        bool ischecked = false;
 305        for (UInt_t icid = 0;icid<cluster_id.size();icid++)
 306          if (cluster_id[icid] == iClus)
 307            ischecked = true; // we've already looked at this uncorrected cluster
 308        if (ischecked)
 309          continue; //no need to go further
 310        AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 311        AliVCluster *cluster = emcluster->GetCluster();
 312        if (! cluster)
 313          continue;
 314        if (! AcceptCluster(cluster))
 315          continue;
 316        Double_t etadiff = 999;
 317        Double_t phidiff = 999;
 318        AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
 319        if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
 320          continue;
 321        TrackMatch++; //this cluster has been matched, let's add it to the list
 322        cluster_id.push_back(iClus);
 323        Int_t ismatch = -1;
 324        //now we need to find its matched corrected cluster if any
 325        for (Int_t ic = 0;ic < nCluster2;ic++){
 326          //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
 327          AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
 328          if (! emcluster2)
 329            continue;
 330          if (! AcceptCluster(emcluster2))
 331            continue;
 332          if (! IsJetCluster(jet,ic,false))
 333            continue;
 334          TLorentzVector nPart;
 335          cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
 336          TLorentzVector nPart2;
 337          emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
 338          float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
 339          if (R < 0.001){//this cluster stayed in the jet!
 340            cluster_id2.push_back(ic);
 341            ismatch = ic;
 342          }
 343        } // end of cluster2 loop
 344        //we only get here if the track was matched to a cluster that hadn't been looked at before
 345        if (ismatch < 0) //this cluster was entirely deleted
 346          Esub+=cluster->E();
 347        else{ //get the corrected cluster
 348          AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
 349          Esub+=(cluster->E() - emcluster2temp->E());
 350        }
 351
 352      } // end of track loop
 353      fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
 354      fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
 355      fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
 356      fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
 357      fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
 358      fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
 359      fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
 360      if (jet->MaxTrackPt()<5.0)
 361        continue;
 362      fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
 363      fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
 364      fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
 365      fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
 366      fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
 367      fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
 368      fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
 369      if (centbin == 0)
 370        fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
 371   }
 372   fHistNjetvsCent->Fill(fCent,NjetAcc);
 373   return kTRUE;
 374 }
 375
 376
 377
 378
 379