PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeConst.cxx

   1 //
   2 // Do subtraction for jet shapes using constituents arXiv:1403.3108
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THnSparse.h>
  11 #include <TF1.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TChain.h>
  16 #include <TSystem.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19 #include <TTree.h>
  20
  21 #include "AliVCluster.h"
  22 #include "AliVTrack.h"
  23 #include "AliEmcalJet.h"
  24 #include "AliRhoParameter.h"
  25 #include "AliLog.h"
  26 #include "AliEmcalParticle.h"
  27 #include "AliMCEvent.h"
  28 #include "AliAODEvent.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliJetContainer.h"
  33 #include "AliParticleContainer.h"
  34
  35 #include "AliAnalysisTaskJetShapeConst.h"
  36
  37 ClassImp(AliAnalysisTaskJetShapeConst)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() :
  41   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
  42   fContainerBase(0),
  43   fContainerSub(1),
  44   fContainerNoEmb(2),
  45   fMinFractionShared(0),
  46   fSingleTrackEmb(kFALSE),
  47   fCreateTree(kFALSE),
  48   fJetMassVarType(kMass),
  49   fTreeJetBkg(),
  50   fJet1Vec(new TLorentzVector()),
  51   fJet2Vec(new TLorentzVector()),
  52   fJetSubVec(new TLorentzVector()),
  53   fArea(0),
  54   fAreaPhi(0),
  55   fAreaEta(0),
  56   fRho(0),
  57   fRhoM(0),
  58   fNConst(0),
  59   fMatch(0),
  60   fh2MSubMatch(0x0),
  61   fh2MSubPtRawAll(0x0),
  62   fh3MSubPtRawDRMatch(0x0),
  63   fh3MSubPtTrueDR(0x0),
  64   fh3MTruePtTrueDR(0x0),
  65   fh3PtTrueDeltaMDR(0x0),
  66   fh3PtTrueDeltaMRelDR(0x0),
  67   fhnMassResponse(0x0)
  68 {
  69   // Default constructor.
  70
  71   fh2MSubMatch         = new TH2F*[fNcentBins];
  72   fh2MSubPtRawAll      = new TH2F*[fNcentBins];
  73   fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
  74   fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
  75   fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
  76   fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
  77   fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
  78   fhnMassResponse      = new THnSparse*[fNcentBins];
  79
  80   for (Int_t i = 0; i < fNcentBins; i++) {
  81     fh2MSubMatch[i]         = 0;
  82     fh2MSubPtRawAll[i]      = 0;
  83     fh3MSubPtRawDRMatch[i]  = 0;
  84     fh3MSubPtTrueDR[i]      = 0;
  85     fh3MTruePtTrueDR[i]     = 0;
  86     fh3PtTrueDeltaMDR[i]    = 0;
  87     fh3PtTrueDeltaMRelDR[i] = 0;
  88     fhnMassResponse[i]      = 0;
  89   }
  90
  91   SetMakeGeneralHistograms(kTRUE);
  92   if(fCreateTree) DefineOutput(2, TTree::Class());
  93 }
  94
  95 //________________________________________________________________________
  96 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
  97   AliAnalysisTaskEmcalJet(name, kTRUE),
  98   fContainerBase(0),
  99   fContainerSub(1),
 100   fContainerNoEmb(2),
 101   fMinFractionShared(0),
 102   fSingleTrackEmb(kFALSE),
 103   fCreateTree(kFALSE),
 104   fJetMassVarType(kMass),
 105   fTreeJetBkg(0),
 106   fJet1Vec(new TLorentzVector()),
 107   fJet2Vec(new TLorentzVector()),
 108   fJetSubVec(new TLorentzVector()),
 109   fArea(0),
 110   fAreaPhi(0),
 111   fAreaEta(0),
 112   fRho(0),
 113   fRhoM(0),
 114   fNConst(0),
 115   fMatch(0),
 116   fh2MSubMatch(0x0),
 117   fh2MSubPtRawAll(0x0),
 118   fh3MSubPtRawDRMatch(0x0),
 119   fh3MSubPtTrueDR(0x0),
 120   fh3MTruePtTrueDR(0x0),
 121   fh3PtTrueDeltaMDR(0x0),
 122   fh3PtTrueDeltaMRelDR(0x0),
 123   fhnMassResponse(0x0)
 124 {
 125   // Standard constructor.
 126
 127   fh2MSubMatch         = new TH2F*[fNcentBins];
 128   fh2MSubPtRawAll      = new TH2F*[fNcentBins];
 129   fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
 130   fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
 131   fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
 132   fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
 133   fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
 134   fhnMassResponse      = new THnSparse*[fNcentBins];
 135
 136   for (Int_t i = 0; i < fNcentBins; i++) {
 137     fh2MSubMatch[i]         = 0;
 138     fh2MSubPtRawAll[i]      = 0;
 139     fh3MSubPtRawDRMatch[i]  = 0;
 140     fh3MSubPtTrueDR[i]      = 0;
 141     fh3MTruePtTrueDR[i]     = 0;
 142     fh3PtTrueDeltaMDR[i]    = 0;
 143     fh3PtTrueDeltaMRelDR[i] = 0;
 144     fhnMassResponse[i]      = 0;
 145   }
 146
 147   SetMakeGeneralHistograms(kTRUE);
 148   if(fCreateTree) DefineOutput(2, TTree::Class());
 149 }
 150
 151 //________________________________________________________________________
 152 AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
 153 {
 154   // Destructor.
 155 }
 156
 157 //________________________________________________________________________
 158 void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
 159 {
 160   // Create user output.
 161
 162   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 163
 164   Bool_t oldStatus = TH1::AddDirectoryStatus();
 165   TH1::AddDirectory(kFALSE);
 166
 167   const Int_t nBinsPt  = 200;
 168   const Double_t minPt = -50.;
 169   const Double_t maxPt = 150.;
 170
 171   Int_t nBinsM  = 100;
 172   Double_t minM = -20.;
 173   Double_t maxM = 80.;
 174   if(fJetMassVarType==kRatMPt) {
 175     nBinsM = 100;
 176     minM   = -0.2;
 177     maxM   = 0.8;
 178   }
 179
 180   Int_t nBinsDM  = 100;
 181   Double_t minDM = -50.;
 182   Double_t maxDM = 50.;
 183   if(fJetMassVarType==kRatMPt) {
 184     nBinsDM = 100;
 185     minDM   = -0.5;
 186     maxDM   = 0.5;
 187   }
 188
 189   const Int_t nBinsDRToLJ  = 20; //distance to leading jet in Pb-Pb only event
 190   const Double_t minDRToLJ = 0.;
 191   const Double_t maxDRToLJ = 1.;
 192
 193   //Binning for THnSparse
 194   const Int_t nBinsSparse0 = 5;
 195   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
 196   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt, minDRToLJ};
 197   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
 198
 199   TString histName = "";
 200   TString histTitle = "";
 201   TString varName = "#it{M}_{jet}";
 202   if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
 203
 204   for (Int_t i = 0; i < fNcentBins; i++) {
 205     histName = Form("fh2MSubMatch_%d",i);
 206     histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
 207     fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
 208     fOutput->Add(fh2MSubMatch[i]);
 209
 210     histName = Form("fh2MSubPtRawAll_%d",i);
 211     histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
 212     fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 213     fOutput->Add(fh2MSubPtRawAll[i]);
 214
 215     histName = Form("fh3MSubPtRawDRMatch_%d",i);
 216     histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
 217     fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 218     fOutput->Add(fh3MSubPtRawDRMatch[i]);
 219
 220     histName = Form("fh3MSubPtTrueDR_%d",i);
 221     histTitle = Form("fh3MSubPtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
 222     fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 223     fOutput->Add(fh3MSubPtTrueDR[i]);
 224
 225     histName = Form("fh3MTruePtTrueDR_%d",i);
 226     histTitle = Form("fh3MTruePtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
 227     fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 228     fOutput->Add(fh3MTruePtTrueDR[i]);
 229
 230     histName = Form("fh3PtTrueDeltaMDR_%d",i);
 231     histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
 232     fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 233     fOutput->Add(fh3PtTrueDeltaMDR[i]);
 234
 235     histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
 236     histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
 237     fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 238     fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
 239
 240     histName = Form("fhnMassResponse_%d",i);
 241     histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#Delta R_{LJ}",i,varName.Data(),varName.Data());
 242     fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 243     fOutput->Add(fhnMassResponse[i]);
 244   }
 245
 246   // =========== Switch on Sumw2 for all histos ===========
 247   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 248     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 249     if (h1){
 250       h1->Sumw2();
 251       continue;
 252     }
 253     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 254     if(hn)hn->Sumw2();
 255   }
 256
 257   TH1::AddDirectory(oldStatus);
 258
 259   // Create a tree.
 260   if(fCreateTree) {
 261     fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
 262     fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
 263     fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
 264     fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
 265     fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
 266     fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
 267     fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
 268     fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
 269     fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
 270     fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
 271   }
 272
 273   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 274   if(fCreateTree) PostData(2, fTreeJetBkg);
 275 }
 276
 277 //________________________________________________________________________
 278 Bool_t AliAnalysisTaskJetShapeConst::Run()
 279 {
 280   // Run analysis code here, if needed. It will be executed before FillHistograms().
 281
 282   return kTRUE;
 283 }
 284
 285 //________________________________________________________________________
 286 Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
 287 {
 288   // Fill histograms.
 289
 290   AliEmcalJet* jet1  = NULL; //AA jet
 291   AliEmcalJet *jet2  = NULL; //Embedded Pythia jet
 292   //  AliEmcalJet *jet1T = NULL; //tagged AA jet
 293   //  AliEmcalJet *jet2T = NULL; //tagged Pythia jet
 294   AliEmcalJet *jetS = NULL;  //subtracted jet
 295   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 296   AliJetContainer *jetContS = GetJetContainer(fContainerSub);
 297
 298   //Get leading jet in Pb-Pb event without embedded objects
 299   AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
 300   AliEmcalJet *jetL = NULL;
 301   if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
 302
 303   AliDebug(11,Form("NJets  Incl: %d  Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
 304   if(jetCont) {
 305     jetCont->ResetCurrentID();
 306     while((jet1 = jetCont->GetNextAcceptJet())) {
 307       jet2  = NULL;
 308       jetS  = NULL;
 309
 310       //Get constituent subtacted version of jet
 311       Int_t ifound = 0;
 312       Int_t ilab = -1;
 313       for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
 314         //if(ifound==1) continue;
 315         jetS = jetContS->GetJet(i);
 316         if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
 317           ifound++;
 318           if(ifound==1) ilab = i;
 319         }
 320       }
 321       if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
 322       if(ifound==0) jetS = 0x0;
 323       else jetS = jetContS->GetJet(ilab);
 324
 325       Double_t mjetS = jetS->M();
 326       Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
 327       Double_t var = mjetS;
 328       if(fJetMassVarType==kRatMPt) {
 329         if(ptjet1>0. || ptjet1<0.) var = mjetS/ptjet1;
 330         else var = -999.;
 331       }
 332
 333       //Fill histograms for all AA jets
 334       fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
 335
 336       Double_t fraction = 0.;
 337       fMatch = 0;
 338       fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
 339       if(fSingleTrackEmb) {
 340         AliVParticle *vp = GetEmbeddedConstituent(jet1);
 341         if(vp) {
 342           fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 343           fMatch = 1;
 344         }
 345       } else {
 346         jet2 = jet1->ClosestJet();
 347         fraction = jetCont->GetFractionSharedPt(jet1);
 348         fMatch = 1;
 349         if(fMinFractionShared>0.) {
 350           if(fraction>fMinFractionShared) {
 351             fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
 352             fMatch = 1;
 353           } else
 354             fMatch = 0;
 355         }
 356       }
 357
 358       //      if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
 359
 360       //Fill histograms for matched jets
 361       fh2MSubMatch[fCentBin]->Fill(var,fMatch);
 362       if(fMatch==1) {
 363         Double_t drToLJ = -1.;
 364         if(jetL) drToLJ = jet1->DeltaR(jetL);
 365         fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
 366         if(jet2) {
 367           Double_t var2 = -999.;
 368           if(jet2->Pt()>0. || jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
 369           fh3MSubPtTrueDR[fCentBin]->Fill(var,jet2->Pt(),drToLJ);
 370           fh3MTruePtTrueDR[fCentBin]->Fill(var2,jet2->Pt(),drToLJ);
 371           fh3PtTrueDeltaMDR[fCentBin]->Fill(jet2->Pt(),var-var2,drToLJ);
 372           if(jet2->M()>0.) fh3PtTrueDeltaMRelDR[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,drToLJ);
 373           Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),drToLJ};
 374           fhnMassResponse[fCentBin]->Fill(varsp);
 375         }
 376       }
 377
 378       if(fCreateTree) {
 379         fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
 380         if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
 381         else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
 382         fArea = (Float_t)jet1->Area();
 383         fAreaPhi = (Float_t)jet1->AreaPhi();
 384         fAreaEta = (Float_t)jet1->AreaEta();
 385         fRho  = (Float_t)jetCont->GetRhoVal();
 386         fRhoM = (Float_t)jetCont->GetRhoMassVal();
 387         fNConst = (Int_t)jet1->GetNumberOfTracks();
 388         fTreeJetBkg->Fill();
 389       }
 390     } //jet1 loop
 391   }//jetCont
 392
 393
 394   return kTRUE;
 395 }
 396
 397 //________________________________________________________________________
 398 AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
 399
 400   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 401   AliVParticle *vp = 0x0;
 402   AliVParticle *vpe = 0x0; //embedded particle
 403   Int_t nc = 0;
 404   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 405     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
 406     if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
 407     if(!vpe) vpe = vp;
 408     else if(vp->Pt()>vpe->Pt()) vpe = vp;
 409     nc++;
 410   }
 411   AliDebug(11,Form("Found %d embedded particles",nc));
 412   return vpe;
 413 }
 414
 415
 416 //________________________________________________________________________
 417 Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
 418   //
 419   // retrieve event objects
 420   //
 421
 422   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 423     return kFALSE;
 424
 425   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 426   jetCont->LoadRhoMass(InputEvent());
 427
 428   return kTRUE;
 429 }
 430
 431 //_______________________________________________________________________
 432 void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
 433 {
 434   // Called once at the end of the analysis.
 435 }
 436