5ff619c4eed5dbd60662be2c8b369e4b16d661a9
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCore.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TMath.h"
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TH3F.h"
7 #include "THnSparse.h"
8 #include "TCanvas.h"
9
10 #include "AliLog.h"
11
12 #include "AliAnalysisTask.h"
13 #include "AliAnalysisManager.h"
14
15 #include "AliVEvent.h"
16 #include "AliESDEvent.h"
17 #include "AliESDInputHandler.h"
18 #include "AliCentrality.h"
19 #include "AliAnalysisHelperJetTasks.h"
20 #include "AliInputEventHandler.h"
21 #include "AliAODJetEventBackground.h"
22 #include "AliAnalysisTaskFastEmbedding.h"
23
24 #include "AliAODEvent.h"
25 #include "AliAODJet.h"
26
27 #include "AliAnalysisTaskJetCore.h"
28
29 ClassImp(AliAnalysisTaskJetCore)
30
31 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
32 AliAnalysisTaskSE(),
33 fESD(0x0),
34 fAOD(0x0),
35 fBackgroundBranch(""),
36 fIsPbPb(kTRUE),
37 fOfflineTrgMask(AliVEvent::kAny),
38 fMinContribVtx(1),
39 fVtxZMin(-8.),
40 fVtxZMax(8.),
41 fEvtClassMin(0),
42 fEvtClassMax(4),
43 fCentMin(0.),
44 fCentMax(100.),
45 fNInputTracksMin(0),
46 fNInputTracksMax(-1),
47 fJetEtaMin(-.5),
48 fJetEtaMax(.5),
49 fJetPtMin(20.),
50 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
51 fJetPtFractionMin(0.5),
52 fNMatchJets(4),
53 fMatchMaxDist(0.8),
54 fKeepJets(kFALSE),
55 fRadioFrac(0.2),
56 fMinDist(0.1),
57 fkNbranches(2),
58 fkEvtClasses(12),
59 fOutputList(0x0),
60 fbEvent(kTRUE),
61 fHistEvtSelection(0x0),
62 fHistJetSelection(0x0),
63 fh2JetSelection(0x0),
64 fh2JetCoreMethod1C10(0x0),
65 fh2JetCoreMethod2C10(0x0),
66 fh2JetCoreMethod3C10(0x0),
67 fh2JetCoreMethod1C30(0x0),
68 fh2JetCoreMethod2C30(0x0),
69 fh2JetCoreMethod3C30(0x0),
70 fh2JetCoreMethod1C60(0x0),
71 fh2JetCoreMethod2C60(0x0),
72 fh2JetCoreMethod3C60(0x0)
73
74  
75 {
76    // default Constructor
77
78    fJetBranchName[0] = "";
79    fJetBranchName[1] = "";
80
81    fListJets[0] = new TList;
82    fListJets[1] = new TList;
83 }
84
85 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
86 AliAnalysisTaskSE(name),
87 fESD(0x0),
88 fAOD(0x0),
89 fBackgroundBranch(""),
90 fIsPbPb(kTRUE),
91 fOfflineTrgMask(AliVEvent::kAny),
92 fMinContribVtx(1),
93 fVtxZMin(-8.),
94 fVtxZMax(8.),
95 fEvtClassMin(0),
96 fEvtClassMax(4),
97 fCentMin(0.),
98 fCentMax(100.),
99 fNInputTracksMin(0),
100 fNInputTracksMax(-1),
101 fJetEtaMin(-.5),
102 fJetEtaMax(.5),
103 fJetPtMin(20.),
104 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
105 fJetPtFractionMin(0.5),
106 fNMatchJets(4),
107 fMatchMaxDist(0.8),
108 fKeepJets(kFALSE),
109 fRadioFrac(0.2),
110 fMinDist(0.1),
111 fkNbranches(2),
112 fkEvtClasses(12),
113 fOutputList(0x0),
114 fbEvent(kTRUE),
115 fHistEvtSelection(0x0),
116 fHistJetSelection(0x0),
117 fh2JetSelection(0x0),
118 fh2JetCoreMethod1C10(0x0),
119 fh2JetCoreMethod2C10(0x0),
120 fh2JetCoreMethod3C10(0x0),
121 fh2JetCoreMethod1C30(0x0),
122 fh2JetCoreMethod2C30(0x0),
123 fh2JetCoreMethod3C30(0x0),
124 fh2JetCoreMethod1C60(0x0),
125 fh2JetCoreMethod2C60(0x0),
126 fh2JetCoreMethod3C60(0x0)
127
128
129  {
130    // Constructor
131
132    fJetBranchName[0] = "";
133    fJetBranchName[1] = "";
134
135    fListJets[0] = new TList;
136    fListJets[1] = new TList;
137
138    DefineOutput(1, TList::Class());
139 }
140
141 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
142 {
143    delete fListJets[0];
144    delete fListJets[1];
145 }
146
147 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
148 {
149    fJetBranchName[0] = branch1;
150    fJetBranchName[1] = branch2;
151 }
152
153 void AliAnalysisTaskJetCore::Init()
154 {
155
156    // check for jet branches
157    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
158       AliError("Jet branch name not set.");
159    }
160
161 }
162
163 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
164 {
165    // Create histograms
166    // Called once
167    OpenFile(1);
168    if(!fOutputList) fOutputList = new TList;
169    fOutputList->SetOwner(kTRUE);
170
171    Bool_t oldStatus = TH1::AddDirectoryStatus();
172    TH1::AddDirectory(kFALSE);
173
174
175    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
176    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
177    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
178    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
179    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
180    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
181    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
182
183    fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
184    fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
185    fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
186    fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
187    fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
188    fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
189    fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
190    fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
191    fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
192
193    fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
194    fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
195    fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
196    fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
197    fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
198    fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
199    fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
200    fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
201    fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
202
203
204    //   UInt_t entries = 0; // bit coded, see GetDimParams() below
205    //   UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
206
207     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
208     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
209     fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5); 
210     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
211     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
212     fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5); 
213     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
214     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
215     fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5); 
216
217
218    fOutputList->Add(fHistEvtSelection);
219    fOutputList->Add(fHistJetSelection);
220    fOutputList->Add(fh2JetSelection);
221   
222   
223       fOutputList->Add(fh2JetCoreMethod1C10);
224       fOutputList->Add(fh2JetCoreMethod2C10);
225       fOutputList->Add(fh2JetCoreMethod3C10);
226       fOutputList->Add(fh2JetCoreMethod1C30);
227       fOutputList->Add(fh2JetCoreMethod2C30);
228       fOutputList->Add(fh2JetCoreMethod3C30);
229       fOutputList->Add(fh2JetCoreMethod1C60);
230       fOutputList->Add(fh2JetCoreMethod2C60);
231       fOutputList->Add(fh2JetCoreMethod3C60);
232      
233      
234
235    // =========== Switch on Sumw2 for all histos ===========
236    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
237       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
238       if (h1){
239          h1->Sumw2();
240          continue;
241       }
242     
243    }
244    TH1::AddDirectory(oldStatus);
245
246    PostData(1, fOutputList);
247 }
248
249 void AliAnalysisTaskJetCore::UserExec(Option_t *)
250 {
251    
252
253    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
254       AliError("Jet branch name not set.");
255       return;
256    }
257
258    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
259    if (!fESD) {
260       AliError("ESD not available");
261       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
262    } else {
263       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
264    }
265    if (!fAOD) {
266       AliError("AOD not available");
267       return;
268    }
269
270    // -- event selection --
271    fHistEvtSelection->Fill(1); // number of events before event selection
272
273    // physics selection
274    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
275    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
276    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
277       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
278       fHistEvtSelection->Fill(2);
279       PostData(1, fOutputList);
280       return;
281    }
282
283    // vertex selection
284    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
285    Int_t nTracksPrim = primVtx->GetNContributors();
286    if ((nTracksPrim < fMinContribVtx) ||
287          (primVtx->GetZ() < fVtxZMin) ||
288          (primVtx->GetZ() > fVtxZMax) ){
289       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
290       fHistEvtSelection->Fill(3);
291       PostData(1, fOutputList);
292       return;
293    }
294
295    // event class selection (from jet helper task)
296    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
297    if(fDebug) Printf("Event class %d", eventClass);
298    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
299       fHistEvtSelection->Fill(4);
300       PostData(1, fOutputList);
301       return;
302    }
303
304    // centrality selection
305    AliCentrality *cent = 0x0;
306    Float_t centValue = 0.; 
307    if(fESD) cent = fESD->GetCentrality();
308    if(cent) centValue = cent->GetCentralityPercentile("V0M");
309    if(fDebug) printf("centrality: %f\n", centValue);
310    if (centValue < fCentMin || centValue > fCentMax){
311       fHistEvtSelection->Fill(4);
312       PostData(1, fOutputList);
313       return;
314    }
315
316
317    // multiplicity due to input tracks
318    Int_t nInputTracks = GetNInputTracks();
319
320    if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
321       fHistEvtSelection->Fill(5);
322       PostData(1, fOutputList);
323       return;
324    }
325
326    
327    fHistEvtSelection->Fill(0); // accepted events  
328    // -- end event selection --
329
330
331    // get background
332    AliAODJetEventBackground* externalBackground = 0;
333    if(!externalBackground&&fBackgroundBranch.Length()){
334       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
335       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
336    }
337    Float_t rho = 0;
338    if(externalBackground)rho = externalBackground->GetBackground(0);
339
340
341    // fetch jets
342    TClonesArray *aodJets[2];
343    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // 
344    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // 
345
346
347
348    for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
349       fListJets[iJetType]->Clear();
350       if (!aodJets[iJetType]) continue;
351
352       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
353       
354       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
355          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
356          if (jet) fListJets[iJetType]->Add(jet);
357       }
358       fListJets[iJetType]->Sort();
359    }
360    
361    Double_t etabig=0;
362    Double_t ptbig=0;
363    Double_t areabig=0;
364    Double_t phibig=0.;
365    Double_t etasmall=0;
366    Double_t ptsmall=0;
367    Double_t areasmall=0;
368    Double_t distr=0.;
369    Double_t phismall=0.;
370    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
371            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
372            etabig  = jetbig->Eta();
373            phibig  = jetbig->Phi();
374            ptbig   = jetbig->Pt();
375            if(ptbig==0) continue; 
376            areabig = jetbig->EffectiveAreaCharged();
377            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
378                    Double_t dismin=100.;
379                    Double_t ptmax=-10.; 
380                    Int_t index1=-1;
381                    Int_t index2=-1;
382                    Double_t fracin=0.; 
383            //compute sum of the pt of the tracks in a concentric cone
384            TRefArray *genTrackList = jetbig->GetRefTracks();
385            Int_t nTracksGenJet = genTrackList->GetEntriesFast();
386            AliAODTrack* genTrack;
387              for(Int_t ir=0; ir<nTracksGenJet; ++ir){
388              genTrack = (AliAODTrack*)(genTrackList->At(ir));
389              Float_t etr=genTrack->Eta();
390              Float_t phir=genTrack->Phi();
391              distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
392              distr=TMath::Sqrt(distr);
393              if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}    
394              if(centValue<10) fh2JetCoreMethod3C10->Fill(ptbig-rho*areabig,fracin/ptbig);
395              if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptbig-rho*areabig,fracin/ptbig);
396              if(centValue>60) fh2JetCoreMethod3C60->Fill(ptbig-rho*areabig,fracin/ptbig); 
397              //cout<<"method3  "<<fracin/ptbig<<" "<<ptbig-rho*areabig<<endl;
398              /////////////////////////////////////////////////////////////
399
400                 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
401                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
402                   etasmall  = jetsmall->Eta();
403                   phismall = jetsmall->Phi();
404                   ptsmall   = jetsmall->Pt();
405                   areasmall = jetsmall->EffectiveAreaCharged();
406
407                   if((etasmall<fJetEtaMin)||(etasmall>fJetEtaMax)) continue;
408                   
409                    Float_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
410                    tmpDeltaR=TMath::Sqrt(tmpDeltaR);
411
412                   if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
413                     index2=j;}  
414   
415                   if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
416                     index1=j;}}
417
418                 //method1:most concentric jet=core 
419                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
420         
421                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
422                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
423                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig); 
424
425                   //cout<<"method1  "<<jetmethod1->Pt()/ptbig<<" "<<ptbig<<endl;   
426                 }
427                 //method2:hardest contained jet=core   
428                 if(index2!=-1){ 
429                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
430                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
431                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
432                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig); 
433
434                   //cout<<"method2  "<<jetmethod2->Pt()/ptbig<<" "<<ptbig<<endl;
435                                }  
436
437    
438       
439         
440    }
441      
442       
443
444
445    PostData(1, fOutputList);
446 }
447
448 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
449 {
450    // Draw result to the screen
451    // Called once at the end of the query
452
453    if (!GetOutputData(1))
454    return;
455 }
456
457 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
458 {
459
460    Int_t nInputTracks = 0;
461
462    TString jbname(fJetBranchName[1]);
463    //needs complete event, use jets without background subtraction
464    for(Int_t i=1; i<=3; ++i){
465       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
466    }
467    // use only HI event
468    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
469    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
470
471    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
472    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
473    if(!tmpAODjets){
474       Printf("Jet branch %s not found", jbname.Data());
475       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
476       return -1;
477    }
478    
479    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
480       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
481       if(!jet) continue;
482       TRefArray *trackList = jet->GetRefTracks();
483       Int_t nTracks = trackList->GetEntriesFast();
484       nInputTracks += nTracks;
485       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
486    }
487    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
488
489    return nInputTracks;  
490 }
491
492
493
494
495