PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMassStructure.cxx

   1 //
   2 // Jet mass structure analysis task.
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <TProfile2D.h>
  11 #include <THnSparse.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 <TMath.h>
  20 #include <TRandom3.h>
  21
  22 #include "AliVCluster.h"
  23 #include "AliVTrack.h"
  24 #include "AliEmcalJet.h"
  25 #include "AliRhoParameter.h"
  26 #include "AliLog.h"
  27 #include "AliEmcalParticle.h"
  28 #include "AliMCEvent.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliMCEvent.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliJetContainer.h"
  34
  35 #include "AliAODEvent.h"
  36
  37 #include "AliAnalysisTaskEmcalJetMassStructure.h"
  38
  39 ClassImp(AliAnalysisTaskEmcalJetMassStructure)
  40
  41 //________________________________________________________________________
  42 AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure() :
  43   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassStructure", kTRUE),
  44   fContainerBase(0),
  45   fMinFractionShared(0),
  46   fJetMassType(kRaw),
  47   fRandom(0),
  48   fEfficiencyFixed(1.),
  49   fh3PtDRMass(0),
  50   fh3PtDRRho(0),
  51   fh3PtDRMassCorr(0),
  52   fh3PtDRRhoCorr(0),
  53   fh2PtMass(0),
  54   fh2PtMassCorr(0),
  55   fhnMassResponse(0),
  56   fhnMassResponseCorr(0),
  57   fh3JetPtDRTrackPt(0)
  58 {
  59   // Default constructor.
  60
  61   fh3PtDRMass                       = new TH3F*[fNcentBins];
  62   fh3PtDRRho                        = new TH3F*[fNcentBins];
  63   fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
  64   fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
  65   fh2PtMass                         = new TH2F*[fNcentBins];
  66   fh2PtMassCorr                     = new TH2F*[fNcentBins];
  67   fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];
  68
  69   for (Int_t i = 0; i < fNcentBins; i++) {
  70     fh3PtDRMass[i]       = 0;
  71     fh3PtDRRho[i]        = 0;
  72     fh3PtDRMassCorr[i]   = 0;
  73     fh3PtDRRhoCorr[i]    = 0;
  74     fh2PtMass[i]         = 0;
  75     fh2PtMassCorr[i]     = 0;
  76     fh3JetPtDRTrackPt[i] = 0;
  77   }
  78
  79   SetMakeGeneralHistograms(kTRUE);
  80 }
  81
  82 //________________________________________________________________________
  83 AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure(const char *name) :
  84   AliAnalysisTaskEmcalJet(name, kTRUE),
  85   fContainerBase(0),
  86   fMinFractionShared(0),
  87   fJetMassType(kRaw),
  88   fRandom(0),
  89   fEfficiencyFixed(1.),
  90   fh3PtDRMass(0),
  91   fh3PtDRRho(0),
  92   fh3PtDRMassCorr(0),
  93   fh3PtDRRhoCorr(0),
  94   fh2PtMass(0),
  95   fh2PtMassCorr(0),
  96   fhnMassResponse(0),
  97   fhnMassResponseCorr(0),
  98   fh3JetPtDRTrackPt(0)
  99 {
 100   // Standard constructor.
 101
 102   fh3PtDRMass                       = new TH3F*[fNcentBins];
 103   fh3PtDRRho                        = new TH3F*[fNcentBins];
 104   fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
 105   fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
 106   fh2PtMass                         = new TH2F*[fNcentBins];
 107   fh2PtMassCorr                     = new TH2F*[fNcentBins];
 108   fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];
 109
 110   for (Int_t i = 0; i < fNcentBins; i++) {
 111     fh3PtDRMass[i]       = 0;
 112     fh3PtDRRho[i]        = 0;
 113     fh3PtDRMassCorr[i]   = 0;
 114     fh3PtDRRhoCorr[i]    = 0;
 115     fh2PtMass[i]         = 0;
 116     fh2PtMassCorr[i]     = 0;
 117     fh3JetPtDRTrackPt[i] = 0;
 118   }
 119
 120   SetMakeGeneralHistograms(kTRUE);
 121 }
 122
 123 //________________________________________________________________________
 124 AliAnalysisTaskEmcalJetMassStructure::~AliAnalysisTaskEmcalJetMassStructure()
 125 {
 126   // Destructor.
 127
 128   delete fRandom;
 129 }
 130
 131 //________________________________________________________________________
 132 void AliAnalysisTaskEmcalJetMassStructure::UserCreateOutputObjects()
 133 {
 134   // Create user output.
 135
 136   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 137
 138   if(!fRandom) fRandom = new TRandom3(0);
 139
 140   Bool_t oldStatus = TH1::AddDirectoryStatus();
 141   TH1::AddDirectory(kFALSE);
 142
 143   const Int_t nBinsPt  = 200;
 144   const Double_t minPt = -50.;
 145   const Double_t maxPt = 150.;
 146   Double_t *binsPt = new Double_t[nBinsPt+1];
 147   for(Int_t i=0; i<=nBinsPt; i++) binsPt[i]=(Double_t)minPt + (maxPt-minPt)/nBinsPt*(Double_t)i ;
 148
 149   const Int_t nBinsM  = 120;
 150   const Double_t minM = -20.;
 151   const Double_t maxM = 40.;
 152   Double_t *binsM = new Double_t[nBinsM+1];
 153   for(Int_t i=0; i<=nBinsM; i++) binsM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;
 154
 155   const Int_t nBinsR  = 80;
 156   const Double_t minR = -0.005;
 157   const Double_t maxR = 0.795;
 158   Double_t *binsR = new Double_t[nBinsR+1];
 159   for(Int_t i=0; i<=nBinsR; i++) binsR[i]=(Double_t)minR + (maxR-minR)/nBinsR*(Double_t)i ;
 160
 161   //track pt
 162   Double_t binWidth1 = 0.1;
 163   Double_t binWidth2 = 1.;
 164   Double_t binWidth3 = 1.;
 165   const Float_t ptmin1 =  0.;
 166   const Float_t ptmax1 =  5.;
 167   const Float_t ptmin2 =  ptmax1;
 168   const Float_t ptmax2 =  10.;
 169   const Float_t ptmin3 =  ptmax2 ;
 170   const Float_t ptmax3 =  100.;
 171   const Int_t nbin11 = (int)((ptmax1-ptmin1)/binWidth1);
 172   const Int_t nbin12 = (int)((ptmax2-ptmin2)/binWidth2)+nbin11;
 173   const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
 174   Int_t nBinsPtTr=nbin13;
 175   //Create array with low edges of each bin
 176   Double_t *binsPtTr=new Double_t[nBinsPtTr+1];
 177   for(Int_t i=0; i<=nBinsPtTr; i++) {
 178     if(i<=nbin11) binsPtTr[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
 179     if(i<=nbin12 && i>nbin11) binsPtTr[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 180     if(i<=nbin13 && i>nbin12) binsPtTr[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 181   }
 182
 183   //Binning for THnSparse
 184   const Int_t nBinsSparse0 = 5;
 185   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
 186   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt};
 187   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt};
 188
 189   TString histName = "";
 190   TString histTitle = "";
 191   for (Int_t i = 0; i < fNcentBins; i++) {
 192     histName = TString::Format("fh3PtDRMass_%d",i);
 193     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
 194     fh3PtDRMass[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 195     fOutput->Add(fh3PtDRMass[i]);
 196
 197     histName = TString::Format("fh3PtDRRho_%d",i);
 198     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
 199     fh3PtDRRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 200     fOutput->Add(fh3PtDRRho[i]);
 201
 202     histName = TString::Format("fh3PtDRMassCorr_%d",i);
 203     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
 204     fh3PtDRMassCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 205     fOutput->Add(fh3PtDRMassCorr[i]);
 206
 207     histName = TString::Format("fh3PtDRRhoCorr_%d",i);
 208     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
 209     fh3PtDRRhoCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 210     fOutput->Add(fh3PtDRRhoCorr[i]);
 211
 212     histName = TString::Format("fh2PtMass_%d",i);
 213     histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
 214     fh2PtMass[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 215     fOutput->Add(fh2PtMass[i]);
 216
 217     histName = TString::Format("fh2PtMassCorr_%d",i);
 218     histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
 219     fh2PtMassCorr[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 220     fOutput->Add(fh2PtMassCorr[i]);
 221
 222     histName = TString::Format("fh3JetPtDRTrackPt_%d",i);
 223     histTitle = TString::Format("%s;#it{p}_{T,jet};r;#it{p}_{T,track}",histName.Data());
 224     fh3JetPtDRTrackPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,binsPt,nBinsR,binsR,nBinsPtTr,binsPtTr);
 225     fOutput->Add(fh3JetPtDRTrackPt[i]);
 226   }
 227
 228   histName = "fhnMassResponse";
 229   histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
 230   fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 231   fOutput->Add(fhnMassResponse);
 232
 233   histName = "fhnMassResponseCorr";
 234   histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
 235   fhnMassResponseCorr = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 236   fOutput->Add(fhnMassResponseCorr);
 237
 238   // =========== Switch on Sumw2 for all histos ===========
 239   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 240     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 241     if (h1){
 242       h1->Sumw2();
 243       continue;
 244     }
 245     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 246     if(hn)hn->Sumw2();
 247   }
 248
 249   TH1::AddDirectory(oldStatus);
 250
 251   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 252
 253   if(binsPt)                delete [] binsPt;
 254   if(binsPtTr)              delete [] binsPtTr;
 255   if(binsM)                 delete [] binsM;
 256   if(binsR)                 delete [] binsR;
 257
 258 }
 259
 260 //________________________________________________________________________
 261 Bool_t AliAnalysisTaskEmcalJetMassStructure::Run()
 262 {
 263   // Run analysis code here, if needed. It will be executed before FillHistograms().
 264
 265   return kTRUE;
 266 }
 267
 268 //________________________________________________________________________
 269 Bool_t AliAnalysisTaskEmcalJetMassStructure::FillHistograms()
 270 {
 271   // Fill histograms.
 272
 273   // study how jet mass builds up as function of radial distance to jet axis
 274
 275   AliEmcalJet* jet1 = NULL;
 276   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 277   if(jetCont) {
 278     jetCont->ResetCurrentID();
 279     while((jet1 = jetCont->GetNextAcceptJet())) {
 280
 281       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 282       Double_t mJet1 = GetJetMass(jet1);
 283
 284       // if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
 285       //        continue;
 286
 287       //Sort array based on distance to jet axis
 288       static Int_t indexes[9999] = {-1};
 289       static Float_t dr[9999] = {0};
 290       if(!GetSortedArray(indexes,dr,jet1)) continue;
 291
 292       AliVParticle *vp;
 293       TLorentzVector sumVec;      sumVec.SetPtEtaPhiM(0.,0.,0.,0.);
 294       TLorentzVector corrVec;     corrVec.SetPtEtaPhiM(0.,0.,0.,0.);
 295       TLorentzVector curVec;
 296       for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
 297         vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
 298         if(!vp) continue;
 299
 300         fh3JetPtDRTrackPt[fCentBin]->Fill(ptJet1,dr[indexes[i]],vp->Pt());
 301
 302         curVec.SetPtEtaPhiM(vp->Pt(),vp->Eta(),vp->Phi(),vp->M());
 303         sumVec+=curVec;
 304         corrVec+=curVec;
 305         //      Printf("%d  %f",i,dr[indexes[i]]);
 306         fh3PtDRMass[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.M()/mJet1);
 307         fh3PtDRRho[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.Pt()/ptJet1);
 308
 309         fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
 310         fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
 311
 312         Double_t eff = GetEfficiency(vp->Pt());
 313         Double_t rnd = fRandom->Uniform(1.);
 314         if(rnd>eff) {//put particle back at random position in annulus; using now zero width for annulus
 315           Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
 316           Double_t rr = dr[indexes[i]];
 317           rr = fRandom->Uniform(0.,jetCont->GetJetRadius());
 318           Double_t deta = rr*TMath::Cos(t);
 319           Double_t dphi = rr*TMath::Sin(t);
 320           curVec.SetPtEtaPhiM(vp->Pt(),deta+jet1->Eta(),dphi+jet1->Phi(),vp->M());
 321           corrVec+=curVec;
 322
 323           fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
 324           fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
 325         }
 326       }
 327
 328       fh2PtMass[fCentBin]->Fill(ptJet1,mJet1);
 329       fh2PtMassCorr[fCentBin]->Fill(corrVec.Pt(), corrVec.M());
 330
 331       //fill detector response
 332       AliEmcalJet *jPart = jet1->ClosestJet();
 333       if(jPart) {
 334         Double_t var[4] = {mJet1,jPart->M(),ptJet1,jPart->Pt()};
 335         fhnMassResponse->Fill(var);
 336
 337         Double_t varCorr[4] = {corrVec.M(),jPart->M(),corrVec.Pt(),jPart->Pt()};
 338         fhnMassResponseCorr->Fill(varCorr);
 339       }
 340     }
 341   }
 342
 343   return kTRUE;
 344 }
 345
 346 //________________________________________________________________________
 347 Double_t AliAnalysisTaskEmcalJetMassStructure::GetEfficiency(Double_t pt) {
 348   pt = 1.*pt;
 349   Double_t eff = 1.;
 350   if(fEfficiencyFixed<1.) return fEfficiencyFixed;
 351
 352   return eff;
 353 }
 354
 355 //________________________________________________________________________
 356 Bool_t AliAnalysisTaskEmcalJetMassStructure::GetSortedArray(Int_t indexes[], Float_t dr[], AliEmcalJet *jet) const
 357 {
 358   // sort array
 359
 360   const Int_t n = (Int_t)jet->GetNumberOfTracks();//(Int_t)array.size();
 361
 362   if (n < 1)
 363     return kFALSE;
 364
 365   AliVParticle *vp;
 366   for (Int_t i = 0; i < n; i++) {
 367     vp = static_cast<AliVParticle*>(jet->TrackAt(i, fTracks));
 368     if(!vp) continue;
 369     dr[i] = jet->DeltaR(vp);
 370   }
 371
 372   TMath::Sort(n, dr, indexes);
 373
 374   return kTRUE;
 375 }
 376
 377
 378 //________________________________________________________________________
 379 Double_t AliAnalysisTaskEmcalJetMassStructure::GetJetMass(AliEmcalJet *jet) {
 380   //calc subtracted jet mass
 381   if(fJetMassType==kRaw)
 382     return jet->M();
 383   else if(fJetMassType==kDeriv)
 384     return jet->GetSecondOrderSubtracted();
 385
 386   return 0;
 387 }
 388
 389 //________________________________________________________________________
 390 Bool_t AliAnalysisTaskEmcalJetMassStructure::RetrieveEventObjects() {
 391   //
 392   // retrieve event objects
 393   //
 394   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 395     return kFALSE;
 396
 397   return kTRUE;
 398 }
 399
 400 //_______________________________________________________________________
 401 void AliAnalysisTaskEmcalJetMassStructure::Terminate(Option_t *)
 402 {
 403   // Called once at the end of the analysis.
 404 }
 405