PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskDcalDijetPerf.cxx

   1 // $Id$
   2 //
   3 // Dcal dijet performance task
   4 //
   5 // Author: R. Reed
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliAODCaloCluster.h"
  17 #include "AliESDCaloCluster.h"
  18 #include "AliVTrack.h"
  19 #include "AliEmcalJet.h"
  20 #include "AliRhoParameter.h"
  21 #include "AliLog.h"
  22 #include "AliJetContainer.h"
  23 #include "AliParticleContainer.h"
  24 #include "AliClusterContainer.h"
  25 #include "AliPicoTrack.h"
  26
  27 #include "AliAnalysisTaskDcalDijetPerf.h"
  28
  29 ClassImp(AliAnalysisTaskDcalDijetPerf)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf() :
  33   AliAnalysisTaskEmcalJet("AliAnalysisTaskDcalDijetPerf", kTRUE),
  34   fHistTracksPt(0),
  35   fHistTracksEtaPhi(0),
  36   fHistClustersPt(0),
  37   fHistLeadingJetPt(0),
  38   fHistJetsPhiEta(0),
  39   fHistJetsPtArea(0),
  40   fHistJetsPtLeadHad(0),
  41   fHistJetsCorrPtArea(0),
  42   fHistJet1(0),
  43   fHistJet1m(0),
  44   fHistJet1nm(0),
  45   fHistJet2(0),
  46   fHistJet1to2(0),
  47   fJetsCont(0),
  48   fJetsCont2(0),
  49   fTracksCont(0),
  50   fCaloClustersCont(0)
  51
  52 {
  53   // Default constructor.
  54
  55   fHistTracksPt       = new TH1*[fNcentBins];
  56   fHistTracksEtaPhi   = new TH2*[fNcentBins];
  57   fHistClustersPt     = new TH1*[fNcentBins];
  58   fHistLeadingJetPt   = new TH1*[fNcentBins];
  59   fHistJetsPhiEta     = new TH2*[fNcentBins];
  60   fHistJetsPtArea     = new TH2*[fNcentBins];
  61   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  62   fHistJetsCorrPtArea = new TH2*[fNcentBins];
  63
  64   for (Int_t i = 0; i < fNcentBins; i++) {
  65     fHistTracksPt[i] = 0;
  66     fHistTracksEtaPhi[i] = 0;
  67     fHistClustersPt[i] = 0;
  68     fHistLeadingJetPt[i] = 0;
  69     fHistJetsPhiEta[i] = 0;
  70     fHistJetsPtArea[i] = 0;
  71     fHistJetsPtLeadHad[i] = 0;
  72     fHistJetsCorrPtArea[i] = 0;
  73   }
  74     fHistJet1 = 0;
  75     fHistJet1m = 0;
  76     fHistJet1nm = 0;
  77     fHistJet2 = 0;
  78     fHistJet1to2 = 0;
  79   SetMakeGeneralHistograms(kTRUE);
  80 }
  81
  82 //________________________________________________________________________
  83 AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf(const char *name) :
  84   AliAnalysisTaskEmcalJet(name, kTRUE),
  85   fHistTracksPt(0),
  86   fHistTracksEtaPhi(0),
  87   fHistClustersPt(0),
  88   fHistLeadingJetPt(0),
  89   fHistJetsPhiEta(0),
  90   fHistJetsPtArea(0),
  91   fHistJetsPtLeadHad(0),
  92   fHistJetsCorrPtArea(0),
  93   fHistJet1(0),
  94   fHistJet1m(0),
  95   fHistJet1nm(0),
  96   fHistJet2(0),
  97   fHistJet1to2(0),
  98   fJetsCont(0),
  99   fJetsCont2(0),
 100   fTracksCont(0),
 101   fCaloClustersCont(0)
 102 {
 103   // Standard constructor.
 104
 105   fHistTracksPt       = new TH1*[fNcentBins];
 106   fHistTracksEtaPhi   = new TH2*[fNcentBins];
 107   fHistClustersPt     = new TH1*[fNcentBins];
 108   fHistLeadingJetPt   = new TH1*[fNcentBins];
 109   fHistJetsPhiEta     = new TH2*[fNcentBins];
 110   fHistJetsPtArea     = new TH2*[fNcentBins];
 111   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
 112   fHistJetsCorrPtArea = new TH2*[fNcentBins];
 113
 114   for (Int_t i = 0; i < fNcentBins; i++) {
 115     fHistTracksPt[i] = 0;
 116     fHistTracksEtaPhi[i] = 0;
 117     fHistClustersPt[i] = 0;
 118     fHistLeadingJetPt[i] = 0;
 119     fHistJetsPhiEta[i] = 0;
 120     fHistJetsPtArea[i] = 0;
 121     fHistJetsPtLeadHad[i] = 0;
 122     fHistJetsCorrPtArea[i] = 0;
 123   }
 124   fHistJet1 = 0;
 125   fHistJet1m = 0;
 126   fHistJet1nm = 0;
 127   fHistJet2 = 0;
 128   fHistJet1to2 = 0;
 129
 130   SetMakeGeneralHistograms(kTRUE);
 131 }
 132
 133 //________________________________________________________________________
 134 AliAnalysisTaskDcalDijetPerf::~AliAnalysisTaskDcalDijetPerf()
 135 {
 136   // Destructor.
 137 }
 138
 139 //________________________________________________________________________
 140 void AliAnalysisTaskDcalDijetPerf::UserCreateOutputObjects()
 141 {
 142   // Create user output.
 143
 144   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 145
 146   fJetsCont           = GetJetContainer(0);
 147   fJetsCont2           = GetJetContainer(1);
 148   if(fJetsCont) { //get particles and clusters connected to jets
 149     fTracksCont       = fJetsCont->GetParticleContainer();
 150     fCaloClustersCont = fJetsCont->GetClusterContainer();
 151   } else {        //no jets, just analysis tracks and clusters
 152     fTracksCont       = GetParticleContainer(0);
 153     fCaloClustersCont = GetClusterContainer(0);
 154   }
 155   fTracksCont->SetClassName("AliVTrack");
 156   //fCaloClustersCont->SetClassName("AliVCluster");
 157   TString histname;
 158
 159   for (Int_t i = 0; i < fNcentBins; i++) {
 160     if (fParticleCollArray.GetEntriesFast()>0) {
 161       histname = "fHistTracksPt_";
 162       histname += i;
 163       fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 164       fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 165       fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
 166       fOutput->Add(fHistTracksPt[i]);
 167       histname = "fHistTracksEtaPhi_";
 168       histname += i;
 169         fHistTracksEtaPhi[i] = new TH2F(histname.Data(), histname.Data(), fNbins, -0.7, 0.7, fNbins, 0, TMath::Pi()*2);
 170       fHistTracksEtaPhi[i]->GetXaxis()->SetTitle("eta");
 171       fHistTracksEtaPhi[i]->GetYaxis()->SetTitle("phi");
 172       fOutput->Add(fHistTracksEtaPhi[i]);
 173     }
 174
 175     if (fClusterCollArray.GetEntriesFast()>0) {
 176       histname = "fHistClustersPt_";
 177       histname += i;
 178       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 179       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 180       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
 181       fOutput->Add(fHistClustersPt[i]);
 182     }
 183
 184     if (fJetCollArray.GetEntriesFast()>0) {
 185       histname = "fHistLeadingJetPt_";
 186       histname += i;
 187       fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 188       fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 189       fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 190       fOutput->Add(fHistLeadingJetPt[i]);
 191
 192       histname = "fHistJetsPhiEta_";
 193       histname += i;
 194       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
 195       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 196       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 197       fOutput->Add(fHistJetsPhiEta[i]);
 198
 199       histname = "fHistJetsPtArea_";
 200       histname += i;
 201       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
 202       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 203       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 204       fOutput->Add(fHistJetsPtArea[i]);
 205
 206       histname = "fHistJetsPtLeadHad_";
 207       histname += i;
 208       fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 209       fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 210       fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
 211       fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
 212       fOutput->Add(fHistJetsPtLeadHad[i]);
 213
 214       if (!(GetJetContainer()->GetRhoName().IsNull())) {
 215         histname = "fHistJetsCorrPtArea_";
 216         histname += i;
 217         fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 218         fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 219         fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 220         fOutput->Add(fHistJetsCorrPtArea[i]);
 221       }
 222     }
 223   }
 224
 225   Int_t nbins[] = {150,100,100,100,150};
 226   Double_t xmin[] = {  0,-0.7,             0, 0,   0};
 227   Double_t xmax[] = {150, 0.7,TMath::TwoPi(), 1, 150};
 228   fHistJet1 = new THnSparseF("Jets1Collection","Jets1Collection",5,nbins,xmin,xmax);
 229   fOutput->Add(fHistJet1);
 230   fHistJet1m = new THnSparseF("Jets1CollectionMatched","Jets1Collection",5,nbins,xmin,xmax);
 231   fOutput->Add(fHistJet1m);
 232   fHistJet1nm = new THnSparseF("Jets1CollectionNotMatched","Jets1Collection",5,nbins,xmin,xmax);
 233   fOutput->Add(fHistJet1nm);
 234   fHistJet2 = new THnSparseF("Jets2Collection","Jets2Collection",5,nbins,xmin,xmax);
 235   fOutput->Add(fHistJet2);
 236
 237   Int_t nbins2[] = {150,100,100,100,150,100,100,100,100};
 238   Double_t xmin2[] = {0,-0.7,0,0,0,-0.7,0,0,0};
 239   Double_t xmax2[] = {150,0.7,6.28,1,150,0.7,6.28,1,0.2};
 240   fHistJet1to2 = new THnSparseF("Jets1to2Collection","Jets1to2Collection",9,nbins2,xmin2,xmax2);
 241   fOutput->Add(fHistJet1to2);
 242
 243   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 244 }
 245
 246 //________________________________________________________________________
 247 Bool_t AliAnalysisTaskDcalDijetPerf::FillHistograms()
 248 {
 249   // Fill histograms.
 250
 251   if (fTracksCont) {
 252     AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 253     while(track) {
 254       fHistTracksPt[fCentBin]->Fill(track->Pt());
 255       fHistTracksEtaPhi[fCentBin]->Fill(track->Eta(),track->Phi());
 256       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 257     }
 258   }
 259
 260   if (fCaloClustersCont) {
 261     AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 262     while(cluster) {
 263       TLorentzVector nPart;
 264       cluster->GetMomentum(nPart, fVertex);
 265       fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 266
 267       cluster = fCaloClustersCont->GetNextAcceptCluster();
 268     }
 269   }
 270
 271     int N1 = 0;
 272   if (fJetsCont) {
 273     AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
 274     while(jet) {
 275         Float_t NEFpT = jet->Pt()*jet->NEF();
 276         N1++;
 277         //cout<<"loop 1 jet 1 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<" and NEFpT "<<NEFpT<<endl;
 278       fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 279       fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 280       Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
 281       fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
 282
 283         //110 degrees in the azimuthal angle |eta|<0.7
 284         //EMCal Tower 0.0143 x 0.0143
 285         //16x16 =  0.2288 x 0.2288 -> R = 0.1144
 286
 287         Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
 288         //cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
 289         fHistJet1->Fill(jetarray);
 290       if (fHistJetsCorrPtArea[fCentBin]) {
 291         Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
 292         fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
 293       }
 294       jet = fJetsCont->GetNextAcceptJet();
 295     }
 296
 297     jet = fJetsCont->GetLeadingJet();
 298     if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 299   }
 300    // cout<<"DONE LOOP 1"<<endl;
 301     int N2 = 0;
 302     if(fJetsCont2){
 303         //cout<<"We have a 2nd collection!!"<<endl;
 304         AliEmcalJet *jet = fJetsCont2->GetNextAcceptJet(0);
 305         while(jet){
 306             Float_t NEFpT = jet->Pt()*jet->NEF();
 307             //cout<<"loop 2 jet 2 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
 308             N2++;
 309             Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
 310             //cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
 311             fHistJet2->Fill(jetarray);
 312           //  cout<<" we have a jet! wiht pt = "<<jet->Pt()<<endl;
 313             jet = fJetsCont2->GetNextAcceptJet();
 314         }
 315     }
 316     //<<"DONE LOOP 2"<<endl;
 317
 318     //cout<<" There are "<<fJetsCont->GetNJets()<<" cont1 jets and "<<fJetsCont2->GetNJets()<<" cont2 jets"<<endl;
 319
 320     int N1N2 = 0;
 321     int N1N2m = 0;
 322    if (fJetsCont&&fJetsCont2) {
 323          AliEmcalJet *jet1 = fJetsCont->GetNextAcceptJet(0);
 324          AliEmcalJet *jet2 = fJetsCont2->GetNextAcceptJet(0);
 325          while(jet1){
 326              bool ismatched = false;
 327              Float_t NEFpT1 = jet1->Pt()*jet1->NEF();
 328              Double_t jetarray1[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),NEFpT1};
 329              while(jet2){
 330                  N1N2++;
 331                  //cout<<"jet 1 has label, pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet1->Pt()<<" "<<jet1->Eta()<<" "<<jet1->Phi()<<" "<<jet1->NEF()<<endl;
 332                  //<<"jet 2 has lable pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet2->Pt()<<" "<<jet2->Eta()<<" "<<jet2->Phi()<<" "<<jet2->NEF()<<endl;
 333                  Double_t deta = jet1->Eta()-jet2->Eta();
 334                  Double_t dphi = RelativePhi(jet1->Phi(),jet2->Phi());
 335                  Double_t deta2 = deta*deta;
 336                  Double_t dphi2 = dphi*dphi;
 337                  Double_t dR = pow(deta2+dphi2,0.5);
 338                  //cout<<"dR is "<<dR<<endl;
 339                  Double_t jetarray[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet2->Pt(),jet2->Eta(),jet2->Phi(),jet2->NEF(),dR};
 340                  if (dR<0.2){
 341                      N1N2m++;
 342                      fHistJet1to2->Fill(jetarray);
 343                      ismatched = true;
 344                  }
 345                  jet2 = fJetsCont2->GetNextAcceptJet();
 346              }
 347              if (ismatched)
 348                  fHistJet1m->Fill(jetarray1);
 349              else
 350                  fHistJet1nm->Fill(jetarray1);
 351              jet1 = fJetsCont->GetNextAcceptJet();
 352              jet2 = fJetsCont2->GetNextAcceptJet(0);
 353          }
 354      }
 355
 356   return kTRUE;
 357 }
 358
 359 //________________________________________________________________________
 360 Float_t AliAnalysisTaskDcalDijetPerf:: RelativePhi(Double_t mphi,Double_t vphi) const
 361 {
 362     if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 363     else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 364     if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 365     else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 366     double dphi = mphi-vphi;
 367     if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 368     else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 369
 370     return dphi;//dphi in [-Pi, Pi]
 371 }
 372
 373
 374 //________________________________________________________________________
 375 void AliAnalysisTaskDcalDijetPerf::ExecOnce() {
 376
 377   AliAnalysisTaskEmcalJet::ExecOnce();
 378
 379   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 380   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 381   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 382
 383 }
 384
 385 //________________________________________________________________________
 386 Bool_t AliAnalysisTaskDcalDijetPerf::Run()
 387 {
 388   // Run analysis code here, if needed. It will be executed before FillHistograms().
 389
 390   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 391 }
 392
 393 //________________________________________________________________________
 394 void AliAnalysisTaskDcalDijetPerf::Terminate(Option_t *)
 395 {
 396   // Called once at the end of the analysis.
 397 }