update to master versions
[u/mrichter/AliRoot.git] / 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 = -25.;
182   Double_t maxDM = 25.;
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,400,-1.,3.,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