avoid declaring TLorentzVector per event; follow change in AliMCAnalysisUtils on...
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaClusterPileUp.cxx
CommitLineData
bc41680b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//_________________________________________________________________________
17//
18// Class for the study of Pile-up effect on
19// Calorimeter clusters.
20// Open time cuts in reader.
21//
22//-- Author: Gustavo Conesa (CNRS-LPSC-Grenoble)
23//////////////////////////////////////////////////////////////////////////////
24
25// --- ROOT system ---
26#include <TH2F.h>
27#include <TClonesArray.h>
28#include <TObjString.h>
29
30// --- Analysis system ---
31#include "AliAnaClusterPileUp.h"
32#include "AliCaloTrackReader.h"
33#include "AliFiducialCut.h"
34#include "AliVCluster.h"
35#include "AliAODEvent.h"
36#include "AliESDEvent.h"
37
38ClassImp(AliAnaClusterPileUp)
39
40//___________________________________________
41AliAnaClusterPileUp::AliAnaClusterPileUp() :
42AliAnaCaloTrackCorrBaseClass(),
0cea6003 43fNCellsCut(0),
bc41680b 44// Histograms
45fhTimePtNoCut(0), fhTimePtSPD(0),
46fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
47fhTimeNPileUpVertContributors(0),
48fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
49fhClusterMultSPDPileUp(), fhClusterMultNoPileUp(),
50fhEtaPhiBC0(0), fhEtaPhiBCPlus(0), fhEtaPhiBCMinus(0),
51fhEtaPhiBC0PileUpSPD(0),
52fhEtaPhiBCPlusPileUpSPD(0), fhEtaPhiBCMinusPileUpSPD(0),
53fhPtNPileUpSPDVtx(0), fhPtNPileUpTrkVtx(0),
54fhPtNPileUpSPDVtxTimeCut(0), fhPtNPileUpTrkVtxTimeCut(0),
55fhPtNPileUpSPDVtxTimeCut2(0), fhPtNPileUpTrkVtxTimeCut2(0)
56{
57 //default ctor
58
59 for(Int_t i = 0; i < 7; i++)
60 {
61 fhPtPileUp [i] = 0;
62 fhPtNeutralPileUp[i] = 0;
63
64 fhLambda0PileUp [i] = 0;
65 fhLambda0NeutralPileUp[i] = 0;
66
67 fhClusterEFracLongTimePileUp [i] = 0;
68
69 fhClusterCellTimePileUp [i] = 0;
70 fhClusterTimeDiffPileUp [i] = 0;
71 fhClusterTimeDiffNeutralPileUp[i] = 0;
72
73 }
74
75 for(Int_t i = 0; i < 4; i++)
76 {
77 fhClusterMultSPDPileUp[i] = 0;
78 fhClusterMultNoPileUp [i] = 0;
79 }
80
81 //Initialize parameters
82 InitParameters();
83
84}
85
86//___________________________________________
87TObjString * AliAnaClusterPileUp::GetAnalysisCuts()
88{
89 //Save parameters used for analysis
90 TString parList ; //this will be list of parameters used for this analysis.
91 const Int_t buffersize = 255;
92 char onePar[buffersize] ;
93
94 snprintf(onePar,buffersize,"--- AliAnaClusterPileUp ---\n") ;
95 parList+=onePar ;
cb67d5f1 96 snprintf(onePar,buffersize,"Calorimeter: %s\n",GetCalorimeterString().Data()) ;
bc41680b 97 parList+=onePar ;
98
99 //Get parameters set in base class.
100 //parList += GetBaseParametersList() ;
101
102 return new TObjString(parList) ;
103}
104
105//________________________________________________________________________
106TList * AliAnaClusterPileUp::GetCreateOutputObjects()
107{
108 // Create histograms to be saved in output file and
109 // store them in outputContainer
110 TList * outputContainer = new TList() ;
111 outputContainer->SetName("PhotonHistos") ;
112
113 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptimecluster = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
114 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
115 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
116 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
117 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
118
119
120 fhTimePtNoCut = new TH2F ("hTimePt_NoCut","time of cluster vs pT of clusters, no event selection", nptbins,ptmin,ptimecluster, ntimebins,timemin,timemax);
121 fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
122 fhTimePtNoCut->SetYTitle("#it{time} (ns)");
123 outputContainer->Add(fhTimePtNoCut);
124
125 fhTimePtSPD = new TH2F ("hTimePt_SPD","time of cluster vs pT of clusters, SPD Pile-up events", nptbins,ptmin,ptimecluster, ntimebins,timemin,timemax);
126 fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
127 fhTimePtSPD->SetYTitle("#it{time} (ns)");
128 outputContainer->Add(fhTimePtSPD);
129
130 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
131 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
132 fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
133 outputContainer->Add(fhTimeNPileUpVertSPD);
134
135 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
136 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
137 fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
138 outputContainer->Add(fhTimeNPileUpVertTrack);
139
140 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
141 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
142 fhTimeNPileUpVertContributors->SetXTitle("#it{time} (ns)");
143 outputContainer->Add(fhTimeNPileUpVertContributors);
144
145 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
146 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
147 fhTimePileUpMainVertexZDistance->SetXTitle("#it{time} (ns)");
148 outputContainer->Add(fhTimePileUpMainVertexZDistance);
149
150 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
151 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
152 fhTimePileUpMainVertexZDiamond->SetXTitle("#it{time} (ns)");
153 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
154
155 fhEtaPhiBC0 = new TH2F ("hEtaPhiBC0","eta-phi for clusters tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
156 fhEtaPhiBC0->SetXTitle("#eta ");
157 fhEtaPhiBC0->SetYTitle("#phi (rad)");
158 outputContainer->Add(fhEtaPhiBC0);
159
160 fhEtaPhiBCPlus = new TH2F ("hEtaPhiBCPlus","eta-phi for clusters tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
161 fhEtaPhiBCPlus->SetXTitle("#eta ");
162 fhEtaPhiBCPlus->SetYTitle("#phi (rad)");
163 outputContainer->Add(fhEtaPhiBCPlus);
164
165 fhEtaPhiBCMinus = new TH2F ("hEtaPhiBCMinus","eta-phi for clusters tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
166 fhEtaPhiBCMinus->SetXTitle("#eta ");
167 fhEtaPhiBCMinus->SetYTitle("#phi (rad)");
168 outputContainer->Add(fhEtaPhiBCMinus);
169
170 fhEtaPhiBC0PileUpSPD = new TH2F ("hEtaPhiBC0PileUpSPD","eta-phi for clusters tof corresponding to BC=0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
171 fhEtaPhiBC0PileUpSPD->SetXTitle("#eta ");
172 fhEtaPhiBC0PileUpSPD->SetYTitle("#phi (rad)");
173 outputContainer->Add(fhEtaPhiBC0PileUpSPD);
174
175 fhEtaPhiBCPlusPileUpSPD = new TH2F ("hEtaPhiBCPlusPileUpSPD","eta-phi for clusters tof corresponding to BC>0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
176 fhEtaPhiBCPlusPileUpSPD->SetXTitle("#eta ");
177 fhEtaPhiBCPlusPileUpSPD->SetYTitle("#phi (rad)");
178 outputContainer->Add(fhEtaPhiBCPlusPileUpSPD);
179
180 fhEtaPhiBCMinusPileUpSPD = new TH2F ("hEtaPhiBCMinusPileUpSPD","eta-phi for clusters tof corresponding to BC<0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
181 fhEtaPhiBCMinusPileUpSPD->SetXTitle("#eta ");
182 fhEtaPhiBCMinusPileUpSPD->SetYTitle("#phi (rad)");
183 outputContainer->Add(fhEtaPhiBCMinusPileUpSPD);
184
185 TString title[] = {"no |t diff| cut","|t diff|<20 ns","|t diff|>20 ns","|t diff|>40 ns"};
186 TString name [] = {"TDiffNoCut","TDiffSmaller20ns","TDiffLarger20ns","TDiffLarger40ns"};
187 for(Int_t i = 0; i < 4; i++)
188 {
189 fhClusterMultSPDPileUp[i] = new TH2F(Form("fhClusterMultSPDPileUp_%s", name[i].Data()),
190 Form("Number of clusters per pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
191 nptbins,ptmin,ptimecluster,100,0,100);
192 fhClusterMultSPDPileUp[i]->SetYTitle("n clusters ");
193 fhClusterMultSPDPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
194 outputContainer->Add(fhClusterMultSPDPileUp[i]) ;
195
196 fhClusterMultNoPileUp[i] = new TH2F(Form("fhClusterMultNoPileUp_%s", name[i].Data()),
197 Form("Number of clusters per non pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
198 nptbins,ptmin,ptimecluster,100,0,100);
199 fhClusterMultNoPileUp[i]->SetYTitle("n clusters ");
200 fhClusterMultNoPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
201 outputContainer->Add(fhClusterMultNoPileUp[i]) ;
202 }
203
204 fhPtNPileUpSPDVtx = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
205 nptbins,ptmin,ptimecluster,20,0,20);
206 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
207 fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
208 outputContainer->Add(fhPtNPileUpSPDVtx);
209
210 fhPtNPileUpTrkVtx = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
211 nptbins,ptmin,ptimecluster, 20,0,20 );
212 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
213 fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
214 outputContainer->Add(fhPtNPileUpTrkVtx);
215
216 fhPtNPileUpSPDVtxTimeCut = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
217 nptbins,ptmin,ptimecluster,20,0,20);
218 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
219 fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
220 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
221
222 fhPtNPileUpTrkVtxTimeCut = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
223 nptbins,ptmin,ptimecluster, 20,0,20 );
224 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
225 fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
226 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
227
228 fhPtNPileUpSPDVtxTimeCut2 = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
229 nptbins,ptmin,ptimecluster,20,0,20);
230 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
231 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
232 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
233
234 fhPtNPileUpTrkVtxTimeCut2 = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
235 nptbins,ptmin,ptimecluster, 20,0,20 );
236 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
237 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
238 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
239
240
cb67d5f1 241 TString pileUpName[] = {"SPD",kEMCAL,"SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
bc41680b 242
243 for(Int_t i = 0 ; i < 7 ; i++)
244 {
245 fhPtPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
246 Form("Cluster #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster);
247 fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
248 outputContainer->Add(fhPtPileUp[i]);
249
250 fhPtNeutralPileUp[i] = new TH1F(Form("hPtNeutralPileUp%s",pileUpName[i].Data()),
251 Form("Neutral clusters #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster);
252 fhPtNeutralPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
253 outputContainer->Add(fhPtNeutralPileUp[i]);
254
255 fhClusterEFracLongTimePileUp[i] = new TH2F(Form("hClusterEFracLongTimePileUp%s",pileUpName[i].Data()),
256 Form("Cluster E vs fraction of cluster energy from large T cells, %s Pile-Up event",pileUpName[i].Data()),
257 nptbins,ptmin,ptimecluster,200,0,1);
258 fhClusterEFracLongTimePileUp[i]->SetXTitle("#it{E} (GeV)");
259 fhClusterEFracLongTimePileUp[i]->SetYTitle("E(large time) / E");
260 outputContainer->Add(fhClusterEFracLongTimePileUp[i]);
261
262 fhClusterCellTimePileUp[i] = new TH2F(Form("hClusterCellTimePileUp%s",pileUpName[i].Data()),
263 Form("Cluster E vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
264 nptbins,ptmin,ptimecluster,ntimebins,timemin,timemax);
265 fhClusterCellTimePileUp[i]->SetXTitle("#it{E} (GeV)");
266 fhClusterCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
267 outputContainer->Add(fhClusterCellTimePileUp[i]);
268
269 fhClusterTimeDiffPileUp[i] = new TH2F(Form("hClusterTimeDiffPileUp%s",pileUpName[i].Data()),
270 Form("Cluster E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
271 nptbins,ptmin,ptimecluster,400,-200,200);
272 fhClusterTimeDiffPileUp[i]->SetXTitle("#it{E} (GeV)");
273 fhClusterTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
274 outputContainer->Add(fhClusterTimeDiffPileUp[i]);
275
276 fhClusterTimeDiffNeutralPileUp[i] = new TH2F(Form("hClusterTimeDiffNeutralPileUp%s",pileUpName[i].Data()),
277 Form("Neutral clusters E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
278 nptbins,ptmin,ptimecluster,400,-200,200);
279 fhClusterTimeDiffNeutralPileUp[i]->SetXTitle("#it{E} (GeV)");
280 fhClusterTimeDiffNeutralPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
281 outputContainer->Add(fhClusterTimeDiffNeutralPileUp[i]);
282
283 fhLambda0PileUp[i] = new TH2F(Form("hLambda0PileUp%s",pileUpName[i].Data()),
284 Form("Cluster E vs #lambda^{2}_{0} in cluster, %s Pile-Up event",pileUpName[i].Data()),
285 nptbins,ptmin,ptimecluster,ssbins,ssmin,ssmax);
286 fhLambda0PileUp[i]->SetXTitle("#it{E} (GeV)");
287 fhLambda0PileUp[i]->SetYTitle("#lambda^{2}_{0}");
288 outputContainer->Add(fhLambda0PileUp[i]);
289
290 fhLambda0NeutralPileUp[i] = new TH2F(Form("hLambda0NeutralPileUp%s",pileUpName[i].Data()),
291 Form("Neutral clusters E vs #lambda^{2}_{0}in cluster, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptimecluster,ssbins,ssmin,ssmax);
292 fhLambda0NeutralPileUp[i]->SetXTitle("#it{E} (GeV)");
293 fhLambda0NeutralPileUp[i]->SetYTitle("#lambda^{2}_{0}");
294 outputContainer->Add(fhLambda0NeutralPileUp[i]);
295
296 }
297
298 return outputContainer ;
299
300}
301
302//_______________________
303void AliAnaClusterPileUp::Init()
304{
305 //Init
306
307 //Do some checks
cb67d5f1 308 if(GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn())
bc41680b 309 AliFatal("You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
310
cb67d5f1 311 if(GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn())
bc41680b 312 AliFatal("You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
313
314 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC)
315 AliFatal("You want to use MC data in analysis but this is not possible in pile-up!!");
316
317}
318
319//____________________________________________________________________________
320void AliAnaClusterPileUp::InitParameters()
321{
322
323 //Initialize the parameters of the analysis.
324 AddToHistogramsName("AnaClusterPileUp_");
325
bc41680b 326 fNCellsCut = 2;
327}
328
329//__________________________________________________________________
330void AliAnaClusterPileUp::MakeAnalysisFillHistograms()
331{
332 // Do cluster analysis
333 // Remember to open time cuts in reader
334
335 //Select the calorimeter
336 TObjArray * pl = 0x0;
337 AliVCaloCells* cells = 0;
cb67d5f1 338 if (GetCalorimeter() == kPHOS )
bc41680b 339 {
340 pl = GetPHOSClusters();
341 cells = GetPHOSCells();
342 }
cb67d5f1 343 else if (GetCalorimeter() == kEMCAL)
bc41680b 344 {
345 pl = GetEMCALClusters();
346 cells = GetEMCALCells();
347 }
348
349 if(!pl)
350 {
cb67d5f1 351 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",GetCalorimeterString().Data());
bc41680b 352 return;
353 }
354
355 AliVEvent * event = GetReader()->GetInputEvent();
356 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
357 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
358
359 //-------------------
360 // N pile up vertices
361 Int_t nVtxSPD = -1;
362 Int_t nVtxTrk = -1;
363
364 if (esdEv)
365 {
366 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
367 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
368 }//ESD
369 else if (aodEv)
370 {
371 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
372 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
373 }//AOD
374
375 //-------------------
376 // Loop on clusters
377 Int_t nCaloClusters = pl->GetEntriesFast();
378
cb67d5f1 379 if(GetDebug() > 0) printf("AliAnaClusterPileUp::MakeAnalysisFillAOD() - input %s cluster entries %d\n", GetCalorimeterString().Data(), nCaloClusters);
bc41680b 380 //Init variables
381 TLorentzVector mom;
382 Int_t idMax = 0;
383 Float_t ptMax = 0;
384 Float_t tMax = 0;
385
386 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
387 {
388 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
389 //printf("calo %d, %f\n",icalo,calo->E());
390
391 if(!calo) continue; // it should not happen, but just in case
392
393 calo->GetMomentum(mom,GetVertex(0)) ;
394
395 Float_t ecluster = mom.E();
396 Float_t ptcluster = mom.Pt();
397 Float_t l0cluster = calo->GetM02();
398 Float_t etacluster = mom.Eta();
399 Float_t phicluster = mom.Phi();
400 if(phicluster < 0) phicluster+=TMath::TwoPi();
401 Float_t tofcluster = calo->GetTOF()*1.e9;
402
403 Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
404
405 //.......................................
406 //If too small or big energy, skip it
407 if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) continue ;
408
409 //.......................................
410 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
411
412 //.......................................
413 //Check acceptance selection
414 if(IsFiducialCutOn())
415 {
1290eee4 416 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom.Eta(),mom.Phi(),GetCalorimeter()) ;
bc41680b 417 if(! in ) continue;
418 }
419
420 // Select highest pt cluster passing the cuts
421 if(ptcluster > ptMax && tofcluster < 30)
422 {
423 ptMax = ptcluster;
424 tMax = tofcluster;
425 idMax = icalo;
426 }
427
428 //-------------------------------------
429 // Cluster timing for different pile-up
430
431 fhTimePtNoCut->Fill(ptcluster,tofcluster);
432 if(GetReader()->IsPileUpFromSPD()) fhTimePtSPD->Fill(ptcluster,tofcluster);
433
434 //----------------------------------------
435 // correlate cluster and number of vertices
436
437 fhPtNPileUpSPDVtx->Fill(ptcluster,nVtxSPD);
438 fhPtNPileUpTrkVtx->Fill(ptcluster,nVtxTrk);
439
440 if(TMath::Abs(tofcluster) < 30)
441 {
442 fhPtNPileUpSPDVtxTimeCut->Fill(ptcluster,nVtxSPD);
443 fhPtNPileUpTrkVtxTimeCut->Fill(ptcluster,nVtxTrk);
444 }
445
446 if(tofcluster < 75 && tofcluster > -30)
447 {
448 fhPtNPileUpSPDVtxTimeCut2->Fill(ptcluster,nVtxSPD);
449 fhPtNPileUpTrkVtxTimeCut2->Fill(ptcluster,nVtxTrk);
450 }
451
452 // Loop on the vertices arrays, correlate with timing
453 // only for sufficiently large cluster energy
454 if(ecluster > 8)
455 {
456 fhTimeNPileUpVertSPD ->Fill(tofcluster,nVtxSPD);
457 fhTimeNPileUpVertTrack->Fill(tofcluster,nVtxTrk);
458
459 Int_t ncont = -1;
460 Float_t z1 = -1, z2 = -1;
461 Float_t diamZ = -1;
462 for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
463 {
464 if (esdEv)
465 {
466 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
467 ncont=pv->GetNContributors();
468 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
469 z2 = pv->GetZ();
470 diamZ = esdEv->GetDiamondZ();
471 }//ESD
472 else if (aodEv)
473 {
474 AliAODVertex *pv=aodEv->GetVertex(iVert);
475 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
476 ncont=pv->GetNContributors();
477 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
478 z2=pv->GetZ();
479 diamZ = aodEv->GetDiamondZ();
480 }// AOD
481
482 Double_t distZ = TMath::Abs(z2-z1);
483 diamZ = TMath::Abs(z2-diamZ);
484
485 fhTimeNPileUpVertContributors ->Fill(tofcluster,ncont);
486 fhTimePileUpMainVertexZDistance->Fill(tofcluster,distZ);
487 fhTimePileUpMainVertexZDiamond ->Fill(tofcluster,diamZ);
488
489 }// vertex loop
490 }
491
492 //------------------------------------
493 // Eta-Phi cluster position depending on timing
494 // Continue only for BC0
495 if (tofcluster > 28)
496 {
497 fhEtaPhiBCPlus ->Fill(etacluster,phicluster);
498 if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCPlusPileUpSPD ->Fill(etacluster,phicluster);
499 continue;
500 }
501 else if (tofcluster <-28)
502 {
503 fhEtaPhiBCMinus->Fill(etacluster,phicluster);
504 if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCMinusPileUpSPD->Fill(etacluster,phicluster);
505 continue ;
506 }
507
508 //--------------------------------------
509 // Fill histograms for clusters in BC=0
510
511 fhEtaPhiBC0->Fill(etacluster,phicluster); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBC0PileUpSPD ->Fill(etacluster,phicluster);
512
513 if(GetReader()->IsPileUpFromSPD()) {fhPtPileUp[0]->Fill(ptcluster); fhLambda0PileUp[0]->Fill(ptcluster,l0cluster); }
514 if(GetReader()->IsPileUpFromEMCal()) {fhPtPileUp[1]->Fill(ptcluster); fhLambda0PileUp[1]->Fill(ptcluster,l0cluster); }
515 if(GetReader()->IsPileUpFromSPDOrEMCal()) {fhPtPileUp[2]->Fill(ptcluster); fhLambda0PileUp[2]->Fill(ptcluster,l0cluster); }
516 if(GetReader()->IsPileUpFromSPDAndEMCal()) {fhPtPileUp[3]->Fill(ptcluster); fhLambda0PileUp[3]->Fill(ptcluster,l0cluster); }
517 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) {fhPtPileUp[4]->Fill(ptcluster); fhLambda0PileUp[4]->Fill(ptcluster,l0cluster); }
518 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) {fhPtPileUp[5]->Fill(ptcluster); fhLambda0PileUp[5]->Fill(ptcluster,l0cluster); }
519 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(ptcluster); fhLambda0PileUp[6]->Fill(ptcluster,l0cluster); }
520
521 if(!matched)
522 {
523 if(GetReader()->IsPileUpFromSPD()) {fhPtNeutralPileUp[0]->Fill(ptcluster); fhLambda0NeutralPileUp[0]->Fill(ptcluster,l0cluster); }
524 if(GetReader()->IsPileUpFromEMCal()) {fhPtNeutralPileUp[1]->Fill(ptcluster); fhLambda0NeutralPileUp[1]->Fill(ptcluster,l0cluster); }
525 if(GetReader()->IsPileUpFromSPDOrEMCal()) {fhPtNeutralPileUp[2]->Fill(ptcluster); fhLambda0NeutralPileUp[2]->Fill(ptcluster,l0cluster); }
526 if(GetReader()->IsPileUpFromSPDAndEMCal()) {fhPtNeutralPileUp[3]->Fill(ptcluster); fhLambda0NeutralPileUp[3]->Fill(ptcluster,l0cluster); }
527 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) {fhPtNeutralPileUp[4]->Fill(ptcluster); fhLambda0NeutralPileUp[4]->Fill(ptcluster,l0cluster); }
528 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) {fhPtNeutralPileUp[5]->Fill(ptcluster); fhLambda0NeutralPileUp[5]->Fill(ptcluster,l0cluster); }
529 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtNeutralPileUp[6]->Fill(ptcluster); fhLambda0NeutralPileUp[6]->Fill(ptcluster,l0cluster); }
530 }
531
532 //----------------------------------------------------------------------------
533 // Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
534 // Get the fraction of the cluster energy that carries the cell with highest
535 // energy and its absId
536
537 Float_t maxCellFraction = 0.;
538 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
539
540 Float_t clusterLongTimePt = 0;
541 Float_t clusterOKTimePt = 0;
542
543 if(cells->GetCellAmplitude(absIdMax) < 0.1) continue ;
544
545 for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
546 {
547 Int_t absId = calo->GetCellsAbsId()[ipos];
548
549 if( absId == absIdMax ) continue ;
550
551 Double_t time = cells->GetCellTime(absId);
552 Float_t amp = cells->GetCellAmplitude(absId);
553 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
554 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,time,cells);
555 time*=1e9;
556
557 Float_t diff = (tofcluster-time);
558
559 if(GetReader()->IsInTimeWindow(time,amp)) clusterOKTimePt += amp;
560 else clusterLongTimePt += amp;
561
562 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
563
564 if(GetReader()->IsPileUpFromSPD())
565 {
566 fhClusterCellTimePileUp[0]->Fill(ptcluster, time);
567 fhClusterTimeDiffPileUp[0]->Fill(ptcluster, diff);
568 if(!matched) fhClusterTimeDiffNeutralPileUp[0]->Fill(ptcluster, diff);
569 }
570
571 if(GetReader()->IsPileUpFromEMCal())
572 {
573 fhClusterCellTimePileUp[1]->Fill(ptcluster, time);
574 fhClusterTimeDiffPileUp[1]->Fill(ptcluster, diff);
575 if(!matched) fhClusterTimeDiffNeutralPileUp[1]->Fill(ptcluster, diff);
576 }
577
578 if(GetReader()->IsPileUpFromSPDOrEMCal())
579 {
580 fhClusterCellTimePileUp[2]->Fill(ptcluster, time);
581 fhClusterTimeDiffPileUp[2]->Fill(ptcluster, diff);
582 if(!matched) fhClusterTimeDiffNeutralPileUp[2]->Fill(ptcluster, diff);
583 }
584
585 if(GetReader()->IsPileUpFromSPDAndEMCal())
586 {
587 fhClusterCellTimePileUp[3]->Fill(ptcluster, time);
588 fhClusterTimeDiffPileUp[3]->Fill(ptcluster, diff);
589 if(!matched) fhClusterTimeDiffNeutralPileUp[3]->Fill(ptcluster, diff);
590 }
591
592 if(GetReader()->IsPileUpFromSPDAndNotEMCal())
593 {
594 fhClusterCellTimePileUp[4]->Fill(ptcluster, time);
595 fhClusterTimeDiffPileUp[4]->Fill(ptcluster, diff);
596 if(!matched) fhClusterTimeDiffNeutralPileUp[4]->Fill(ptcluster, diff);
597 }
598
599 if(GetReader()->IsPileUpFromEMCalAndNotSPD())
600 {
601 fhClusterCellTimePileUp[5]->Fill(ptcluster, time);
602 fhClusterTimeDiffPileUp[5]->Fill(ptcluster, diff);
603 if(!matched) fhClusterTimeDiffNeutralPileUp[5]->Fill(ptcluster, diff);
604 }
605
606 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
607 {
608 fhClusterCellTimePileUp[6]->Fill(ptcluster, time);
609 fhClusterTimeDiffPileUp[6]->Fill(ptcluster, diff);
610 if(!matched) fhClusterTimeDiffNeutralPileUp[6]->Fill(ptcluster, diff);
611 }
612 }//loop
613
614 Float_t frac = 0;
615 if(clusterLongTimePt+clusterOKTimePt > 0.001)
616 frac = clusterLongTimePt/(clusterLongTimePt+clusterOKTimePt);
617 //printf("E long %f, E OK %f, Fraction large time %f, E %f\n",clusterLongTimePt,clusterOKTimePt,frac,ptcluster);
618
619 if(GetReader()->IsPileUpFromSPD()) fhClusterEFracLongTimePileUp[0]->Fill(ptcluster,frac);
620 if(GetReader()->IsPileUpFromEMCal()) fhClusterEFracLongTimePileUp[1]->Fill(ptcluster,frac);
621 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhClusterEFracLongTimePileUp[2]->Fill(ptcluster,frac);
622 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhClusterEFracLongTimePileUp[3]->Fill(ptcluster,frac);
623 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhClusterEFracLongTimePileUp[4]->Fill(ptcluster,frac);
624 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhClusterEFracLongTimePileUp[5]->Fill(ptcluster,frac);
625 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhClusterEFracLongTimePileUp[6]->Fill(ptcluster,frac);
626 }//loop
627
628 //----------------------------------------------------------------------------------------------
629 // Loop again on clusters to compare this max cluster t and the rest of the clusters, if E > 0.3
630 Int_t n20 = 0;
631 Int_t n40 = 0;
632 Int_t n = 0;
633 Int_t nOK = 0;
634
635 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
636 {
637 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
638
639 if(!calo || calo->E() < 0.3 || icalo == idMax) continue;
640
641 Float_t tdiff = TMath::Abs(tMax-calo->GetTOF()*1e9);
642 n++;
643 if(tdiff < 25) nOK++;
644 else
645 {
646 n20++;
647 if(tdiff > 40 ) n40++;
648 }
649 }
650
651 // Check pile-up and fill histograms depending on the different cluster multiplicities
652 if(GetReader()->IsPileUpFromSPD())
653 {
654 fhClusterMultSPDPileUp[0]->Fill(ptMax,n );
655 fhClusterMultSPDPileUp[1]->Fill(ptMax,nOK);
656 fhClusterMultSPDPileUp[2]->Fill(ptMax,n20);
657 fhClusterMultSPDPileUp[3]->Fill(ptMax,n40);
658 }
659 else
660 {
661 fhClusterMultNoPileUp[0]->Fill(ptMax,n );
662 fhClusterMultNoPileUp[1]->Fill(ptMax,nOK);
663 fhClusterMultNoPileUp[2]->Fill(ptMax,n20);
664 fhClusterMultNoPileUp[3]->Fill(ptMax,n40);
665 }
666
667
668 if(GetDebug() > 1) printf("AliAnaClusterPileUp::MakeAnalysisFillHistograms() End fill histograms\n");
669
670}
671
672
673
674//__________________________________________________________________
675void AliAnaClusterPileUp::Print(const Option_t * opt) const
676{
677 //Print some relevant parameters set for the analysis
678
679 if(! opt)
680 return;
681
682 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
683 AliAnaCaloTrackCorrBaseClass::Print(" ");
684
cb67d5f1 685 printf("Calorimeter = %s\n", GetCalorimeterString().Data()) ;
bc41680b 686 printf(" \n") ;
687
688}