]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/totEt/AliAnalysisEtReconstructed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Operating mode changes: Regular/Light/SuperLight Cascade output storage for Pb-Pb
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtReconstructed.cxx
CommitLineData
cf6522d1 1//_________________________________________________________________________
2// Utility Class for transverse energy studies
3// Base class for ESD analysis
4// - reconstruction output
5// implementation file
6//
7//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8//_________________________________________________________________________
2fbf38ac 9
10#include "AliAnalysisEtReconstructed.h"
11#include "AliAnalysisEtCuts.h"
12#include "AliESDtrack.h"
ba136eb4 13#include "AliEMCALTrack.h"
2fbf38ac 14#include "AliESDCaloCluster.h"
15#include "TVector3.h"
ba136eb4 16#include "TGeoGlobalMagField.h"
17#include "AliMagF.h"
2fbf38ac 18#include "AliVEvent.h"
19#include "AliESDEvent.h"
b5821c13 20#include "AliESDtrackCuts.h"
2fbf38ac 21#include "AliVParticle.h"
87efb15c 22#include "TDatabasePDG.h"
23#include "TList.h"
b5821c13 24#include "AliESDpid.h"
2fbf38ac 25#include <iostream>
09fcb185 26#include "TH3F.h"
2fbf38ac 27#include "TH2F.h"
ef647350 28#include "TH2I.h"
29#include "TH1I.h"
30#include "TFile.h"
964c8159 31#include "AliAnalysisHadEtCorrections.h"
ef647350 32#include "AliAnalysisEtSelector.h"
0f6416f3 33#include "AliLog.h"
e9da35da 34#include "AliCentrality.h"
ef647350 35#include "AliPHOSGeoUtils.h"
36#include "AliPHOSGeometry.h"
d3ce32b8 37#include "AliAnalysisEtRecEffCorrection.h"
9a365626 38#include "AliESDpid.h"
0f6416f3 39
2fbf38ac 40
16abb579 41using namespace std;
42
43ClassImp(AliAnalysisEtReconstructed);
44
45
2fbf38ac 46AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
f61cec2f 47 AliAnalysisEt()
3919b0d0 48 ,fQATree(0)
49 ,fMakeQATree(0)
50 ,fClusterMultiplicity(0)
51 ,fTrackMultiplicity(0)
52 ,fEventID(0)
f61cec2f 53 ,fCorrections(0)
54 ,fPidCut(0)
fb116385 55 ,nChargedHadronsMeasured(0)
56 ,nChargedHadronsTotal(0)
f61cec2f 57 ,fHistChargedPionEnergyDeposit(0)
58 ,fHistProtonEnergyDeposit(0)
59 ,fHistAntiProtonEnergyDeposit(0)
60 ,fHistChargedKaonEnergyDeposit(0)
61 ,fHistMuonEnergyDeposit(0)
62 ,fHistRemovedEnergy(0)
63 ,fGeomCorrection(1.0)
b2c10007 64 ,fEMinCorrection(1.0/0.687)
f61cec2f 65 ,fRecEffCorrection(1.0)
43dd5a38 66 ,fClusterPositionAccepted(0)
67 ,fClusterPositionAll(0)
68 ,fClusterPositionAcceptedEnergy(0)
69 ,fClusterPositionAllEnergy(0)
53302bfe 70 ,fClusterEnergy(0)
9a365626 71 ,fClusterEnergyCent(0)
d386205e 72 ,fClusterEnergyModifiedTrackMatchesCent(0)
9a365626 73 ,fClusterEnergyCentMatched(0)
74 ,fClusterEnergyCentNotMatched(0)
53302bfe 75 ,fClusterEt(0)
f61cec2f 76 ,fHistChargedEnergyRemoved(0)
77 ,fHistNeutralEnergyRemoved(0)
78 ,fHistGammaEnergyAdded(0)
ac610b08 79 ,fHistMatchedTracksEvspTvsCent(0)
80 ,fHistMatchedTracksEvspTvsCentEffCorr(0)
459e9c44 81 ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
3af79c0f 82 ,fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr(0)
5881f036 83 ,fHistMatchedTracksEvspTvsCentEffTMCorr500MeV(0)
6a152780 84 ,fHistFoundHadronsvsCent(0)
85 ,fHistNotFoundHadronsvsCent(0)
86 ,fHistFoundHadronsEtvsCent(0)
87 ,fHistNotFoundHadronsEtvsCent(0)
43dd5a38 88 ,fHistFoundHadronsvsCent500MeV(0)
89 ,fHistNotFoundHadronsvsCent500MeV(0)
90 ,fHistFoundHadronsEtvsCent500MeV(0)
91 ,fHistNotFoundHadronsEtvsCent500MeV(0)
d3ce32b8 92 ,fHistNominalRawEt(0)
93 ,fHistNominalNonLinHighEt(0)
94 ,fHistNominalNonLinLowEt(0)
95 ,fHistNominalEffHighEt(0)
96 ,fHistNominalEffLowEt(0)
43dd5a38 97 ,fHistTotRawEtEffCorr(0)
3e9c52ca 98 ,fHistTotRawEt(0)
43dd5a38 99 ,fHistTotRawEtEffCorr500MeV(0)
100 ,fHistTotAllRawEt(0)
101 ,fHistTotAllRawEtEffCorr(0)
eefad699 102 ,fHistNUsedClusters(0)
9a365626 103 ,fHistNClustersPhosVsEmcal(0)
104 ,fHistClusterSizeVsCent(0)
105 ,fHistMatchedClusterSizeVsCent(0)
106 ,fHistTotAllRawEtVsTotalPt(0)
107 ,fHistTotAllRawEtVsTotalPtVsCent(0)
108 ,fHistTotMatchedRawEtVsTotalPtVsCent(0)
109 ,fHistPIDProtonsTrackMatchedDepositedVsNch(0)
110 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNch(0)
0fcbed20 111 ,fHistPIDProtonsTrackMatchedDepositedVsNcl(0)
112 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNcl(0)
9a365626 113 ,fHistPiKPTrackMatchedDepositedVsNch(0)
fb116385 114 //,
115 ,fHistPIDProtonsTrackMatchedDepositedVsNchNoEff(0)
116 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff(0)
117 ,fHistPIDProtonsTrackMatchedDepositedVsNclNoEff(0)
118 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff(0)
119 ,fHistPiKPTrackMatchedDepositedVsNchNoEff(0)
9a365626 120 ,fHistCentVsNchVsNclReco(0)
fb116385 121 ,fHistRawSignalReco(0)
122 ,fHistEffCorrSignalReco(0)
02d47689 123 ,fHistRecoRCorrVsPtVsCent(0)
2fbf38ac 124{
125
126}
127
e9da35da 128AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
f61cec2f 129{//destructor
3919b0d0 130 if(fMakeQATree){
131 //fQATree->Clear();
132 delete fQATree;
133 }
e9da35da 134 delete fCorrections;
ef647350 135 delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
136 delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
137 delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
138 delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
d0c22dcc 139 delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
140
141 delete fHistRemovedEnergy; // removed energy
43dd5a38 142 delete fClusterPositionAccepted;
143 delete fClusterPositionAll;
144 delete fClusterPositionAcceptedEnergy;
145 delete fClusterPositionAllEnergy;
53302bfe 146 delete fClusterEnergy;
9a365626 147 delete fClusterEnergyCent;
d386205e 148 delete fClusterEnergyModifiedTrackMatchesCent;
9a365626 149 delete fClusterEnergyCentMatched;
150 delete fClusterEnergyCentNotMatched;
53302bfe 151 delete fClusterEt;
311c6540 152 delete fHistChargedEnergyRemoved;
153 delete fHistNeutralEnergyRemoved;
154 delete fHistGammaEnergyAdded;
ac610b08 155 delete fHistMatchedTracksEvspTvsCent;
156 delete fHistMatchedTracksEvspTvsCentEffCorr;
459e9c44 157 delete fHistMatchedTracksEvspTvsCentEffTMCorr;
3af79c0f 158 delete fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr;
5881f036 159 delete fHistMatchedTracksEvspTvsCentEffTMCorr500MeV;
6a152780 160 delete fHistFoundHadronsvsCent;
161 delete fHistNotFoundHadronsvsCent;
162 delete fHistFoundHadronsEtvsCent;
163 delete fHistNotFoundHadronsEtvsCent;
43dd5a38 164 delete fHistFoundHadronsvsCent500MeV;
165 delete fHistNotFoundHadronsvsCent500MeV;
166 delete fHistFoundHadronsEtvsCent500MeV;
167 delete fHistNotFoundHadronsEtvsCent500MeV;
d3ce32b8 168 delete fHistNominalRawEt;
169 delete fHistNominalNonLinHighEt;
170 delete fHistNominalNonLinLowEt;
171 delete fHistNominalEffHighEt;
172 delete fHistNominalEffLowEt;
43dd5a38 173 delete fHistTotRawEtEffCorr;
3e9c52ca 174 delete fHistTotRawEt;
43dd5a38 175 delete fHistTotAllRawEt;
176 delete fHistTotAllRawEtEffCorr;
177 delete fHistTotRawEtEffCorr500MeV;
eefad699 178 delete fHistNUsedClusters;
9a365626 179 delete fHistNClustersPhosVsEmcal;
180 delete fHistClusterSizeVsCent;
181 delete fHistMatchedClusterSizeVsCent;
182 delete fHistTotAllRawEtVsTotalPt;
183 delete fHistTotAllRawEtVsTotalPtVsCent;
184 delete fHistTotMatchedRawEtVsTotalPtVsCent;
185 delete fHistPIDProtonsTrackMatchedDepositedVsNch;
186 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNch;
0fcbed20 187 delete fHistPIDProtonsTrackMatchedDepositedVsNcl;
188 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNcl;
9a365626 189 delete fHistPiKPTrackMatchedDepositedVsNch;
fb116385 190 delete fHistPIDProtonsTrackMatchedDepositedVsNchNoEff;
191 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff;
192 delete fHistPIDProtonsTrackMatchedDepositedVsNclNoEff;
193 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff;
194 delete fHistPiKPTrackMatchedDepositedVsNchNoEff;
9a365626 195 delete fHistCentVsNchVsNclReco;
fb116385 196 delete fHistRawSignalReco;
197 delete fHistEffCorrSignalReco;
02d47689 198 delete fHistRecoRCorrVsPtVsCent;
cf6522d1 199}
200
201Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
e9da35da 202{
f61cec2f 203
204 //AliAnalysisEt::AnalyseEvent(ev);
e9da35da 205 // analyse ESD event
2fbf38ac 206 ResetEventValues();
e9da35da 207 if (!ev) {
208 AliFatal("ERROR: Event does not exist");
209 return 0;
210 }
211
2fbf38ac 212 AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
e9da35da 213 if (!event) {
214 AliFatal("ERROR: ESD Event does not exist");
215 return 0;
ec956c46 216 }
2aab9269 217 if(!fSelector){
218 AliFatal("ERROR: fSelector does not exist");
219 return 0;
220 }
ef647350 221 fSelector->SetEvent(event);
f61cec2f 222
743ce29f 223 Int_t cent = -1;
d3ce32b8 224 fCentrality = event->GetCentrality();
f6b36c54 225 if (fCentrality && cent)
743ce29f 226 {
d3ce32b8 227 cent = fCentrality->GetCentralityClass5("V0M");
228 fCentClass = fCentrality->GetCentralityClass5("V0M");
743ce29f 229 }
230
9a365626 231
232 //for PID
3af79c0f 233 //AliESDpid *pID = new AliESDpid();
234 //pID->MakePID(event);
9a365626 235 Float_t etPIDProtons = 0.0;
236 Float_t etPIDAntiProtons = 0.0;
237 Float_t etPiKPMatched = 0.0;
fb116385 238 Float_t etPIDProtonsNoEff = 0.0;
239 Float_t etPIDAntiProtonsNoEff = 0.0;
240 Float_t etPiKPMatchedNoEff = 0.0;
9a365626 241 Float_t multiplicity = fEsdtrackCutsTPC->GetReferenceMultiplicity(event,kTRUE);
3919b0d0 242 fTrackMultiplicity = multiplicity;
243 fEventID = event->GetPeriodNumber();//This is not the event id
9a365626 244
245 Float_t totalMatchedPt = 0.0;
246 Float_t totalPt = 0.0;
247 TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
248 Int_t nGoodTracks = list->GetEntries();
249 for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++){
250 AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
251 if (!track)
252 {
253 Printf("ERROR: Could not get track %d", iTrack);
254 continue;
255 }
256 else{
257 totalPt +=track->Pt();
3af79c0f 258 //pID->MakeITSPID(track);
9a365626 259
260
261 }
262 }
263
5881f036 264 //TRefArray *caloClusters = fSelector->GetClusters();//just gets the correct set of clusters - does not apply any cuts
265 //Float_t fClusterMult = caloClusters->GetEntries();
09fcb185 266
d3ce32b8 267 Float_t nominalRawEt = 0;
5881f036 268 Float_t totEt500MeV = 0;
d3ce32b8 269 Float_t nonlinHighRawEt = 0;
270 Float_t nonlinLowRawEt = 0;
271 Float_t effHighRawEt = 0;
272 Float_t effLowRawEt = 0;
3e9c52ca 273 Float_t uncorrEt = 0;
f2a956c2 274 Float_t rawSignal = 0;
275 Float_t effCorrSignal = 0;
d3ce32b8 276
fb116385 277 nChargedHadronsMeasured = 0.0;
278 nChargedHadronsTotal = 0.0;
6a152780 279 Float_t nChargedHadronsEtMeasured = 0.0;
280 Float_t nChargedHadronsEtTotal = 0.0;
43dd5a38 281 Float_t nChargedHadronsMeasured500MeV = 0.0;
282 Float_t nChargedHadronsTotal500MeV = 0.0;
283 Float_t nChargedHadronsEtMeasured500MeV = 0.0;
284 Float_t nChargedHadronsEtTotal500MeV = 0.0;
285 Float_t fTotAllRawEt = 0.0;
286 Float_t fTotRawEt = 0.0;
f809f93b 287 Float_t fTotRawEtEffCorr = 0.0;
43dd5a38 288 Float_t fTotAllRawEtEffCorr = 0.0;
9a365626 289 Int_t nPhosClusters = 0;
290 Int_t nEmcalClusters = 0;
d3ce32b8 291
eefad699 292 Int_t nUsedClusters = 0;
fb116385 293
294 TRefArray *caloClusters = fSelector->GetClusters();
295 Int_t nCluster = caloClusters->GetEntries();
3919b0d0 296 fClusterMultiplicity = nCluster;
297 //if we are making the QA tree and the cluster multiplicity (for PHOS) is less than expected, fill the QA tree so we know what event it was
298 if(fMakeQATree && fClusterMultiplicity < 5e-3*fTrackMultiplicity-1.5) fQATree->Fill();
fb116385 299
300 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
2fbf38ac 301 {
fb116385 302 AliESDCaloCluster* cluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
2fbf38ac 303 if (!cluster)
304 {
e9da35da 305 AliError(Form("ERROR: Could not get cluster %d", iCluster));
2fbf38ac 306 continue;
307 }
eefad699 308// if(!fSelector->CutGeometricalAcceptance(*cluster)){
309// Float_t pos[3];
310// cluster->GetPosition(pos);
311// TVector3 cp(pos);
312// cout<<"rejected cluster phi "<<cp.Phi()*180.0/TMath::Pi()<<" eta "<<cp.Eta()<<endl;
313// }
ef647350 314 int x = 0;
eefad699 315 fCutFlow->Fill(x++);//fills 0
9a365626 316 if(cluster->IsEMCAL()) nEmcalClusters++;
317 else nPhosClusters++;
86e7d5db 318 if(!fSelector->IsDetectorCluster(*cluster)) continue;
eefad699 319 fCutFlow->Fill(x++);//fills 1
86e7d5db 320 if(!fSelector->PassMinEnergyCut(*cluster)) continue;
eefad699 321 fCutFlow->Fill(x++);//fills 2
86e7d5db 322 if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
eefad699 323 fCutFlow->Fill(x++);//fills 3
5881f036 324 if (!fSelector->CutGeometricalAcceptance(*cluster)) continue;
325 //fCutFlow->Fill(x++);
2fbf38ac 326 Float_t pos[3];
e9da35da 327
2fbf38ac 328 cluster->GetPosition(pos);
e9da35da 329 TVector3 cp(pos);
43dd5a38 330 fClusterPositionAll->Fill(cp.Phi(), cp.PseudoRapidity());
02d47689 331 Float_t fReconstructedE = cluster->E();
6280e578 332 Float_t lostEnergy = 0.0;
333 Float_t lostTrackPt = 0.0;
02d47689 334 fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
87efb15c 335
5881f036 336 //if(TMath::Abs(cp.Eta())> fCuts->fCuts->GetGeometryEmcalEtaAccCut() || cp.Phi() > fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180. || cp.Phi() > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.) continue;//Do not accept if cluster is not in the acceptance
02d47689 337 fTotAllRawEt += TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
27d85daa 338 fTotAllRawEtEffCorr +=GetCorrectionModification(*cluster,0,0,cent)* CorrectForReconstructionEfficiency(*cluster,cent);
5881f036 339
02d47689 340 fClusterEnergyCent->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
09fcb185 341 Bool_t matched = kTRUE;//default to no track matched
6280e578 342 Bool_t countasmatched = kFALSE;
96e7f328 343 Bool_t correctedcluster = kFALSE;
6280e578 344
09fcb185 345 Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
584f2478 346 matched = !(fSelector->PassTrackMatchingCut(*cluster));//PassTrackMatchingCut is false if there is a matched track
347 if(matched){//if the track match is good (, is the track good?
348 if(trackMatchedIndex < 0) matched=kFALSE;//If the index is bad, don't count it
349 if(matched){
350 AliESDtrack *track = event->GetTrack(trackMatchedIndex);
351 //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
352 matched = fEsdtrackCutsTPC->AcceptTrack(track);//If the track is bad, don't count it
02d47689 353 if(matched){//if it is still matched see if the track p was less than the energy
354 Float_t rcorr = TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
355 fHistRecoRCorrVsPtVsCent->Fill(rcorr,track->Pt(), fCentClass);
6280e578 356 if(fSelector->PassMinEnergyCut( (fReconstructedE - fsub* track->P())*TMath::Sin(cp.Theta()) ) ){//if more energy was deposited than the momentum of the track and more than one particle led to the cluster
7ee3992c 357 // if(fReconstructedE - fsub* track->P() > 0.0){
96e7f328 358 //cout<<"match corrected"<<endl;
02d47689 359 fReconstructedE = fReconstructedE - fsub* track->P();
360 matched = kFALSE;
6280e578 361 countasmatched = kTRUE;
96e7f328 362 correctedcluster = kTRUE;
6280e578 363 lostEnergy = fsub* track->P();
364 lostTrackPt = track->Pt();
d386205e 365 fClusterEnergyModifiedTrackMatchesCent->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
02d47689 366 }
96e7f328 367// else{
368// cerr<<"match passed ";
369// cerr<<"E "<<fReconstructedE<<" fsubmeanhad "<<fsub<<" p "<< track->P();
370// if(correctedcluster) cout<<" corrected";
371// cerr<<endl;
372// }
02d47689 373 }
584f2478 374 }
09fcb185 375 }
e9da35da 376
e9da35da 377
378 if (matched)
379 {
f61cec2f 380
e9da35da 381 if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
382 {
383 AliVTrack *track = event->GetTrack(trackMatchedIndex);
384 if (!track) {
385 AliError("Error: track does not exist");
386 }
387 else {
9a365626 388 totalMatchedPt +=track->Pt();
02d47689 389 fClusterEnergyCentMatched->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
9a365626 390 fHistMatchedClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
391
5881f036 392 float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cent);
584f2478 393 if(TMath::Abs(eff)<1e-5) eff = 1.0;
9a365626 394 //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
6a152780 395 nChargedHadronsMeasured++;
584f2478 396 nChargedHadronsTotal += 1/eff;
02d47689 397 Double_t effCorrEt = GetCorrectionModification(*cluster,0,0,cent) * CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cent);
398 nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
3e9c52ca 399 //One efficiency is the gamma efficiency and the other is the track matching efficiency.
02d47689 400 nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
401 //cout<<"nFound "<<1<<" nFoundTotal "<<1/eff<<" etMeas "<<TMath::Sin(cp.Theta())*fReconstructedE<<" ET total "<< 1/eff *TMath::Sin(cp.Theta())*fReconstructedE<<endl;
9a365626 402
403 Float_t nSigmaPion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
404 Float_t nSigmaProton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
405 bool isProton = (nSigmaPion>3.0 && nSigmaProton<3.0 && track->Pt()<0.9);
406 //cout<<"NSigmaProton "<<nSigmaProton<<endl;
407 etPiKPMatched += effCorrEt;
02d47689 408 etPiKPMatchedNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 409 if(isProton){
410 if(track->Charge()>0){
411 etPIDProtons += effCorrEt;
02d47689 412 etPIDProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 413 }
414 else{
02d47689 415 etPIDAntiProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 416 etPIDAntiProtons += effCorrEt;
417 }
418 }
02d47689 419 if(TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE>0.5){
43dd5a38 420 nChargedHadronsMeasured500MeV++;
421 nChargedHadronsTotal500MeV += 1/eff;
02d47689 422 nChargedHadronsEtMeasured500MeV+= TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
423 nChargedHadronsEtTotal500MeV+= 1/eff *TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
43dd5a38 424 }
5881f036 425 cluster->GetPosition(pos);
426 TVector3 p2(pos);
02d47689 427 uncorrEt += TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
96e7f328 428 if(correctedcluster || fReconstructedE <fsubmeanhade* track->P() ){//if more energy was deposited than the momentum of the track and more than one particle led to the cluster and the corrected energy is greater than zero
6280e578 429 fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
430 fHistMatchedTracksEvspTvsCentEffCorr->Fill(track->P(),effCorrEt,cent);
431 //Weighed by the number of tracks we didn't find
432 fHistMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
433 if(cent<16 && cent>11){//centralities 60-80% where false track matches are low
434 for(int cbtest = 0; cbtest<20; cbtest++){//then we calculate the deposit matched to hadrons with different centrality bins' efficiencies
435 float efftest = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cbtest);
436 if(TMath::Abs(efftest)<1e-5) efftest = 1.0;
437 Double_t effCorrEttest = GetCorrectionModification(*cluster,0,0,cent)*CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cbtest);
438 fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEttest,cbtest, (1/efftest-1) );
439 }
440 }
441 if(uncorrEt>=0.5) fHistMatchedTracksEvspTvsCentEffTMCorr500MeV->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
442 }
6280e578 443 const Double_t *pidWeights = track->PID();
444
445 Double_t maxpidweight = 0;
446 Int_t maxpid = 0;
447
448 if (pidWeights)
e9da35da 449 {
6280e578 450 for (Int_t p =0; p < AliPID::kSPECIES; p++)
e9da35da 451 {
6280e578 452 if (pidWeights[p] > maxpidweight)
e9da35da 453 {
6280e578 454 maxpidweight = pidWeights[p];
455 maxpid = p;
e9da35da 456 }
457 }
6280e578 458 if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
e9da35da 459 {
6280e578 460 fEnergyDeposited =GetCorrectionModification(*cluster,0,0,cent)* fReconstructedE;
461 fMomentumTPC = track->P();
462 fCharge = track->Charge();
463 fParticlePid = maxpid;
464 fPidProb = maxpidweight;
465 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
466 if (!esdTrack) {
467 AliError("Error: track does not exist");
468 }
469 else {
470 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
471 fDepositTree->Fill();
472 }
e9da35da 473 }
6280e578 474
475 if (maxpidweight > fPidCut)
e9da35da 476 {
6280e578 477 //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
478
479 //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
480
02d47689 481 //Float_t et = fReconstructedE * TMath::Sin(theta);
e9da35da 482 if (maxpid == AliPID::kProton)
483 {
484
485 if (track->Charge() == 1)
486 {
02d47689 487 fHistProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 488 }
489 else if (track->Charge() == -1)
490 {
02d47689 491 fHistAntiProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 492 }
493 }
494 else if (maxpid == AliPID::kPion)
495 {
02d47689 496 fHistChargedPionEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 497 }
498 else if (maxpid == AliPID::kKaon)
499 {
02d47689 500 fHistChargedKaonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 501 }
502 else if (maxpid == AliPID::kMuon)
503 {
02d47689 504 fHistMuonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 505 }
506 }
507 }
508 }
509 }
510 //continue;
ba136eb4 511 } // distance
6a152780 512 else{//these are clusters which were not track matched
f61cec2f 513 fCutFlow->Fill(x++);
514 //std::cout << x++ << std::endl;
6a152780 515
02d47689 516 //if (fReconstructedE > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
517 //if (fReconstructedE < fClusterEnergyCut) continue;
6a152780 518 cluster->GetPosition(pos);
519
6280e578 520 TVector3 p2(pos);
521 if(countasmatched){//These are tracks where we partially subtracted the energy but we subtracted some energy
522 float eff = fTmCorrections->TrackMatchingEfficiency(lostTrackPt,cent);
523 if(TMath::Abs(eff)<1e-5) eff = 1.0;
524 //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
525 nChargedHadronsMeasured++;
526 nChargedHadronsTotal += 1/eff;
527 //Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,lostEnergy,cent);
528 nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*lostEnergy;
529 //One efficiency is the gamma efficiency and the other is the track matching efficiency.
530 nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*lostEnergy;
531 }
532 fClusterPositionAccepted->Fill(p2.Phi(), p2.PseudoRapidity());
533 fClusterPositionAcceptedEnergy->Fill(p2.Phi(), p2.PseudoRapidity(),GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
534 fClusterEnergy->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
535 fClusterEnergyCentNotMatched->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
536 fHistClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
02c87d4e 537 fClusterEt->Fill(TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
6280e578 538 uncorrEt += TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
539 float myuncorrEt = TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
540 fTotRawEt += myuncorrEt;
02c87d4e 541 nUsedClusters++;
6280e578 542
543 Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cent)*GetCorrectionModification(*cluster,0,0,cent);
544 rawSignal += myuncorrEt;
545 effCorrSignal +=effCorrEt;
546 //cout<<"cluster energy "<<fReconstructedE<<" eff corr Et "<<effCorrEt<<endl;
547 fTotRawEtEffCorr += effCorrEt;
548 fTotNeutralEt += effCorrEt;
549 nominalRawEt += effCorrEt;
550 if(myuncorrEt>=0.5){
551 totEt500MeV += effCorrEt;
552 //cout<<"test "<<myuncorrEt<<"> 0.5"<<endl;
553 }
554 else{
555 //cout<<"test "<<myuncorrEt<<"< 0.5"<<endl;
556 }
557 nonlinHighRawEt += effCorrEt*GetCorrectionModification(*cluster,1,0,cent);
558 nonlinLowRawEt += effCorrEt*GetCorrectionModification(*cluster,-1,0,cent);
559 effHighRawEt += effCorrEt*GetCorrectionModification(*cluster,0,1,cent);
560 effLowRawEt += effCorrEt*GetCorrectionModification(*cluster,0,-1,cent);
561 fNeutralMultiplicity++;
e9da35da 562 }
ef647350 563 fMultiplicity++;
564 }
f61cec2f 565
eefad699 566 fHistNUsedClusters->Fill(nUsedClusters,cent);
567 // cout<<" cent "<<cent<<" total "<<nCluster<<" used "<<nUsedClusters<<endl;
fb116385 568
569 fHistRawSignalReco->Fill(rawSignal);
570 fHistEffCorrSignalReco->Fill(effCorrSignal);
571
9a365626 572 fHistNClustersPhosVsEmcal->Fill(nPhosClusters,nEmcalClusters,cent);
ef647350 573 fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
574 fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
575 fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
576 fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
f61cec2f 577
ef647350 578 fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
579 fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
2fbf38ac 580
8529c4c1 581 //Double_t removedEnergy = GetChargedContribution(cent) + GetNeutralContribution(cent) + GetGammaContribution(cent) + GetSecondaryContribution(cent);//fNeutralMultiplicity
582 Double_t removedEnergy = GetHadronContribution(cent) + GetNeutronContribution(cent) + GetKaonContribution(cent) + GetSecondaryContribution(cent);//fNeutralMultiplicity
583 //cout<<" centbin "<<cent<<" removed energy ch "<< GetHadronContribution(cent)<<" n " << GetNeutronContribution(cent) <<" kaon "<< GetKaonContribution(cent)<<" secondary "<< GetSecondaryContribution(cent);//fNeutralMultiplicity
584 //cout<<" test min et "<<fTmCorrections->GetMinEtCorrection(cent);
585 //cout<<" test neutral "<<fTmCorrections->NeutralContr(cent);
586 //cout<<" test neutron "<<fTmCorrections->NeutralContr(cent);
e9da35da 587 fHistRemovedEnergy->Fill(removedEnergy);
f61cec2f 588
459e9c44 589 fTotNeutralEtAcc = fTotNeutralEt;
f809f93b 590 //fHistTotRawEtEffCorr->Fill(fTotNeutralEt,cent);
591 fHistTotRawEtEffCorr->Fill(fTotRawEtEffCorr,cent);
43dd5a38 592 fHistTotRawEt->Fill(fTotRawEt,cent);
593 fHistTotAllRawEt->Fill(fTotAllRawEt,cent);
9a365626 594 fHistTotAllRawEtVsTotalPt->Fill(fTotAllRawEt,totalPt);
595 fHistTotAllRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalPt,cent);
596 fHistTotMatchedRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalMatchedPt,cent);
43dd5a38 597 fHistTotAllRawEtEffCorr->Fill(fTotAllRawEtEffCorr,cent);
f809f93b 598 //cout<<"fTotAllRawEtEffCorr "<<fTotAllRawEtEffCorr<<" fTotAllRawEt "<<fTotAllRawEt<<" fTotRawEtEffCorr "<<fTotRawEtEffCorr<<"("<<fTotNeutralEt<<")"<<" fTotRawEt "<<fTotRawEt<<endl;
3e9c52ca 599 //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
8529c4c1 600 //cout<<" raw "<<fTotNeutralEt<<" removed "<<removedEnergy<<" etmin "<<GetMinEtCorrection(cent)<<" final ";
601 if(GetMinEtCorrection(cent)>0) fTotNeutralEt = (fTotNeutralEt - removedEnergy)/GetMinEtCorrection(cent);
602 //cout<<fTotNeutralEt<<endl;
3e9c52ca 603 //cout<<" tot corr "<<fTotNeutralEt<<endl;
2fbf38ac 604 fTotEt = fTotChargedEt + fTotNeutralEt;
f61cec2f 605// Fill the histograms...0
2fbf38ac 606 FillHistograms();
b2c10007 607 //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
d3ce32b8 608 //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
609 fHistNominalRawEt->Fill(nominalRawEt,cent);
43dd5a38 610 fHistTotRawEtEffCorr500MeV->Fill(totEt500MeV,cent);
d3ce32b8 611 fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
612 fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
613 fHistNominalEffHighEt->Fill(effHighRawEt,cent);
614 fHistNominalEffLowEt->Fill(effLowRawEt,cent);
6a152780 615 fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
616 fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
617 fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
618 fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
9a365626 619 //cout<<"found "<<nChargedHadronsMeasured<<" total "<<nChargedHadronsTotal<<" not found "<<nChargedHadronsTotal-nChargedHadronsMeasured<<" found "<< nChargedHadronsMeasured500MeV<<" not found "<< nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV <<" total "<<nChargedHadronsTotal500MeV<<endl;
43dd5a38 620 fHistFoundHadronsvsCent500MeV->Fill(nChargedHadronsMeasured500MeV,cent);
621 fHistNotFoundHadronsvsCent500MeV->Fill(nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV,cent);
622 fHistFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtMeasured500MeV,cent);
623 fHistNotFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtTotal500MeV-nChargedHadronsEtMeasured500MeV,cent);
6a152780 624// cout<<"Number of hadrons measured: "<<nChargedHadronsMeasured<<" Estimated total number of hadrons "<<nChargedHadronsTotal<<" ET in track matched hadrons "<<
625// nChargedHadronsEtMeasured;
626// if(nChargedHadronsMeasured>0)cout<<" ("<<nChargedHadronsEtMeasured/nChargedHadronsMeasured<<") ";
627// cout<<" ET in all hadrons ";
628// cout<<nChargedHadronsEtTotal;
629// if(nChargedHadronsTotal>0) cout<<" ("<<nChargedHadronsEtTotal/nChargedHadronsTotal<<") ";
630// cout<<endl;
9a365626 631 fHistPIDProtonsTrackMatchedDepositedVsNch->Fill(etPIDProtons,multiplicity);
632 fHistPIDAntiProtonsTrackMatchedDepositedVsNch->Fill(etPIDAntiProtons,multiplicity);
fb116385 633 fHistPIDProtonsTrackMatchedDepositedVsNcl->Fill(etPIDProtons,nCluster);
634 fHistPIDAntiProtonsTrackMatchedDepositedVsNcl->Fill(etPIDAntiProtons,nCluster);
635 fHistPIDProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDProtonsNoEff,multiplicity);
636 fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDAntiProtonsNoEff,multiplicity);
637 fHistPIDProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDProtonsNoEff,nCluster);
638 fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDAntiProtonsNoEff,nCluster);
639 fHistCentVsNchVsNclReco->Fill(cent,multiplicity,nCluster);
9a365626 640 fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
fb116385 641 fHistPiKPTrackMatchedDepositedVsNchNoEff->Fill(etPiKPMatchedNoEff,multiplicity);
3af79c0f 642 //delete pID;
2fbf38ac 643 return 0;
644}
645
646bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
f61cec2f 647{ // check vertex
2fbf38ac 648
649 Float_t bxy = 999.;
650 Float_t bz = 999.;
e9da35da 651 if (!track) {
652 AliError("ERROR: no track");
653 return kFALSE;
ec956c46 654 }
1b8c3d66 655 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
e9da35da 656 if (!esdTrack) {
657 AliError("ERROR: no track");
658 return kFALSE;
1b8c3d66 659 }
660 esdTrack->GetImpactParametersTPC(bxy,bz);
ec956c46 661
2fbf38ac 662
e9da35da 663 bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
664 (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
665 (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
666 (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
667 (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
2fbf38ac 668
4998becf 669 return status;
2fbf38ac 670}
671
672void AliAnalysisEtReconstructed::Init()
f61cec2f 673{ // Init
2fbf38ac 674 AliAnalysisEt::Init();
83d0f02c 675 fPidCut = fCuts->GetReconstructedPidCut();
ba136eb4 676 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
e9da35da 677 if (!fCorrections) {
678 cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl;
0f97be4c 679 }
2fbf38ac 680}
681
682bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
f61cec2f 683{ // propagate track to detector radius
e9da35da 684 if (!track) {
685 cout<<"Warning: track empty"<<endl;
686 return kFALSE;
687 }
688 AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
689 if (!esdTrack) {
690 AliError("ERROR: no ESD track");
691 return kFALSE;
692 }
693 // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac 694
695 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
696
cf6522d1 697 // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
f61cec2f 698 return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
2fbf38ac 699}
700
87efb15c 701void AliAnalysisEtReconstructed::FillOutputList(TList* list)
f61cec2f 702{ // add some extra histograms to the ones from base class
87efb15c 703 AliAnalysisEt::FillOutputList(list);
704
3919b0d0 705 if(fMakeQATree){
706 list->Add(fQATree);
707 }
87efb15c 708 list->Add(fHistChargedPionEnergyDeposit);
709 list->Add(fHistProtonEnergyDeposit);
710 list->Add(fHistAntiProtonEnergyDeposit);
711 list->Add(fHistChargedKaonEnergyDeposit);
712 list->Add(fHistMuonEnergyDeposit);
e9da35da 713
714 list->Add(fHistRemovedEnergy);
43dd5a38 715 list->Add(fClusterPositionAccepted);
716 list->Add(fClusterPositionAll);
717 list->Add(fClusterPositionAcceptedEnergy);
718 list->Add(fClusterPositionAllEnergy);
53302bfe 719 list->Add(fClusterEnergy);
9a365626 720 list->Add(fClusterEnergyCent);
d386205e 721 list->Add(fClusterEnergyModifiedTrackMatchesCent);
9a365626 722 list->Add(fClusterEnergyCentMatched);
723 list->Add(fClusterEnergyCentNotMatched);
53302bfe 724 list->Add(fClusterEt);
f61cec2f 725
ef647350 726 list->Add(fHistChargedEnergyRemoved);
727 list->Add(fHistNeutralEnergyRemoved);
728 list->Add(fHistGammaEnergyAdded);
ac610b08 729 list->Add(fHistMatchedTracksEvspTvsCent);
730 list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
459e9c44 731 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
3af79c0f 732 list->Add(fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr);
5881f036 733 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
6a152780 734 list->Add(fHistFoundHadronsvsCent);
735 list->Add(fHistNotFoundHadronsvsCent);
736 list->Add(fHistFoundHadronsEtvsCent);
737 list->Add(fHistNotFoundHadronsEtvsCent);
43dd5a38 738 list->Add(fHistFoundHadronsvsCent500MeV);
739 list->Add(fHistNotFoundHadronsvsCent500MeV);
740 list->Add(fHistFoundHadronsEtvsCent500MeV);
741 list->Add(fHistNotFoundHadronsEtvsCent500MeV);
d3ce32b8 742 list->Add(fHistNominalRawEt);
743 list->Add(fHistNominalNonLinHighEt);
744 list->Add(fHistNominalNonLinLowEt);
745 list->Add(fHistNominalEffHighEt);
746 list->Add(fHistNominalEffLowEt);
43dd5a38 747 list->Add(fHistTotRawEtEffCorr);
748 list->Add(fHistTotRawEtEffCorr500MeV);
749 list->Add(fHistTotAllRawEtEffCorr);
3e9c52ca 750 list->Add(fHistTotRawEt);
43dd5a38 751 list->Add(fHistTotAllRawEt);
eefad699 752 list->Add(fHistNUsedClusters);
9a365626 753 list->Add(fHistNClustersPhosVsEmcal);
754 list->Add(fHistClusterSizeVsCent);
755 list->Add(fHistMatchedClusterSizeVsCent);
756 list->Add(fHistTotAllRawEtVsTotalPt);
757 list->Add(fHistTotAllRawEtVsTotalPtVsCent);
758 list->Add(fHistTotMatchedRawEtVsTotalPtVsCent);
759 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNch);
760 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNch);
0fcbed20 761 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNcl);
762 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNcl);
9a365626 763 list->Add(fHistPiKPTrackMatchedDepositedVsNch);
fb116385 764 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNchNoEff);
765 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff);
766 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNclNoEff);
767 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff);
768 list->Add(fHistPiKPTrackMatchedDepositedVsNchNoEff);
9a365626 769 list->Add(fHistCentVsNchVsNclReco);
fb116385 770 list->Add(fHistRawSignalReco);
771 list->Add(fHistEffCorrSignalReco);
02d47689 772 list->Add(fHistRecoRCorrVsPtVsCent);
87efb15c 773}
774
775void AliAnalysisEtReconstructed::CreateHistograms()
f61cec2f 776{ // add some extra histograms to the ones from base class
87efb15c 777 AliAnalysisEt::CreateHistograms();
3919b0d0 778 if(fMakeQATree){
779 fQATree = new TTree("fQATree", "fQATree");
780 fQATree->Branch("fClusterMultiplicity", &fClusterMultiplicity, "fClusterMultiplicity/I");
781 fQATree->Branch("fTrackMultiplicity", &fTrackMultiplicity, "fTrackMultiplicity/I");
782 fQATree->Branch("fEventID",&fEventID,"fEventID/I");
783 fQATree->Branch("fCentClass",&fCentClass,"fCentClass/I");
784 }
87efb15c 785
6fe6205f 786 Float_t scale = 1;//scale up histograms if EMCal 2011 so we have the right range
787 if(fDataSet==2011 && !fHistogramNameSuffix.Contains("P")){
788 scale = 2.5;
789 }
790 Int_t nbinsEt = 1000*scale;
0fa8c632 791 Double_t minEt = 0;
6fe6205f 792 Double_t maxEt = 10*scale;
0fa8c632 793
794 // possibly change histogram limits
53302bfe 795// if (fCuts) {
796// nbinsEt = fCuts->GetHistNbinsParticleEt();
797// minEt = fCuts->GetHistMinParticleEt();
798// maxEt = fCuts->GetHistMaxParticleEt();
799// }
0fa8c632 800
87efb15c 801 TString histname;
802 histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 803 fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 804 fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
805 fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
e9da35da 806
87efb15c 807 histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 808 fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 809 fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
810 fHistProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 811
87efb15c 812 histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 813 fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 814 fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
815 fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 816
87efb15c 817 histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 818 fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 819 fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
820 fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 821
87efb15c 822 histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 823 fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 824 fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
825 fHistMuonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 826
827 histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
828 fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
829 //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
830 //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
831
43dd5a38 832 histname = "fClusterPositionAccepted" + fHistogramNameSuffix;
833 fClusterPositionAccepted = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
834 fClusterPositionAccepted->SetXTitle("#phi");
835 fClusterPositionAccepted->SetYTitle("#eta");
836
837 histname = "fClusterPositionAll" + fHistogramNameSuffix;
838 fClusterPositionAll = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
839 fClusterPositionAll->SetXTitle("#phi");
840 fClusterPositionAll->SetYTitle("#eta");
841
842 histname = "fClusterPositionAcceptedEnergy" + fHistogramNameSuffix;
843 fClusterPositionAcceptedEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
844 fClusterPositionAcceptedEnergy->SetXTitle("#phi");
845 fClusterPositionAcceptedEnergy->SetYTitle("#eta");
846
847 histname = "fClusterPositionAllEnergy" + fHistogramNameSuffix;
848 fClusterPositionAllEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
849 fClusterPositionAllEnergy->SetXTitle("#phi");
850 fClusterPositionAllEnergy->SetYTitle("#eta");
ef647350 851
53302bfe 852 histname = "fClusterEnergy" + fHistogramNameSuffix;
853 fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
9a365626 854 fClusterEnergy->SetYTitle("Number of clusters");
855 fClusterEnergy->SetXTitle("Energy of cluster");
856
857 histname = "fClusterEnergyCent" + fHistogramNameSuffix;
858 fClusterEnergyCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
859 fClusterEnergyCent->SetXTitle("Energy of cluster");
860 fClusterEnergyCent->SetYTitle("Centrality Bin");
861 fClusterEnergyCent->SetZTitle("Number of clusters");
862
d386205e 863 histname = "fClusterEnergyModifiedTrackMatchesCent" + fHistogramNameSuffix;
864 fClusterEnergyModifiedTrackMatchesCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
865 fClusterEnergyModifiedTrackMatchesCent->SetXTitle("Energy of cluster");
866 fClusterEnergyModifiedTrackMatchesCent->SetYTitle("Centrality Bin");
867 fClusterEnergyModifiedTrackMatchesCent->SetZTitle("Number of clusters");
868
9a365626 869 histname = "fClusterEnergyCentMatched" + fHistogramNameSuffix;
870 fClusterEnergyCentMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
871 fClusterEnergyCentMatched->SetXTitle("Energy of cluster");
872 fClusterEnergyCentMatched->SetYTitle("Centrality Bin");
873 fClusterEnergyCentMatched->SetZTitle("Number of Clusters");
874
875 histname = "fClusterEnergyCentNotMatched" + fHistogramNameSuffix;
876 fClusterEnergyCentNotMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
877 fClusterEnergyCentNotMatched->SetXTitle("Energy of cluster");
878 fClusterEnergyCentNotMatched->SetYTitle("Centrality Bin");
879 fClusterEnergyCentNotMatched->SetZTitle("Number of clusters");
53302bfe 880
881 histname = "fClusterEt" + fHistogramNameSuffix;
02c87d4e 882 fClusterEt = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
53302bfe 883 fClusterEt->SetXTitle("Number of clusters");
884 fClusterEt->SetYTitle("E_{T} of cluster");
ef647350 885
886 histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
887 fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
888
889 histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
890 fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
891
892 histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
893 fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
e9da35da 894
5881f036 895 fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
896 fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
897 fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
3af79c0f 898 fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr", "fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
5881f036 899 fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
9a365626 900
6fe6205f 901 float max = 200*scale;
9a365626 902 if(fHistogramNameSuffix.Contains("P")){max = 100;}
903 fHistFoundHadronsvsCent = new TH2F("fHistFoundHadronsvsCent","fHistFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
904 fHistNotFoundHadronsvsCent = new TH2F("fHistNotFoundHadronsvsCent","fHistNotFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
905 fHistFoundHadronsEtvsCent = new TH2F("fHistFoundHadronsEtvsCent","fHistFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
906 fHistNotFoundHadronsEtvsCent = new TH2F("fHistNotFoundHadronsEtvsCent","fHistNotFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
907 fHistFoundHadronsvsCent500MeV = new TH2F("fHistFoundHadronsvsCent500MeV","fHistFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
908 fHistNotFoundHadronsvsCent500MeV = new TH2F("fHistNotFoundHadronsvsCent500MeV","fHistNotFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
909 fHistFoundHadronsEtvsCent500MeV = new TH2F("fHistFoundHadronsEtvsCent500MeV","fHistFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
910 fHistNotFoundHadronsEtvsCent500MeV = new TH2F("fHistNotFoundHadronsEtvsCent500MeV","fHistNotFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
5881f036 911
eefad699 912 Int_t nbinsForDistributions = 384;//2^7*3 - allows for easy rebinning
913 fHistTotRawEtEffCorr = new TH2F("fHistTotRawEtEffCorr","fHistTotRawEtEffCorr",nbinsForDistributions,0,250*scale,20,-0.5,19.5);
914 fHistTotRawEt = new TH2F("fHistTotRawEt","fHistTotRawEt",nbinsForDistributions,0,250*scale,20,-0.5,19.5);
915 fHistTotRawEtEffCorr500MeV = new TH2F("fHistTotRawEtEffCorr500MeV","fHistTotRawEtEffCorr500MeV",nbinsForDistributions,0,250*scale,20,-0.5,19.5);
916 fHistTotAllRawEt = new TH2F("fHistTotAllRawEt","fHistTotAllRawEt",nbinsForDistributions,0,250*scale,20,-0.5,19.5);
917 fHistTotAllRawEtEffCorr = new TH2F("fHistTotAllRawEtEffCorr","fHistTotAllRawEtEffCorr",nbinsForDistributions,0,250*scale,20,-0.5,19.5);
02c87d4e 918 fHistNUsedClusters = new TH2F("fHistNUsedClusters","fHistNUsedClusters",500,0,500,20,-0.5,19);
6fe6205f 919 fHistNClustersPhosVsEmcal = new TH3F("fHistNClustersPhosVsEmcal","fHistNClustersPhosVsEmcal",50,0,50,250*scale,0,250*scale,20,-0.5,19);
9a365626 920 fHistClusterSizeVsCent = new TH2F("fHistClusterSizeVsCent","fHistClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
921 fHistMatchedClusterSizeVsCent = new TH2F("fHistMatchedClusterSizeVsCent","fHistMatchedClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
6fe6205f 922 fHistTotAllRawEtVsTotalPt = new TH2F("fHistTotAllRawEtVsTotalPt","fHistTotAllRawEtVsTotalPt",125,0,250*scale,200,0,2000);
923 fHistTotAllRawEtVsTotalPtVsCent = new TH3F("fHistTotAllRawEtVsTotalPtVsCent","fHistTotAllRawEtVsTotalPtVsCent",125,0,250*scale,200,0,2000,20,-0.5,19.5);
eefad699 924 fHistTotMatchedRawEtVsTotalPtVsCent = new TH3F("fHistTotMatchedRawEtVsTotalPtVsCent","fHistTotMatchedRawEtVsTotalPtVsCent",nbinsForDistributions,0,250*scale,100,0,200,20,-0.5,19.5);
d3ce32b8 925
6fe6205f 926 maxEt = 500*scale;
d3ce32b8 927 histname = "fHistNominalRawEt" + fHistogramNameSuffix;
928 fHistNominalRawEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
929 histname = "fHistNominalNonLinHighEt" + fHistogramNameSuffix;
930 fHistNominalNonLinHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
931 histname = "fHistNominalNonLinLowEt" + fHistogramNameSuffix;
932 fHistNominalNonLinLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
933 histname = "fHistNominalEffHighEt" + fHistogramNameSuffix;
934 fHistNominalEffHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
935 histname = "fHistNominalEffLowEt" + fHistogramNameSuffix;
936 fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
937
6fe6205f 938 Float_t maxEtRange = 25*scale;
939 Float_t maxEtRangeHigh = 125*scale;
9a365626 940 Float_t minEtRange = 0;
941 Int_t nbinsMult = 100;
942 Float_t maxMult = 3000;
943 Float_t minMult = 0;
6fe6205f 944 Int_t nbinsCl = 250*scale;
945 Float_t maxCl = 500*scale;
9a365626 946 Float_t minCl = 0;
947 fHistPIDProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNch","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
948 fHistPIDAntiProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNch","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
0fcbed20 949 fHistPIDProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNcl","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
950 fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
9a365626 951 fHistPiKPTrackMatchedDepositedVsNch = new TH2F("fHistPiKPTrackMatchedDepositedVsNch","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
fb116385 952
953 fHistPIDProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNchNoEff","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
954 fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
955 fHistPIDProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNclNoEff","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
956 fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
957 fHistPiKPTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPiKPTrackMatchedDepositedVsNchNoEff","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
958
959
9a365626 960 fHistCentVsNchVsNclReco = new TH3F("fHistCentVsNchVsNclReco","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
fb116385 961
962 fHistRawSignalReco = new TH1F("fHistRawSignalReco","fHistRawSignalReco",20,-0.5,19.5);
963 fHistEffCorrSignalReco = new TH1F("fHistEffCorrSignalReco","fHistEffCorrSignalReco",20,-0.5,19.5);
02d47689 964 fHistRecoRCorrVsPtVsCent = new TH3F("fHistRecoRCorrVsPtVsCent","fHistRecoRCorrVsPtVsCent",72,0,2,50,0,10,20,-0.5,19.5);
fb116385 965
d3ce32b8 966}
5881f036 967Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)
d3ce32b8 968{
969 Float_t pos[3];
970 cluster.GetPosition(pos);
971 TVector3 cp(pos);
5881f036 972 Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),cent);
d3ce32b8 973
5881f036 974 Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,cent);
d3ce32b8 975
43dd5a38 976 cout<<"Warning: This function should not get called!"<<endl;
d3ce32b8 977 //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
978 return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
979}
87efb15c 980
5881f036 981Double_t AliAnalysisEtReconstructed::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low
d3ce32b8 982 if(nonLinCorr==0){
983 cout<<"Warning: This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
984 }
985 if(effCorr==0){
986 cout<<"Warning: This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
987 }
5881f036 988 return cluster.E()*cent;
2aab9269 989}