| 2cc66a40 |
1 | #ifndef ALIANALYSISTASKSE_H |
| 2 | |
| 3 | #include <Riostream.h> |
| 4 | #include <TROOT.h> |
| 5 | #include <TFile.h> |
| 6 | #include <TChain.h> |
| 7 | #include <TTree.h> |
| 8 | #include <TKey.h> |
| 9 | #include <TProfile.h> |
| 10 | #include <TProfile2D.h> |
| 11 | #include <TH1.h> |
| 12 | #include <TH1F.h> |
| 13 | #include <TH2F.h> |
| 14 | #include <TH1D.h> |
| 15 | #include <TH1I.h> |
| b4e22a71 |
16 | #include <TArrayF.h> |
| 2cc66a40 |
17 | #include <THnSparse.h> |
| 18 | #include <TCanvas.h> |
| 19 | #include <TList.h> |
| 20 | #include <TClonesArray.h> |
| 21 | #include <TObject.h> |
| 22 | #include <TMath.h> |
| 23 | #include <TSystem.h> |
| 24 | #include <TInterpreter.h> |
| 25 | #include "AliAnalysisTask.h" |
| 26 | #include "AliCentrality.h" |
| 27 | #include "AliStack.h" |
| 28 | #include "AliESDEvent.h" |
| 29 | #include "AliESDInputHandler.h" |
| 30 | #include "AliAODEvent.h" |
| 31 | #include "AliAODHandler.h" |
| 32 | #include "AliAnalysisManager.h" |
| 33 | #include "AliAnalysisTaskSE.h" |
| c8ab9633 |
34 | #include "AliAnalysisHelperJetTasks.h" |
| 35 | #include "AliInputEventHandler.h" |
| 2cc66a40 |
36 | #endif |
| 37 | |
| 38 | #include <time.h> |
| 39 | #include <TRandom3.h> |
| 9d6abaad |
40 | #include "AliGenEventHeader.h" |
| 2cc66a40 |
41 | #include "AliGenPythiaEventHeader.h" |
| 9d6abaad |
42 | #include "AliGenHijingEventHeader.h" |
| 2cc66a40 |
43 | #include "AliAODMCHeader.h" |
| 44 | #include "AliMCEvent.h" |
| 45 | #include "AliLog.h" |
| 46 | #include <AliEmcalJet.h> |
| 47 | #include <AliPicoTrack.h> |
| 48 | #include "AliVEventHandler.h" |
| 49 | #include "AliVParticle.h" |
| 50 | #include "AliAODMCParticle.h" |
| 51 | #include "AliAnalysisUtils.h" |
| 52 | #include "AliRhoParameter.h" |
| 53 | #include "TVector3.h" |
| b4e22a71 |
54 | #include "AliVVertex.h" |
| 2cc66a40 |
55 | |
| 56 | #include <stdio.h> |
| 57 | #include <stdlib.h> |
| 58 | #include <vector> |
| 59 | |
| 60 | #include "AliAnalysisTaskHJetSpectra.h" |
| 61 | using std::min; |
| 62 | |
| 63 | // ANALYSIS OF HIGH PT HADRON TRIGGER ASSOCIATED SPECTRUM OF RECOIL JETS IN P+PB |
| 64 | // Author Filip Krizek (17.May. 2014) |
| 65 | |
| 66 | //TODO: Not accessing the particles when using MC |
| 67 | //TODO: FillHistogram can be done better with virtual TH1(?) |
| 68 | ClassImp(AliAnalysisTaskHJetSpectra) |
| 69 | //________________________________________________________________________________________ |
| 70 | |
| 71 | AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra(): |
| 9d6abaad |
72 | AliAnalysisTaskSE(), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), |
| 2cc66a40 |
73 | fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fRhoTaskName(), |
| 74 | fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0), |
| 75 | fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5), fNumberOfCentralityBins(20), fCentralityType("V0A"), |
| 9d6abaad |
76 | fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0), |
| 2cc66a40 |
77 | fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0), |
| 78 | fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0), |
| 79 | fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0), |
| 80 | fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0), |
| 81 | fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0), |
| 82 | fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0), |
| 83 | fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0), |
| 9d6abaad |
84 | fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0), |
| 85 | fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0), |
| 2cc66a40 |
86 | fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0), |
| 9d6abaad |
87 | fhDphiTriggerJetMinBias(0x0),fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0), |
| 2cc66a40 |
88 | fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0), |
| 89 | fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize), |
| 9d6abaad |
90 | fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0), |
| 91 | fNofRandomCones(1), |
| 92 | fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR),fnRCones(16) |
| 2cc66a40 |
93 | { |
| 94 | //default constructor |
| 9d6abaad |
95 | for(Int_t k=0; k<50; k++){ |
| 96 | fRConePhi[k] = 0.0; |
| 97 | fRConeEta[k] = 0.0; |
| 98 | } |
| 2cc66a40 |
99 | } |
| 100 | |
| 101 | //________________________________________________________________________ |
| 102 | AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : |
| 9d6abaad |
103 | AliAnalysisTaskSE(name), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), |
| 2cc66a40 |
104 | fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fRhoTaskName(), |
| 105 | fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0), |
| 106 | fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5), fNumberOfCentralityBins(20), fCentralityType("V0A"), |
| 9d6abaad |
107 | fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0), |
| 2cc66a40 |
108 | fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0), |
| 109 | fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0), |
| 110 | fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0), |
| 111 | fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0), |
| 112 | fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0), |
| 113 | fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0), |
| 114 | fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0), |
| 9d6abaad |
115 | fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0), |
| 116 | fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0), |
| 2cc66a40 |
117 | fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0), |
| 9d6abaad |
118 | fhDphiTriggerJetMinBias(0x0), fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0), |
| 2cc66a40 |
119 | fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0), |
| 120 | fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize), |
| 9d6abaad |
121 | fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0), |
| 122 | fNofRandomCones(1), |
| 123 | fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR), fnRCones(16) |
| 2cc66a40 |
124 | { |
| 125 | //constructor that is called |
| 126 | //LIST OF TRACKS |
| 127 | fTrackArrayName = new TString(trackArrayName); |
| 9d6abaad |
128 | if((fTrackArrayName->Contains("MC") && fTrackArrayName->Contains("Particles")) || |
| 129 | (fTrackArrayName->Contains("mc") && fTrackArrayName->Contains("particles"))){ |
| 2cc66a40 |
130 | fIsKinematics = kTRUE; |
| 9d6abaad |
131 | } |
| 2cc66a40 |
132 | |
| 133 | //LIST of JETS |
| 134 | fJetArrayName = new TString(jetArrayName); |
| 135 | if(strcmp(fJetArrayName->Data(),"") == 0){ |
| 136 | AliError(Form("%s: Jet branch missing !", GetName())); |
| 137 | } |
| 138 | |
| 139 | //LIST OF JETS TO BE IGNORED WHILE RHO ESTIMATE |
| 140 | fBackgroundJetArrayName = new TString(backgroundJetArrayName); //jets to be removed from cell median rho estimate |
| 141 | if(strcmp(fBackgroundJetArrayName->Data(),"") == 0){ |
| 142 | AliError(Form("%s: Bg Jet branch missing !", GetName())); |
| 143 | } |
| 144 | |
| 9d6abaad |
145 | for(Int_t k=0; k<50; k++){ |
| 146 | fRConePhi[k] = 0.0; |
| 147 | fRConeEta[k] = 0.0; |
| 148 | } |
| 149 | |
| 2cc66a40 |
150 | DefineOutput(1, TList::Class()); |
| 151 | } |
| 9d6abaad |
152 | |
| 153 | //________________________________________________________________________ |
| 5e303be0 |
154 | void AliAnalysisTaskHJetSpectra::SetAnalyzeMC(Int_t val){ |
| 9d6abaad |
155 | if(val==1){ |
| 156 | fAnalyzePythia = kTRUE; |
| 157 | fAnalyzeHijing = kFALSE; |
| 158 | return; |
| 159 | } |
| 160 | if(val==2){ |
| 161 | fAnalyzeHijing = kTRUE; |
| 162 | fAnalyzePythia = kFALSE; |
| 163 | return; |
| 164 | } |
| 165 | |
| 166 | fAnalyzeHijing = kFALSE; |
| 167 | fAnalyzePythia = kFALSE; |
| 168 | return; |
| 169 | } |
| 2cc66a40 |
170 | //________________________________________________________________________ |
| 5e303be0 |
171 | Double_t AliAnalysisTaskHJetSpectra::GetConePt(Double_t eta, Double_t phi, Double_t radius){ |
| 2cc66a40 |
172 | //sum up pt inside a cone |
| 173 | Double_t tmpConePt = 0.0; |
| 174 | Double_t dphi = 0.0; |
| 175 | Double_t deta = 0.0; |
| 176 | Double_t radiussquared = radius*radius; |
| 177 | |
| 178 | for(Int_t i = 0; i < fTrackArray->GetEntries(); i++){ |
| 179 | AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i)); |
| 180 | if(!tmpTrack) continue; |
| 181 | if(IsTrackInAcceptance(tmpTrack)){ |
| 182 | dphi = RelativePhi(tmpTrack->Phi(),phi); |
| 183 | deta = tmpTrack->Eta() - eta; |
| 184 | if( dphi*dphi + deta*deta < radiussquared ){ |
| 185 | tmpConePt = tmpConePt + tmpTrack->Pt(); |
| 186 | } |
| 187 | } |
| 188 | } |
| 189 | return tmpConePt; |
| 190 | } |
| 191 | |
| 192 | |
| 193 | //________________________________________________________________________ |
| 5e303be0 |
194 | Double_t AliAnalysisTaskHJetSpectra::GetPtHard(){ |
| 9d6abaad |
195 | //get pt hard from pythia |
| b4e22a71 |
196 | AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader()); |
| 197 | if(MCEvent()){ |
| 198 | if(!pythiaHeader){ |
| 199 | // Check if AOD |
| 200 | AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); |
| 201 | |
| 202 | if(aodMCH){ |
| 203 | for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){ |
| 204 | pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i)); |
| 205 | if(pythiaHeader) break; |
| 206 | } |
| 207 | } |
| 208 | } |
| 209 | } |
| 210 | if(pythiaHeader){ |
| 9d6abaad |
211 | |
| c8ab9633 |
212 | TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); |
| 213 | |
| 214 | if(tree){ |
| 215 | TFile *curfile = tree->GetCurrentFile(); |
| 216 | if(!curfile) { |
| 217 | Error("Notify","No current file"); |
| 218 | return kFALSE; |
| 219 | } |
| 220 | Float_t xsection = 0.0; |
| 221 | Float_t trials = 1.0; |
| 222 | |
| 223 | AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials); |
| 224 | |
| 225 | fCrossSection = (Double_t) xsection;//pythiaHeader->GetXsection(); |
| 226 | |
| 227 | if(fCrossSection>0.){ //save cross-section and the number of trials |
| 228 | fTrials = (Double_t) trials; //pythiaHeader->Trials(); |
| 229 | fh1Xsec->Fill("<#sigma>", fCrossSection); |
| 230 | fh1Trials->Fill("#sum{ntrials}",fTrials); |
| 231 | } |
| 232 | } |
| b4e22a71 |
233 | return pythiaHeader->GetPtHard(); |
| 234 | } |
| 235 | AliWarning(Form("In task %s: GetPtHard() failed!", GetName())); |
| 236 | return -1.0; |
| 237 | } |
| b4e22a71 |
238 | //________________________________________________________________________ |
| 5e303be0 |
239 | Double_t AliAnalysisTaskHJetSpectra::GetImpactParameter(){ |
| 9d6abaad |
240 | //get impact parameter from hijing |
| 241 | AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()); |
| b4e22a71 |
242 | if(MCEvent()){ |
| 9d6abaad |
243 | if(!hijingHeader){ |
| b4e22a71 |
244 | // Check if AOD |
| 245 | AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); |
| 246 | |
| 247 | if(aodMCH){ |
| 9d6abaad |
248 | for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){ |
| 249 | hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i)); |
| 250 | if(hijingHeader) break; |
| b4e22a71 |
251 | } |
| 252 | } |
| 2cc66a40 |
253 | } |
| b4e22a71 |
254 | } |
| 9d6abaad |
255 | if(hijingHeader){ |
| 256 | return (Double_t) (hijingHeader->ImpactParameter()); |
| 257 | } |
| 258 | AliWarning(Form("In task %s: GetImpactParameter() failed!", GetName())); |
| 259 | return -1.0; |
| 260 | } |
| 261 | //________________________________________________________________________ |
| 5e303be0 |
262 | Double_t AliAnalysisTaskHJetSpectra::GetSimPrimaryVertex(){ |
| 9d6abaad |
263 | //get generator level primary vertex |
| 264 | AliGenEventHeader* mcHeader = NULL; |
| 265 | AliAODMCHeader* aodMCH = NULL; |
| 266 | if(MCEvent()){ |
| 267 | if(fAnalyzePythia){ |
| 268 | mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader()); |
| 269 | if(!mcHeader){ |
| 270 | // Check if AOD |
| 271 | aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); |
| 272 | |
| 273 | if(aodMCH){ |
| 274 | for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){ |
| 275 | mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i)); |
| 276 | if(mcHeader) break; |
| 277 | } |
| 278 | } |
| 279 | } |
| 280 | } |
| 281 | |
| 282 | if(fAnalyzeHijing){ |
| 283 | mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()); |
| 284 | if(!mcHeader){ |
| 285 | // Check if AOD |
| 286 | aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); |
| 287 | |
| 288 | if(aodMCH){ |
| 289 | for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){ |
| 290 | mcHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i)); |
| 291 | if(mcHeader) break; |
| 292 | } |
| 293 | } |
| 294 | } |
| 295 | } |
| 296 | } |
| 297 | if(mcHeader){ |
| 298 | |
| b4e22a71 |
299 | TArrayF pyVtx(3); |
| 9d6abaad |
300 | mcHeader->PrimaryVertex(pyVtx); |
| b4e22a71 |
301 | return (Double_t) (pyVtx[2]); |
| 302 | } |
| 303 | AliWarning(Form("In task %s: Pythia Vertex failed!", GetName())); |
| 304 | return 9999.0; |
| 305 | } |
| 2cc66a40 |
306 | |
| 2cc66a40 |
307 | |
| 2cc66a40 |
308 | |
| 309 | //________________________________________________________________________ |
| 5e303be0 |
310 | /*Double_t AliAnalysisTaskHJetSpectra::GetPythiaTrials() |
| 2cc66a40 |
311 | { |
| 312 | #ifdef DEBUGMODE |
| 313 | AliInfo("Starting GetPythiaTrials."); |
| 314 | #endif |
| 315 | AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader()); |
| 316 | if (MCEvent()) |
| 317 | if (!pythiaHeader) |
| 318 | { |
| 319 | // Check if AOD |
| 320 | AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); |
| 321 | |
| 322 | if (aodMCH) |
| 323 | { |
| 324 | for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) |
| 325 | { |
| 326 | pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i)); |
| 327 | if (pythiaHeader) break; |
| 328 | } |
| 329 | } |
| 330 | } |
| 331 | |
| 332 | #ifdef DEBUGMODE |
| 333 | AliInfo("Ending GetPythiaTrials."); |
| 334 | #endif |
| 335 | if (pythiaHeader) |
| 336 | return pythiaHeader->Trials(); |
| 337 | |
| 338 | AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName())); |
| 339 | return -1.0; |
| 340 | } |
| 341 | */ |
| 342 | |
| 343 | //________________________________________________________________________ |
| 344 | Double_t AliAnalysisTaskHJetSpectra::GetExternalRho(){ |
| 345 | // Get rho from event using CMS approach |
| 346 | |
| 347 | AliRhoParameter *rho = 0; |
| 348 | if(!fRhoTaskName.IsNull()) { |
| 349 | rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data())); |
| 350 | if(!rho){ |
| 351 | AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data())); |
| 352 | return 0.0; |
| 353 | } |
| 354 | }else return 0.0; |
| 355 | |
| 356 | return (rho->GetVal()); |
| 357 | } |
| 358 | |
| 359 | //________________________________________________________________________ |
| 5e303be0 |
360 | Bool_t AliAnalysisTaskHJetSpectra::IsEventInAcceptance(AliVEvent* event){ |
| 2cc66a40 |
361 | //EVENT SELECTION |
| 362 | |
| 363 | |
| 364 | if(!event) return kFALSE; |
| b4e22a71 |
365 | |
| 366 | //___________________________________________________ |
| 367 | |
| 9d6abaad |
368 | if(fAnalyzePythia || fAnalyzeHijing){ //PURE MC |
| b4e22a71 |
369 | if(!MCEvent()) return kFALSE; |
| 370 | |
| 371 | //BEFORE VERTEX CUT |
| 9d6abaad |
372 | Double_t vtxMC = GetSimPrimaryVertex(); |
| b4e22a71 |
373 | fhVertexZ->Fill(vtxMC); |
| 374 | |
| 375 | if(TMath::Abs(vtxMC) > 10.0){ |
| 376 | fHistEvtSelection->Fill(3); //count events rejected by vertex cut |
| 377 | return kFALSE; |
| 378 | } |
| 379 | fhVertexZAccept->Fill(vtxMC); |
| 380 | |
| 381 | return kTRUE; |
| 382 | } |
| 2cc66a40 |
383 | //___________________________________________________ |
| 384 | //TEST PILE UP |
| 385 | if(fUsePileUpCut){ |
| 386 | if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){ |
| 387 | fHistEvtSelection->Fill(2); //count events rejected by pileup |
| 388 | return kFALSE; |
| 389 | } |
| 390 | } |
| 391 | //___________________________________________________ |
| 392 | //BEFORE VERTEX CUT |
| 393 | fhVertexZ->Fill(event->GetPrimaryVertex()->GetZ()); |
| 394 | |
| 395 | if(fUseDefaultVertexCut){ |
| 396 | if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){ |
| 397 | fHistEvtSelection->Fill(3); //count events rejected by vertex cut |
| 398 | return kFALSE; |
| 399 | } |
| 400 | }else{ |
| 401 | if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > 10.0){ |
| 402 | fHistEvtSelection->Fill(3); //count events rejected by vertex cut |
| 403 | return kFALSE; |
| 404 | } |
| 405 | } |
| 406 | //___________________________________________________ |
| 407 | //AFTER VERTEX CUT |
| 408 | fhVertexZAccept->Fill(event->GetPrimaryVertex()->GetZ()); |
| 409 | |
| 410 | return kTRUE; |
| 411 | } |
| 412 | |
| 413 | //________________________________________________________________________ |
| 5e303be0 |
414 | Bool_t AliAnalysisTaskHJetSpectra::IsTrackInAcceptance(AliVParticle* track){ |
| 2cc66a40 |
415 | // Check if the track pt and eta range |
| 416 | if(track != 0){ |
| 417 | if(fIsKinematics){ |
| 418 | // TODO: Only working for AOD MC |
| 419 | if((!track->Charge()) || (!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) ) |
| 420 | return kFALSE; |
| 421 | } |
| 422 | if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT |
| 423 | if(track->Pt() >= fMinTrackPt){ //APPLY TRACK CUT |
| 424 | return kTRUE; |
| 425 | } |
| 426 | } |
| 427 | } |
| 428 | return kFALSE; |
| 429 | } |
| 430 | |
| 431 | //________________________________________________________________________ |
| 5e303be0 |
432 | Bool_t AliAnalysisTaskHJetSpectra::IsBackgroundJetInAcceptance(AliEmcalJet *jet){ |
| 2cc66a40 |
433 | //find jets to be removed from bg calculation |
| 434 | if(jet != 0){ |
| 435 | if(TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow){ |
| 436 | if(jet->Pt() >= fBackgroundJetPtMin){ //accept only hard jets |
| 437 | return kTRUE; |
| 438 | } |
| 439 | } |
| 440 | } |
| 441 | return kFALSE; |
| 442 | } |
| 443 | |
| 444 | //________________________________________________________________________ |
| 5e303be0 |
445 | Bool_t AliAnalysisTaskHJetSpectra::IsSignalJetInAcceptance(AliEmcalJet *jet){ |
| 2cc66a40 |
446 | //select jets in acceptance |
| 447 | if(jet == 0) return kFALSE; |
| 448 | if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){ |
| 449 | if(jet->Pt() >= fMinTrackPt){ |
| 450 | if(jet->Area() >= fMinJetArea){ |
| 451 | return kTRUE; |
| 452 | } |
| 453 | } |
| 454 | } |
| 455 | return kFALSE; |
| 456 | } |
| 457 | |
| 458 | |
| 459 | //________________________________________________________________________ |
| 460 | void AliAnalysisTaskHJetSpectra::ExecOnce(){ |
| 461 | //Read arrays of jets and tracks |
| 462 | |
| 463 | |
| 464 | fInitialized = kTRUE; //change flag to skip this function next time when processing UserExec |
| 465 | |
| 9d6abaad |
466 | |
| 467 | fnRCones = TMath::Nint(fRandConeRadiusSquared/fRConesRSquared); //the number of small R=0.1 random cones |
| 468 | |
| 2cc66a40 |
469 | // Check for track array |
| 470 | if(strcmp(fTrackArrayName->Data(), "") != 0){ |
| 471 | fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data())); |
| 472 | if(!fTrackArray){ |
| 473 | AliWarning(Form("%s: Could not retrieve tracks %s!", GetName(), fTrackArrayName->Data())); |
| 474 | }else{ |
| 475 | TClass *cl = fTrackArray->GetClass(); |
| 476 | if(!cl->GetBaseClass("AliVParticle")){ |
| 477 | AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data())); |
| 478 | fTrackArray = 0; |
| 479 | } |
| 480 | } |
| 481 | } |
| 482 | |
| 483 | // Check for jet array |
| 484 | if(strcmp(fJetArrayName->Data(), "") != 0){ |
| 485 | fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data())); |
| 486 | |
| 487 | if(!fJetArray){ |
| 488 | AliWarning(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrayName->Data())); |
| 489 | }else{ |
| 490 | if(!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")){ |
| 491 | AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data())); |
| 492 | fJetArray = 0; |
| 493 | } |
| 494 | } |
| 495 | } |
| 496 | |
| 497 | // Check for list of jets to be removed from background |
| 498 | if(strcmp(fBackgroundJetArrayName->Data(), "") != 0){ |
| 499 | fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data())); |
| 500 | if(!fBackgroundJetArray){ |
| 501 | AliInfo(Form("%s: Could not retrieve background jets %s!", GetName(), fBackgroundJetArrayName->Data())); |
| 502 | }else{ |
| 503 | if(!fBackgroundJetArray->GetClass()->GetBaseClass("AliEmcalJet")){ |
| 504 | AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fBackgroundJetArrayName->Data())); |
| 505 | fBackgroundJetArray = 0; |
| 506 | } |
| 507 | } |
| 508 | } |
| 509 | |
| 510 | // Look, if initialization is OK |
| 511 | |
| 512 | // Initialize helper class (for vertex selection & pile up correction) |
| 513 | fHelperClass = new AliAnalysisUtils(); |
| 514 | fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm |
| 515 | |
| 516 | return; |
| 517 | } |
| 518 | |
| 519 | |
| 520 | //________________________________________________________________________ |
| 521 | void AliAnalysisTaskHJetSpectra::GetDeltaPt(Double_t rho1, Double_t &dpt1, Double_t rho2, Double_t &dpt2, |
| 9d6abaad |
522 | Double_t rho3, Double_t &dpt3, |
| 523 | Double_t &rcPhi, Double_t &rcEta, |
| 524 | Double_t leadingJetExclusionProbability){ |
| 2cc66a40 |
525 | |
| 526 | //delta pt = random cone - rho |
| 527 | |
| 528 | // Define an invalid delta pt |
| 529 | dpt1 = -10000.0; |
| 530 | dpt2 = -10000.0; |
| 531 | dpt3 = -10000.0; |
| 532 | |
| 533 | // Define eta range |
| 534 | Double_t etaMin, etaMax; |
| 535 | etaMin = -(fTrackEtaWindow-fRandConeRadius); |
| 536 | etaMax = +(fTrackEtaWindow-fRandConeRadius); |
| 537 | |
| 538 | // Define random cone Eta+Phi |
| 539 | Bool_t coneValid = kTRUE; |
| 540 | Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin); |
| 541 | Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi(); |
| 542 | |
| 543 | // if there is a jet, check for overlap if demanded |
| 544 | if(leadingJetExclusionProbability){ |
| 545 | AliEmcalJet* tmpLeading = NULL; |
| 546 | Double_t lpt = -1.0; |
| 547 | // Get leading jet (regardless of pT) |
| 548 | for(Int_t i = 0; i<fJetArray->GetEntries(); i++){ |
| 549 | AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i)); |
| 550 | if(!tmpJet) continue; |
| 551 | if((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea)){ |
| 552 | if(tmpJet->Pt() > lpt){ |
| 553 | tmpLeading = tmpJet; |
| 554 | lpt = tmpJet->Pt(); |
| 555 | } |
| 556 | } |
| 557 | } |
| 558 | if(tmpLeading){ |
| 559 | Double_t excludedJetPhi = tmpLeading->Phi(); |
| 560 | Double_t tmpDeltaPhi = RelativePhi(tmpRandConePhi, excludedJetPhi); |
| 561 | Double_t excludedJetEta = tmpLeading->Eta()-tmpRandConeEta; |
| 562 | |
| 563 | // Check, if cone has overlap with jet |
| 564 | if(tmpDeltaPhi*tmpDeltaPhi + excludedJetEta*excludedJetEta <= fRandConeRadiusSquared){ |
| 565 | // Define probability to exclude the RC |
| 566 | Double_t probability = leadingJetExclusionProbability; |
| 567 | |
| 568 | // Only exclude cone with a given probability |
| 569 | if(fRandom->Rndm()<=probability) coneValid = kFALSE; |
| 570 | } |
| 571 | } |
| 572 | } |
| 573 | |
| 9d6abaad |
574 | rcPhi = 9999.0; |
| 575 | rcEta = 9999.0; |
| 2cc66a40 |
576 | |
| 577 | // Get the cones' pt and calculate delta pt |
| 578 | if(coneValid){ |
| 9d6abaad |
579 | rcPhi = tmpRandConePhi; |
| 580 | rcEta = tmpRandConeEta; |
| 2cc66a40 |
581 | Double_t conePt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius); |
| 582 | dpt1 = conePt - (rho1*fRandConeRadiusSquared*TMath::Pi()); |
| 583 | dpt2 = conePt - (rho2*fRandConeRadiusSquared*TMath::Pi()); |
| 584 | dpt3 = conePt - (rho3*fRandConeRadiusSquared*TMath::Pi()); |
| 585 | } |
| 586 | |
| 587 | } |
| 588 | |
| 589 | |
| 590 | //________________________________________________________________________ |
| 591 | void AliAnalysisTaskHJetSpectra::Calculate(AliVEvent* event){ |
| 592 | //Analyze the event and Fill histograms |
| 593 | |
| b4e22a71 |
594 | if(fAnalyzePythia){ |
| 595 | fh1PtHard->Fill(GetPtHard()); |
| 596 | } |
| 597 | |
| 9d6abaad |
598 | if(fAnalyzeHijing){ |
| 599 | fImpParam = GetImpactParameter(); |
| 600 | fhImpactParameter->Fill(fImpParam); |
| 601 | } |
| 2cc66a40 |
602 | //_________________________________________________________________ |
| 603 | // FILL EVENT STATISTICS |
| 604 | fHistEvtSelection->Fill(1); //Count input event |
| 605 | |
| 606 | if(!IsEventInAcceptance(event)) return; //post data is in UserExec |
| 607 | |
| 608 | |
| 609 | // Get centrality |
| 610 | AliCentrality* tmpCentrality = event->GetCentrality(); |
| 611 | if(!tmpCentrality){ |
| 612 | fHistEvtSelection->Fill(4); |
| 613 | return; //post data is in UserExec |
| 614 | } |
| 615 | Double_t centralityPercentile = -1.0; |
| 616 | Double_t centralityPercentileV0A = -1.0; |
| 617 | Double_t centralityPercentileV0C = -1.0; |
| 618 | Double_t centralityPercentileV0M = -1.0; |
| 619 | Double_t centralityPercentileZNA = -1.0; |
| 620 | if(tmpCentrality != NULL){ |
| 621 | centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data()); |
| 622 | centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A"); |
| 623 | centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C"); |
| 624 | centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M"); |
| 625 | centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA"); |
| 626 | } |
| 627 | |
| 628 | if((centralityPercentile < 0.0) || (centralityPercentile > 100.0)){ |
| 629 | AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile)); |
| 630 | fHistEvtSelection->Fill(4); |
| 631 | return; |
| 632 | } |
| 633 | fhCentrality->Fill(centralityPercentile); |
| 634 | fhCentralityV0M->Fill(centralityPercentileV0M); |
| 635 | fhCentralityV0A->Fill(centralityPercentileV0A); |
| 636 | fhCentralityV0C->Fill(centralityPercentileV0C); |
| 637 | fhCentralityZNA->Fill(centralityPercentileZNA); |
| 638 | |
| 639 | fHistEvtSelection->Fill(0); //Count input event |
| 640 | |
| 641 | // END EVENT SELECTION |
| 642 | //___________________________________________________________ |
| 643 | |
| 644 | //LOOP OVER TRACKS SEARCH FOR TRIGGER |
| 645 | std::vector<Int_t> trigTracks; //list pf trigger particle indices |
| 9d6abaad |
646 | //Bool_t bContainesHighPtTrack = kFALSE; |
| 2cc66a40 |
647 | |
| 648 | Int_t nTracks = fTrackArray->GetEntries(); |
| 649 | |
| 650 | for(Int_t i = 0; i < nTracks; i++){ |
| 651 | AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i)); |
| 652 | |
| 653 | if(!track) continue; |
| 654 | |
| 655 | if(IsTrackInAcceptance(track)){ |
| 656 | fhTrackPhi->Fill(track->Pt(), RelativePhi(track->Phi(),0.0)); // phi = -pi az pi |
| 657 | fhTrackEta->Fill(track->Pt(), track->Eta()); |
| 658 | fhTrackCentVsPt->Fill(track->Pt(), centralityPercentile); |
| 659 | |
| 660 | if(fTTlow <= track->Pt() && track->Pt() < fTThigh){ |
| 661 | trigTracks.push_back(i); //trigger candidates |
| 662 | } |
| 663 | |
| 9d6abaad |
664 | //if(track->Pt()>=8.0) bContainesHighPtTrack = kTRUE; |
| 2cc66a40 |
665 | } |
| 666 | } |
| 667 | |
| 668 | Int_t ntriggers = (Int_t) trigTracks.size(); |
| 669 | Int_t indexSingleRndTrig = -1; //index of single random trigger |
| 670 | Double_t areaJet, pTJet; |
| 671 | Double_t tmpArray[3]; |
| 672 | Double_t rhoFromCellMedian = 0.0; //UE density cell median |
| 673 | Double_t rhoCone = 0.0; //UE density perp cone |
| 674 | Double_t rhoCMS = 0.0; //UE density ala CMS |
| 9d6abaad |
675 | Double_t deltaptCellMedian, deltaptCone, deltaptCMS, randConePhi, randConeEta; |
| 676 | Double_t distanceFromTrigger; |
| 677 | |
| 678 | if(ntriggers>0){ |
| 679 | if(fTTtype==0){ //select single inclusive trigger |
| 680 | indexSingleRndTrig = fRandom->Integer(ntriggers); //Integer 0 ... ntriggers-1 |
| 681 | } |
| 682 | } |
| 2cc66a40 |
683 | |
| 684 | rhoFromCellMedian = EstimateBgRhoMedian(); |
| 685 | rhoCone = EstimateBgCone(); |
| 686 | rhoCMS = GetExternalRho(); |
| 687 | |
| b4e22a71 |
688 | fhRhoCellMedianIncl->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile); |
| 689 | fhRhoConeIncl->Fill( (Float_t) rhoCone, (Float_t) centralityPercentile); |
| 690 | fhRhoCMSIncl->Fill( (Float_t) rhoCMS, (Float_t) centralityPercentile); |
| 9d6abaad |
691 | |
| 2cc66a40 |
692 | for(Int_t irc=0; irc<fNofRandomCones; irc++){ //generate 4 random cones per event |
| 9d6abaad |
693 | GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0); |
| 2cc66a40 |
694 | |
| 695 | fhDeltaPtMedianIncl->Fill(deltaptCellMedian, (Double_t) centralityPercentile); |
| 696 | fhDeltaPtConeIncl->Fill( deltaptCone, (Double_t) centralityPercentile); |
| 697 | fhDeltaPtCMSIncl->Fill( deltaptCMS, (Double_t) centralityPercentile); |
| b4e22a71 |
698 | |
| 9d6abaad |
699 | if(ntriggers>0){ |
| 700 | //fill delta pt histograms near side + away side |
| 2cc66a40 |
701 | fhDeltaPtMedian->Fill( deltaptCellMedian, (Double_t) centralityPercentile); |
| 9d6abaad |
702 | fhDeltaPtCone->Fill( deltaptCone, (Double_t) centralityPercentile); |
| 703 | fhDeltaPtCMS->Fill( deltaptCMS, (Double_t) centralityPercentile); |
| 704 | |
| 705 | if(indexSingleRndTrig>-1){ |
| 706 | AliVTrack* triggHad = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[indexSingleRndTrig])); |
| 707 | Double_t dphiTrigRC = RelativePhi(triggHad->Phi(), randConePhi); |
| 708 | Double_t detaTrigRC = triggHad->Eta()- randConeEta; |
| 709 | if(TMath::Abs(dphiTrigRC)< TMath::Pi()/2){ //near side |
| 710 | fhDeltaPtMedianNearSide->Fill( deltaptCellMedian, (Double_t) centralityPercentile); |
| 711 | fhDeltaPtCMSNearSide->Fill( deltaptCMS, (Double_t) centralityPercentile); |
| 712 | }else{ //away side |
| 713 | fhDeltaPtMedianAwaySide->Fill( deltaptCellMedian, (Double_t) centralityPercentile); |
| 714 | fhDeltaPtCMSAwaySide->Fill( deltaptCMS, (Double_t) centralityPercentile); |
| 715 | } |
| 716 | |
| 717 | distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC); |
| 718 | while(distanceFromTrigger<0.5 + fRandConeRadius){ |
| 719 | GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0); |
| 720 | dphiTrigRC = RelativePhi(triggHad->Phi(), randConePhi); |
| 721 | detaTrigRC = triggHad->Eta()- randConeEta; |
| 722 | distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC); |
| 723 | } |
| 724 | if(distanceFromTrigger>0.5 + fRandConeRadius){ |
| 725 | fhDeltaPtMedianExclTrigCone->Fill( deltaptCellMedian, (Double_t) centralityPercentile); |
| 726 | fhDeltaPtCMSExclTrigCone->Fill( deltaptCMS, (Double_t) centralityPercentile); |
| 727 | } |
| 728 | } |
| 2cc66a40 |
729 | } |
| 730 | } |
| 9d6abaad |
731 | |
| 732 | //_______________________________________ |
| 733 | Int_t idxLeadingJetAwaySide = -1; |
| 734 | Double_t ptLeadingJetAwaySide = -1.0; |
| 735 | Double_t phiTrigger = -1.0; //-pi,pi |
| 2cc66a40 |
736 | |
| 737 | if(ntriggers>0){ |
| 738 | //Estimate UE density |
| 739 | //Fill once per event |
| 740 | fhRhoCellMedian->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile); |
| 741 | fhRhoCone->Fill( (Float_t) rhoCone, (Float_t) centralityPercentile); |
| 742 | fhRhoCMS->Fill( (Float_t) rhoCMS, (Float_t) centralityPercentile); |
| 743 | |
| 744 | |
| 2cc66a40 |
745 | //TRIGGER PARTICLE LOOP |
| 746 | for(Int_t it=0; it<ntriggers; it++){ //loop over trigger configurations |
| 9d6abaad |
747 | |
| 2cc66a40 |
748 | if(fTTtype==0){ |
| 749 | if(it != indexSingleRndTrig) continue; |
| 750 | } |
| 751 | |
| 752 | AliVTrack* triggerHadron = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[it])); |
| 753 | if(!triggerHadron) continue; |
| 754 | |
| 755 | fh2Ntriggers->Fill((Float_t) centralityPercentile, (Float_t) triggerHadron->Pt()); //trigger p |
| 756 | |
| 9d6abaad |
757 | if(fAnalyzeHijing){ //impact parameter for triggered events |
| 758 | fhImpactParameterTT->Fill(fImpParam); |
| 759 | } |
| 760 | |
| 761 | phiTrigger = RelativePhi(triggerHadron->Phi(),0.0); //-pi,pi |
| 762 | |
| 2cc66a40 |
763 | //JET LOOP |
| 764 | for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){ |
| 765 | AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij)); |
| 766 | if(!jet){ |
| 767 | AliError(Form("%s: Could not receive jet %d", GetName(), ij)); |
| 768 | continue; |
| 769 | } |
| 770 | if(!IsSignalJetInAcceptance(jet)) continue; |
| 771 | |
| 772 | areaJet = jet->Area(); |
| 773 | pTJet = jet->Pt(); |
| 774 | |
| 775 | if(it==0 || it == indexSingleRndTrig){ |
| 776 | fhJetPhi->Fill( pTJet, RelativePhi(jet->Phi(),0.0)); |
| 777 | fhJetEta->Fill( pTJet, jet->Eta()); |
| 778 | } |
| 779 | |
| 780 | Double_t dphi = RelativePhi(triggerHadron->Phi(), jet->Phi()); |
| 781 | |
| 782 | Double_t dfi = dphi; //-0.5*pi to 1.5*Pi |
| 783 | if(dfi<-0.5*TMath::Pi()) dfi += 2*TMath::Pi(); |
| 784 | if(dfi> 1.5*TMath::Pi()) dfi -= 2*TMath::Pi(); |
| 9d6abaad |
785 | fhDphiTriggerJetMinBias->Fill((Float_t) jet->Pt(),(Float_t) dfi); |
| 786 | if(centralityPercentile<20.) fhDphiTriggerJetCent20->Fill((Float_t) jet->Pt(),(Float_t) dfi); |
| 2cc66a40 |
787 | //------------------------- |
| 788 | |
| 789 | if(TMath::Abs(dphi) < fDphiCut) continue; //Dphi cut between trigger and assoc |
| 790 | fhDphiTriggerJetAccept->Fill(dfi); //Accepted |
| 791 | |
| 9d6abaad |
792 | if(pTJet > ptLeadingJetAwaySide){ //search for the leading away side jet |
| 793 | idxLeadingJetAwaySide = ij; |
| 794 | ptLeadingJetAwaySide = pTJet; |
| 795 | } |
| 796 | |
| 2cc66a40 |
797 | //Centrality, A, pTjet |
| 798 | tmpArray[0] = centralityPercentile; |
| 799 | tmpArray[1] = areaJet; |
| 800 | tmpArray[2] = pTJet; |
| 801 | fHJetSpec->Fill(tmpArray); |
| 802 | |
| 803 | //Subtract cell median |
| 804 | tmpArray[2] = pTJet - areaJet*rhoFromCellMedian; |
| 805 | fHJetSpecSubUeMedian->Fill(tmpArray); |
| 806 | fARhoCellMedian->Fill((Float_t) (areaJet*rhoFromCellMedian)); |
| 807 | |
| 808 | //Subtract perp cone |
| 809 | tmpArray[2] = pTJet - areaJet*rhoCone; |
| 810 | fHJetSpecSubUeCone->Fill(tmpArray); |
| 811 | fARhoCone->Fill((Float_t) (areaJet*rhoCone)); |
| 812 | |
| 813 | //Subtract CMS bg |
| 814 | tmpArray[2] = pTJet - areaJet*rhoCMS; |
| 815 | fHJetSpecSubUeCMS->Fill(tmpArray); |
| 816 | fARhoCMS->Fill((Float_t) (areaJet*rhoCMS)); |
| 817 | |
| 818 | }//JET LOOP |
| 819 | } |
| 820 | } |
| 821 | |
| 9d6abaad |
822 | //_______________________________________ |
| 823 | // Get delta phi from small R=0.1 cones |
| 824 | if(ntriggers>0 && indexSingleRndTrig>-1){ |
| 825 | |
| 826 | AliEmcalJet* jet = NULL; |
| 827 | Double_t phiExclJet =0., etaExclJet = 9999., rExclJet = 0.0; |
| 828 | |
| 829 | if(idxLeadingJetAwaySide>-1){ |
| 830 | jet = static_cast<AliEmcalJet*>(fJetArray->At(idxLeadingJetAwaySide)); |
| 831 | if(!jet){ |
| 832 | AliError(Form("%s: Could not receive leading jet %d", GetName(), idxLeadingJetAwaySide)); |
| 833 | }else{ |
| 834 | phiExclJet = jet->Phi(); |
| 835 | etaExclJet = jet->Eta(); |
| 836 | rExclJet = TMath::Sqrt(jet->Area()/TMath::Pi()); |
| 837 | } |
| 838 | } |
| 839 | |
| 840 | Int_t countPlacedRcones = 0; |
| 841 | Int_t inwhile=0; |
| 842 | Double_t dphiMaxFromPiForRcones = TMath::Pi() - fDphiCut + fRandConeRadius - fRConesR; |
| 843 | Double_t detaMaxForRcones = fTrackEtaWindow - fRConesR; |
| 844 | Double_t rcphi,rceta; |
| 845 | Bool_t goodrc; |
| 846 | |
| 847 | while(countPlacedRcones < fnRCones){ //generate fnRCones cones of radius R=0.1 |
| 848 | inwhile++; |
| 849 | if(inwhile>500){ |
| 850 | AliError(Form("%s: Small space where to put another random cone %d", GetName(), idxLeadingJetAwaySide)); |
| 851 | break; |
| 852 | } |
| 853 | rcphi = RelativePhi(phiTrigger + TMath::Pi() + dphiMaxFromPiForRcones*fRandom->Uniform(-1.0,1.0), 0.0); //-pi,pi |
| 854 | rceta = detaMaxForRcones*fRandom->Uniform(-1.0,1.0); |
| 855 | if(jet){//do not merge random cones with the leading jet |
| 856 | if(!DistantCones(rcphi,rceta,fRConesR, phiExclJet, etaExclJet, rExclJet)) continue; |
| 857 | } |
| 858 | |
| 859 | goodrc = kTRUE; //generate disjoint random cones |
| 860 | for(int k=0; k<countPlacedRcones; k++){ |
| 861 | if(!DistantCones(rcphi, rceta, fRConesR, (Double_t) fRConePhi[k], (Double_t) fRConeEta[k], fRConesR)){ |
| 862 | goodrc = kFALSE; |
| 863 | break; |
| 864 | } |
| 865 | }//end loop over already placed cones |
| 866 | |
| 867 | if(!goodrc) continue; |
| 868 | fRConePhi[countPlacedRcones] = rcphi; |
| 869 | fRConeEta[countPlacedRcones] = rceta; |
| 870 | countPlacedRcones++; |
| 871 | }//end of loop generating small R=0.1 random cones |
| 872 | //sum track pT in the random cones |
| 873 | Double_t sumPtofRandomCones = 0.0; |
| 874 | |
| 875 | for(Int_t itrk = 0; itrk < nTracks; itrk++){ |
| 876 | AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(itrk)); |
| 877 | if(!track) continue; |
| 878 | |
| 879 | if(IsTrackInAcceptance(track)){ |
| 880 | for(int k=0; k<countPlacedRcones; k++){ |
| 881 | |
| 882 | rcphi = RelativePhi(track->Phi(), fRConePhi[k]); // phi = -pi az pi |
| 883 | rceta = track->Eta() - fRConeEta[k]; |
| 884 | if(rcphi*rcphi + rceta*rceta < fRConesRSquared){ |
| 885 | sumPtofRandomCones += track->Pt(); //track is in the cone |
| 886 | break; |
| 887 | } |
| 888 | }//loop over cones |
| 889 | } |
| 890 | }//loop over tracks |
| 891 | Double_t totarea = countPlacedRcones*TMath::Pi()*fRConesRSquared; |
| 892 | fhDeltaPtMedianExclAwayJet->Fill( sumPtofRandomCones - rhoFromCellMedian*totarea, (Double_t) centralityPercentile ); |
| 893 | fhDeltaPtCMSExclAwayJet->Fill( sumPtofRandomCones - rhoCMS*totarea , (Double_t) centralityPercentile ); |
| 894 | |
| 895 | }//end delta phi from small R=0.1 random cones |
| 896 | |
| 2cc66a40 |
897 | return; |
| 898 | } |
| 899 | |
| 900 | //________________________________________________________________________ |
| 901 | Bool_t AliAnalysisTaskHJetSpectra::UserNotify(){ |
| 902 | // Implemented Notify() to read the cross sections |
| 903 | // and number of trials from pyxsec.root |
| 9d6abaad |
904 | /* |
| 905 | if(fAnalyzePythia){ |
| 906 | |
| 2cc66a40 |
907 | TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); |
| 908 | TFile *currFile = tree->GetCurrentFile(); |
| 909 | |
| 910 | TString file(currFile->GetName()); |
| 911 | |
| 912 | if(file.Contains("root_archive.zip#")){ |
| 913 | Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact); |
| 914 | Ssiz_t pos = file.Index("#",1,pos1,TString::kExact); |
| 915 | file.Replace(pos+1,20,""); |
| 916 | } |
| 917 | else { |
| 918 | // not an archive take the basename.... |
| 919 | file.ReplaceAll(gSystem->BaseName(file.Data()),""); |
| 920 | } |
| 921 | |
| 922 | TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root |
| 923 | if(!fxsec){ |
| 924 | // next trial fetch the histgram file |
| 925 | fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")); |
| 926 | if(!fxsec){ |
| 927 | // not a severe condition but inciate that we have no information |
| 928 | return kFALSE; |
| 929 | } |
| 930 | else{ |
| 931 | // find the tlist we want to be independtent of the name so use the Tkey |
| 932 | TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); |
| 933 | if(!key){ |
| 934 | fxsec->Close(); |
| 935 | return kFALSE; |
| 936 | } |
| 937 | TList *list = dynamic_cast<TList*>(key->ReadObj()); |
| 938 | if(!list){ |
| 939 | fxsec->Close(); |
| 940 | return kFALSE; |
| 941 | } |
| 942 | fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1); |
| 943 | fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1); |
| 944 | fxsec->Close(); |
| 945 | } |
| 946 | } // no tree pyxsec.root |
| 947 | else { |
| 948 | TTree *xtree = (TTree*)fxsec->Get("Xsection"); |
| 949 | if(!xtree){ |
| 950 | fxsec->Close(); |
| 951 | return kFALSE; |
| 952 | } |
| 953 | UInt_t ntrials = 0; |
| 954 | Double_t xsection = 0; |
| 955 | xtree->SetBranchAddress("xsection",&xsection); |
| 956 | xtree->SetBranchAddress("ntrials",&ntrials); |
| 957 | xtree->GetEntry(0); |
| 958 | fTrials = ntrials; |
| 959 | fCrossSection = xsection; |
| 960 | fxsec->Close(); |
| 961 | } |
| 962 | |
| 9d6abaad |
963 | |
| 2cc66a40 |
964 | fh1Xsec->Fill("<#sigma>", fCrossSection); |
| 965 | fh1Trials->Fill("#sum{ntrials}",fTrials); |
| 9d6abaad |
966 | |
| 2cc66a40 |
967 | } |
| 9d6abaad |
968 | */ |
| 2cc66a40 |
969 | return kTRUE; |
| 970 | } |
| 971 | |
| 972 | //________________________________________________________________________ |
| 973 | void AliAnalysisTaskHJetSpectra::Terminate(Option_t *){ |
| 974 | //Treminate |
| 975 | PostData(1, fOutputList); |
| 976 | |
| 977 | // Mandatory |
| 978 | fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot |
| 979 | if(!fOutputList) { |
| 980 | printf("ERROR: Output list not available\n"); |
| 981 | return; |
| 982 | } |
| 983 | } |
| 984 | |
| 985 | //________________________________________________________________________ |
| 986 | AliAnalysisTaskHJetSpectra::~AliAnalysisTaskHJetSpectra(){ |
| 987 | // Destructor. Clean-up the output list, but not the histograms that are put inside |
| 988 | // (the list is owner and will clean-up these histograms). Protect in PROOF case. |
| 989 | if(fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { |
| 990 | delete fOutputList; |
| 991 | } |
| 992 | delete fRandom; |
| b4e22a71 |
993 | delete fTrackArrayName; |
| 994 | delete fJetArrayName; |
| 995 | delete fBackgroundJetArrayName; |
| 996 | delete fHelperClass; |
| 997 | |
| 2cc66a40 |
998 | } |
| 999 | |
| 1000 | //________________________________________________________________________ |
| 1001 | void AliAnalysisTaskHJetSpectra::UserCreateOutputObjects(){ |
| 1002 | // called once to create user defined output objects like histograms, plots etc. |
| 1003 | // and to put it on the output list. |
| 1004 | // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects. |
| 1005 | |
| 1006 | fRandom = new TRandom3(0); |
| 1007 | |
| 1008 | fOutputList = new TList(); |
| 1009 | fOutputList->SetOwner(); // otherwise it produces leaks in merging |
| 1010 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
| 1011 | TH1::AddDirectory(kFALSE); |
| 1012 | |
| 1013 | //__________________________________________________________ |
| 1014 | // Event statistics |
| 1015 | fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5); |
| 1016 | fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); |
| 1017 | fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN"); |
| 1018 | fHistEvtSelection->GetXaxis()->SetBinLabel(3,"pile up (rejected)"); |
| 1019 | fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)"); |
| 1020 | fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)"); |
| 1021 | |
| 1022 | fOutputList->Add(fHistEvtSelection); |
| 1023 | //___________________________________________________________ |
| 1024 | // Hard trigger counter |
| 1025 | fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers", |
| 1026 | fNumberOfCentralityBins,0.0,100.0,50,0.0,50.0); |
| 1027 | fOutputList->Add(fh2Ntriggers); |
| 1028 | //___________________________________________________________ |
| 1029 | // trigger associated jet spectra (jet pT not corrected for UE) |
| 1030 | Int_t bw = (fUseDoubleBinPrecision==0) ? 1 : 2; //make larger bin width |
| 1031 | |
| 1032 | //jet associated to given TT |
| 1033 | //Centrality, A, pTjet |
| 1034 | const Int_t dimSpec = 3; |
| 1035 | const Int_t nBinsSpec[dimSpec] = {fNumberOfCentralityBins, 50, bw*110}; |
| 1036 | const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20.0}; |
| 1037 | const Double_t hiBinSpec[dimSpec] = {100.0, 1.5, 200.0}; |
| 1038 | fHJetSpec = new THnSparseF("fHJetSpec", |
| 1039 | "Recoil jet spectrum [cent,A,pTjet]", |
| 1040 | dimSpec,nBinsSpec,lowBinSpec,hiBinSpec); |
| 1041 | fOutputList->Add(fHJetSpec); |
| 1042 | //___________________________________________________________ |
| 1043 | //jet associated to given TT (jet pT corrected with rho from cell median) |
| 1044 | fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian"); |
| 1045 | fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]"); |
| 1046 | fOutputList->Add(fHJetSpecSubUeMedian); |
| 1047 | //___________________________________________________________ |
| 1048 | //jet associated to given TT (jet pT corrected with rho from perp cone) |
| 1049 | fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone"); |
| 1050 | fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]"); |
| 1051 | fOutputList->Add(fHJetSpecSubUeCone); |
| 1052 | //___________________________________________________________ |
| 1053 | //jet associated to given TT (jet pT corrected with rho from CMS approach) |
| 1054 | fHJetSpecSubUeCMS = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCMS"); |
| 1055 | fHJetSpecSubUeCMS->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]"); |
| 1056 | fOutputList->Add(fHJetSpecSubUeCMS); |
| 1057 | |
| 1058 | //____________________________________________________________________ |
| 1059 | //UE from cell median [Centrality, rho, pTUe ] |
| 1060 | |
| 1061 | fhRhoCellMedian = new TH2F("fhRhoCellMedian","Rho",40, 0.0, 20.0, fNumberOfCentralityBins, 0.0, 100.); |
| 1062 | fOutputList->Add(fhRhoCellMedian); |
| 1063 | |
| 1064 | fhRhoCone = (TH2F*) fhRhoCellMedian->Clone("fhRhoCone"); |
| 1065 | fOutputList->Add(fhRhoCone); |
| 1066 | |
| 1067 | fhRhoCMS = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMS"); |
| 1068 | fOutputList->Add(fhRhoCMS); |
| 1069 | |
| 1070 | fhRhoCellMedianIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCellMedianIncl"); |
| 1071 | fOutputList->Add(fhRhoCellMedianIncl); |
| 1072 | |
| 1073 | fhRhoConeIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoConeIncl"); |
| 1074 | fOutputList->Add(fhRhoConeIncl); |
| 1075 | |
| 1076 | fhRhoCMSIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMSIncl"); |
| 1077 | fOutputList->Add(fhRhoCMSIncl); |
| 1078 | |
| 1079 | //_______________________________________________________________________ |
| 1080 | // rho times area |
| 1081 | fARhoCellMedian = new TH1F("fARhoCellMedian","Area times rho",40, 0.0, 20.0); |
| 1082 | fOutputList->Add(fARhoCellMedian); |
| 1083 | |
| 1084 | fARhoCone = (TH1F*) fARhoCellMedian->Clone("fARhoCone"); |
| 1085 | fOutputList->Add(fARhoCone); |
| 1086 | |
| 1087 | fARhoCMS = (TH1F*) fARhoCellMedian->Clone("fARhoCMS"); |
| 1088 | fOutputList->Add(fARhoCMS); |
| 1089 | |
| 1090 | //_______________________________________________________________________ |
| 1091 | // Delta pt distributions |
| 1092 | fhDeltaPtMedian = new TH2D("fhDeltaPtMedian","DeltaPt", bw*110, -20, 200, fNumberOfCentralityBins,0.0,100.0); |
| 1093 | fOutputList->Add(fhDeltaPtMedian); |
| 1094 | |
| 1095 | fhDeltaPtCone = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCone"); |
| 1096 | fOutputList->Add(fhDeltaPtCone); |
| 1097 | |
| 1098 | fhDeltaPtCMS = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMS"); |
| 1099 | fOutputList->Add(fhDeltaPtCMS); |
| 1100 | |
| 1101 | fhDeltaPtMedianIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianIncl"); |
| 1102 | fOutputList->Add(fhDeltaPtMedianIncl); |
| 1103 | |
| 1104 | fhDeltaPtConeIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtConeIncl"); |
| 1105 | fOutputList->Add(fhDeltaPtConeIncl); |
| 1106 | |
| 1107 | fhDeltaPtCMSIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSIncl"); |
| 1108 | fOutputList->Add(fhDeltaPtCMSIncl); |
| 9d6abaad |
1109 | |
| 1110 | fhDeltaPtMedianNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianNearSide"); |
| 1111 | fOutputList->Add(fhDeltaPtMedianNearSide); |
| 1112 | |
| 1113 | fhDeltaPtMedianAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianAwaySide"); |
| 1114 | fOutputList->Add(fhDeltaPtMedianAwaySide); |
| 1115 | |
| 1116 | fhDeltaPtCMSNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSNearSide"); |
| 1117 | fOutputList->Add(fhDeltaPtCMSNearSide); |
| 1118 | |
| 1119 | fhDeltaPtCMSAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSAwaySide"); |
| 1120 | fOutputList->Add(fhDeltaPtCMSAwaySide); |
| 1121 | |
| 1122 | fhDeltaPtMedianExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclTrigCone"); |
| 1123 | fOutputList->Add(fhDeltaPtMedianExclTrigCone); |
| 1124 | |
| 1125 | fhDeltaPtCMSExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclTrigCone"); |
| 1126 | fOutputList->Add(fhDeltaPtCMSExclTrigCone); |
| 1127 | |
| 1128 | fhDeltaPtMedianExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclAwayJet"); |
| 1129 | fOutputList->Add(fhDeltaPtMedianExclAwayJet); |
| 1130 | |
| 1131 | fhDeltaPtCMSExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclAwayJet"); |
| 1132 | fOutputList->Add(fhDeltaPtCMSExclAwayJet); |
| 1133 | |
| 2cc66a40 |
1134 | //_______________________________________________________________________ |
| 1135 | //inclusive azimuthal and pseudorapidity histograms |
| 1136 | fhJetPhi = new TH2F("fhJetPhi","Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi()); |
| 1137 | fOutputList->Add(fhJetPhi); |
| 1138 | //------------------------- |
| 1139 | fhTrackPhi = new TH2F("fhTrackPhi","azim dist trig had vs pT,trk", 50, 0, 50, 50,-TMath::Pi(),TMath::Pi()); |
| 1140 | fOutputList->Add(fhTrackPhi); |
| 1141 | //------------------------- |
| 1142 | fhJetEta = new TH2F("fhJetEta","Eta dist jets vs pTjet", 50,0, 100, 40,-0.9,0.9); |
| 1143 | fOutputList->Add(fhJetEta); |
| 1144 | //------------------------- |
| 1145 | fhTrackEta = new TH2F("fhTrackEta","Eta dist trig had vs pT,trk", 50, 0, 50, 40,-0.9,0.9); |
| 1146 | fOutputList->Add(fhTrackEta); |
| 1147 | //------------------------- |
| 1148 | fhTrackCentVsPt = new TH2F("fhTrackCentVsPt","pT,trk vs centrality", 50, 0, 50, fNumberOfCentralityBins,0,100); |
| 1149 | fOutputList->Add(fhTrackCentVsPt); |
| 1150 | //------------------------- |
| 1151 | fhVertexZ = new TH1F("fhVertexZ","z vertex",40,-20,20); |
| 1152 | fOutputList->Add(fhVertexZ); |
| 1153 | //------------------------- |
| 1154 | fhVertexZAccept = new TH1F("fhVertexZAccept","z vertex after cut",40,-20,20); |
| 1155 | fOutputList->Add(fhVertexZAccept); |
| 1156 | //------------------------- |
| 1157 | //fhContribVtx = new TH1F("fhContribVtx","contrib to vtx",200,0,200); |
| 1158 | //fOutputList->Add(fhContribVtx); |
| 1159 | //------------------------- |
| 1160 | //fhContribVtxAccept = new TH1F("fhContribVtxAccept","contrib to vtx after cut",200,0,200); |
| 1161 | //fOutputList->Add(fhContribVtxAccept); |
| 1162 | //------------------------- |
| 9d6abaad |
1163 | fhDphiTriggerJetMinBias = new TH2F("fhDphiTriggerJetMinBias","Deltaphi trig-jet",50,0,100, 100, -0.5*TMath::Pi(),1.5*TMath::Pi()); |
| 1164 | fOutputList->Add(fhDphiTriggerJetMinBias); |
| 1165 | |
| 1166 | fhDphiTriggerJetCent20 = (TH2F*) fhDphiTriggerJetMinBias->Clone("fhDphiTriggerJetCent20"); |
| 1167 | fOutputList->Add(fhDphiTriggerJetCent20); |
| 2cc66a40 |
1168 | //------------------------- |
| 9d6abaad |
1169 | |
| 2cc66a40 |
1170 | fhDphiTriggerJetAccept = new TH1F("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -0.5*TMath::Pi(),1.5*TMath::Pi()); |
| 1171 | fOutputList->Add(fhDphiTriggerJetAccept); |
| 1172 | //------------------------- |
| 1173 | fhCentrality = new TH1F("fhCentrality","Centrality",100,0,100); |
| 1174 | fOutputList->Add(fhCentrality); |
| 1175 | //------------------------- |
| 1176 | fhCentralityV0M = new TH1F("hCentralityV0M","hCentralityV0M",100,0,100); |
| 1177 | fOutputList->Add(fhCentralityV0M); |
| 1178 | //------------------------- |
| 1179 | fhCentralityV0A = new TH1F("hCentralityV0A","hCentralityV0A",100,0,100); |
| 1180 | fOutputList->Add(fhCentralityV0A); |
| 1181 | //------------------------- |
| 1182 | fhCentralityV0C = new TH1F("hCentralityV0C","hCentralityV0C",100,0,100); |
| 1183 | fOutputList->Add(fhCentralityV0C); |
| 1184 | //------------------------- |
| 1185 | fhCentralityZNA = new TH1F("hCentralityZNA","hCentralityZNA",100,0,100); |
| 1186 | fOutputList->Add(fhCentralityZNA); |
| 1187 | |
| 1188 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); |
| 1189 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); |
| 1190 | fOutputList->Add(fh1Xsec); |
| 1191 | |
| 1192 | fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); |
| 1193 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); |
| 1194 | fOutputList->Add(fh1Trials); |
| 1195 | |
| 1196 | fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",300,0,300); |
| 1197 | fOutputList->Add(fh1PtHard); |
| 1198 | |
| 9d6abaad |
1199 | fhImpactParameter = new TH1D("fhImpactParameter","impact parameter distribution from HIJING",50,0,10); |
| 1200 | fOutputList->Add(fhImpactParameter); |
| 1201 | |
| 1202 | fhImpactParameterTT = new TH1D("fhImpactParameterTT","b versus TT",50,0,10); |
| 1203 | fOutputList->Add(fhImpactParameterTT); |
| 2cc66a40 |
1204 | // =========== Switch on Sumw2 for all histos =========== |
| 1205 | for(Int_t i=0; i<fOutputList->GetEntries(); i++){ |
| 1206 | TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i)); |
| 1207 | if(h1){ |
| 1208 | h1->Sumw2(); |
| 1209 | continue; |
| 1210 | } |
| 1211 | THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i)); |
| 1212 | if(hn){ |
| 1213 | hn->Sumw2(); |
| 1214 | } |
| 1215 | } |
| 1216 | TH1::AddDirectory(oldStatus); |
| 1217 | |
| 1218 | |
| 1219 | PostData(1, fOutputList); |
| 1220 | } |
| 1221 | |
| 1222 | //________________________________________________________________________ |
| 1223 | void AliAnalysisTaskHJetSpectra::UserExec(Option_t *){ |
| 1224 | //executed in each event |
| 1225 | |
| 1226 | if(!InputEvent()){ |
| 1227 | AliError("??? Event pointer == 0 ???"); |
| 1228 | return; |
| 1229 | } |
| 1230 | |
| 1231 | //Execute only once: Get tracks, jets, background from arrays if not already given |
| 1232 | if(!fInitialized) ExecOnce(); |
| 1233 | if(fJetArray && fTrackArray && fBackgroundJetArray){ |
| 1234 | Calculate(InputEvent()); |
| 1235 | } |
| 1236 | PostData(1, fOutputList); |
| 1237 | } |
| 1238 | //________________________________________________________________________ |
| 1239 | |
| 1240 | Double_t AliAnalysisTaskHJetSpectra::RelativePhi(Double_t mphi,Double_t vphi){ |
| 1241 | //Get relative azimuthal angle of two particles -pi to pi |
| 1242 | if (vphi < -TMath::Pi()) vphi += TMath::TwoPi(); |
| 1243 | else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi(); |
| 1244 | |
| 1245 | if (mphi < -TMath::Pi()) mphi += TMath::TwoPi(); |
| 1246 | else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi(); |
| 1247 | |
| 1248 | Double_t dphi = mphi - vphi; |
| 1249 | if (dphi < -TMath::Pi()) dphi += TMath::TwoPi(); |
| 1250 | else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi(); |
| 1251 | |
| 1252 | return dphi;//dphi in [-Pi, Pi] |
| 1253 | } |
| 1254 | |
| 1255 | //________________________________________________________________________ |
| 1256 | Double_t AliAnalysisTaskHJetSpectra::EstimateBgRhoMedian(){ |
| 1257 | //Estimate background rho by means of integrating track pT outside identified jet cones |
| 1258 | Double_t rhoMedian = 0.0; |
| 1259 | |
| 1260 | //phi,eta and R2 of jets to be removed |
| 1261 | std::vector<Double_t> jphi; |
| 1262 | std::vector<Double_t> jeta; |
| 1263 | std::vector<Double_t> jRsquared; |
| 1264 | |
| 1265 | if(!fBackgroundJetArray) return 0.0; |
| 1266 | if(!fTrackArray) return 0.0; |
| 1267 | |
| 1268 | for(Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++){ |
| 1269 | AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i)); |
| 1270 | |
| 1271 | if(!backgroundJet){ |
| 1272 | AliError(Form("%s: Could not receive jet %d", GetName(), i)); |
| 1273 | continue; |
| 1274 | } |
| 1275 | if(!IsBackgroundJetInAcceptance(backgroundJet)) continue; //apply minimum pT cut on jet to be removed from bg |
| 1276 | jphi.push_back(RelativePhi(backgroundJet->Phi(),0.0)); //-pi,pi |
| 1277 | jeta.push_back(backgroundJet->Eta()); |
| 1278 | jRsquared.push_back(1.78*backgroundJet->Area()/TMath::Pi()); //1.78 = JetArea_R04/JetArea_R03 |
| 1279 | } |
| 1280 | |
| 1281 | |
| 1282 | static Double_t nOutCone[10][4]; |
| 1283 | static Double_t sumPtOutOfCone[10][4]; |
| 1284 | Double_t rndphi, rndeta; |
| 1285 | Double_t rndphishift, rndetashift; |
| 1286 | Double_t dphi, deta; |
| 1287 | Bool_t bIsInCone; |
| 1288 | |
| 1289 | |
| 1290 | for(Int_t ie=0; ie < fnEta; ie++){ |
| 1291 | for(Int_t ip=0; ip < fnPhi; ip++){ |
| 1292 | nOutCone[ip][ie] = 0.0; //initialize counter |
| 1293 | sumPtOutOfCone[ip][ie] = 0.0; |
| 1294 | } |
| 1295 | } |
| 1296 | |
| 1297 | //get area in cells out of identified jet cones |
| 1298 | if(jphi.size()==0){ //no jet to be removed from the bg => all areas have their nominal area |
| 1299 | for(Int_t ie=0; ie < fnEta; ie++){ |
| 1300 | for(Int_t ip=0; ip < fnPhi; ip++){ |
| 1301 | nOutCone[ip][ie] = fNofRndTrials; |
| 1302 | } |
| 1303 | } |
| 1304 | }else{ |
| 1305 | for(Int_t it=0; it<fNofRndTrials; it++){ |
| 1306 | |
| 1307 | rndphi = fRandom->Uniform(0, fPhiSize); |
| 1308 | rndeta = fRandom->Uniform(0, fEtaSize); |
| 1309 | |
| 1310 | for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell |
| 1311 | rndphishift = rndphi + ip*fPhiSize - TMath::Pi(); |
| 1312 | for(Int_t ie=0; ie<fnEta; ie++){ |
| 1313 | rndetashift = rndeta + ie*fEtaSize - fEtaSize; |
| 1314 | |
| 1315 | bIsInCone = 0; //tag if trial is in the jet cone |
| 1316 | for(Int_t ij=0; ij< (Int_t) jRsquared.size(); ij++){ |
| 1317 | deta = jeta[ij] - rndetashift; |
| 1318 | dphi = RelativePhi(rndphishift,jphi[ij]); |
| 1319 | if((dphi*dphi + deta*deta) < jRsquared[ij]){ |
| 1320 | bIsInCone = 1; |
| 1321 | break; |
| 1322 | } |
| 1323 | } |
| 1324 | if(!bIsInCone) nOutCone[ip][ie]++; |
| 1325 | } |
| 1326 | } |
| 1327 | } |
| 1328 | } |
| 1329 | |
| 1330 | Int_t phicell,etacell; |
| 1331 | for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){ |
| 1332 | AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip)); |
| 1333 | if(!part) continue; |
| 1334 | if(!IsTrackInAcceptance((AliVParticle*) part)) continue; |
| 1335 | |
| 1336 | bIsInCone = 0; //init |
| 1337 | for(Int_t ij=0; ij<(Int_t) jRsquared.size(); ij++){ |
| 1338 | dphi = RelativePhi(jphi[ij], part->Phi()); |
| 1339 | deta = jeta[ij] - part->Eta(); |
| 1340 | if((dphi*dphi + deta*deta) < jRsquared[ij]){ |
| 1341 | bIsInCone = 1; |
| 1342 | break; |
| 1343 | } |
| 1344 | } |
| 1345 | if(!bIsInCone){ |
| 1346 | phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize)); |
| 1347 | etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize)); |
| 1348 | sumPtOutOfCone[phicell][etacell]+= part->Pt(); |
| 1349 | } |
| 1350 | } |
| 1351 | // Calculate rho |
| 1352 | static Double_t rhoInCells[20]; |
| 1353 | Double_t relativeArea; |
| 1354 | Int_t nCells=0; |
| 1355 | Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac |
| 1356 | for(Int_t ip=0; ip<fnPhi; ip++){ |
| 1357 | for(Int_t ie=0; ie<fnEta; ie++){ |
| 1358 | relativeArea = nOutCone[ip][ie]/fNofRndTrials; |
| 1359 | |
| 1360 | bufferArea += relativeArea; |
| 1361 | bufferPt += sumPtOutOfCone[ip][ie]; |
| 1362 | if(bufferArea > fJetFreeAreaFrac){ |
| 1363 | rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea); |
| 1364 | |
| 1365 | bufferArea = 0.0; |
| 1366 | bufferPt = 0.0; |
| 1367 | nCells++; |
| 1368 | } |
| 1369 | } |
| 1370 | } |
| 1371 | |
| 1372 | if(nCells>0){ |
| 1373 | rhoMedian = TMath::Median(nCells, rhoInCells); |
| 1374 | } |
| 1375 | return rhoMedian; |
| 1376 | } |
| 1377 | |
| 1378 | //________________________________________________________________________ |
| 1379 | Double_t AliAnalysisTaskHJetSpectra::EstimateBgCone(){ |
| 1380 | //Estimate background rho by means of integrating track pT outside identified jet cones |
| 1381 | Double_t rhoPerpCone = 0.0; |
| 1382 | |
| 1383 | Double_t pTleading = -1.0; |
| 1384 | Double_t phiLeading = 1000.; |
| 1385 | Double_t etaLeading = 1000.; |
| 1386 | |
| 1387 | if(!fJetArray) return 0.0; |
| 1388 | |
| 1389 | for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){ |
| 1390 | AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij)); |
| 1391 | if(!jet){ |
| 1392 | AliError(Form("%s: Could not receive jet %d", GetName(), ij)); |
| 1393 | continue; |
| 1394 | } |
| 1395 | if(!IsSignalJetInAcceptance(jet)) continue; |
| 1396 | |
| 1397 | if(pTleading < jet->Pt()){ |
| 1398 | pTleading = jet->Pt(); |
| 1399 | phiLeading = jet->Phi(); |
| 1400 | etaLeading = jet->Eta(); |
| 1401 | } |
| 1402 | } |
| 1403 | if(pTleading < 0.0) return 0.0; |
| 1404 | |
| 1405 | Double_t phileftcone = phiLeading + TMath::Pi()/2; |
| 1406 | Double_t phirightcone = phiLeading - TMath::Pi()/2; |
| 1407 | |
| 1408 | /* Double_t dp, de; |
| 1409 | |
| 1410 | for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){ |
| 1411 | |
| 1412 | AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip)); |
| 1413 | if(!part) continue; |
| 1414 | if(!IsTrackInAcceptance((AliVParticle*) part)) continue; |
| 1415 | |
| 1416 | |
| 1417 | dp = RelativePhi(phileftcone, part->Phi()); |
| 1418 | de = etaLeading - part->Eta(); |
| 1419 | if( dp*dp + de*de < fRandConeRadiusSquared ) rhoPerpCone += part->Pt(); |
| 1420 | |
| 1421 | dp = RelativePhi(phirightcone, part->Phi()); |
| 1422 | if( dp*dp + de*de < fRandConeRadiusSquared) rhoPerpCone += part->Pt(); |
| 1423 | |
| 1424 | }*/ |
| 1425 | |
| 1426 | rhoPerpCone += GetConePt(etaLeading, phileftcone, fRandConeRadius); |
| 1427 | rhoPerpCone += GetConePt(etaLeading, phirightcone, fRandConeRadius); |
| 1428 | |
| 1429 | //normalize total pT by two times cone are |
| 1430 | rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fRandConeRadiusSquared); |
| 1431 | |
| 1432 | |
| 1433 | |
| 1434 | return rhoPerpCone; |
| 1435 | } |
| 9d6abaad |
1436 | //________________________________________________________________________ |
| 1437 | Bool_t AliAnalysisTaskHJetSpectra::DistantCones(Double_t phi1, Double_t eta1, Double_t r1, Double_t phi2, Double_t eta2, Double_t r2){ |
| 1438 | //checks if the two cones are farther away than the sum of their radii |
| 2cc66a40 |
1439 | |
| 9d6abaad |
1440 | Double_t dphi = RelativePhi(phi1,phi2); |
| 1441 | Double_t deta = eta1-eta2; |
| 1442 | Double_t d = r1+r2; |
| 1443 | if( dphi*dphi + deta*deta < d*d ) return kFALSE; |
| 1444 | |
| 1445 | return kTRUE; |
| 1446 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame