]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/DstCommonMacros.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge remote-tracking branch 'origin/master' into mergingFlat
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / macrosJPSI / RaaPbPb2010 / DstCommonMacros.C
CommitLineData
eac83c86 1//
2// Common macros used to analyze dst trees
3// Author: Ionut-Cristian Arsene, 2012/03/04
4// email: i.c.arsene@gsi.de
5//
6#include <iostream>
7#include <fstream>
8using namespace std;
9
10#include <TObjArray.h>
11#include <TChain.h>
12#include <TMath.h>
13#include <TFile.h>
14#include <TDirectory.h>
15#include <THashList.h>
16#include <TH1F.h>
17#include <TH2F.h>
18#include <TH3F.h>
19#include <TProfile.h>
20#include <TProfile2D.h>
21#include <TIterator.h>
22#include <TKey.h>
23#include <TDirectory.h>
24#include <TAxis.h>
25
26#ifndef ALICORRELATIONREDUCEDEVENT_H
27#include "AliCorrelationReducedEvent.h"
28#endif
29
30namespace DstCommonMacros {
31
32 enum ParticleId {
33 kUnknown = -1,
34 kElectron = 0,
35 kPion,
36 kKaon,
37 kProton
38 };
39
40 enum Variables {
41 kNothing = -1,
42 // Event wise variables
43 kRunNo = 0, // run number
44 kBC, // bunch crossing
45 kTriggerMask, // trigger mask
46 kOfflineTrigger, // offline trigger
47 kOfflineTriggerFired, // offline trigger fired
48 kOfflineTriggerFired2, // offline trigger if fired, -1 if not fired
49 kIsPhysicsSelection, // physics selection
50 kVtxX, // vtx X
51 kVtxY, // vtx Y
52 kVtxZ, // vtx Z
53 kVtxXtpc, // vtx X from tpc
54 kVtxYtpc, // vtx Y from tpc
55 kVtxZtpc, // vtx Z from tpc
56 kCentVZERO, // centrality from VZERO
57 kCentSPD, // centrality from SPD
58 kCentTPC, // centrality from TPC
59 kCentZDC, // centrality from ZDC
60 kCentQuality, // centrality quality
61 kNV0total, // total number of V0s in the esd
62 kNV0selected, // number of V0s selected
63 kNdielectrons, // number of dielectron pairs
64 kNpairsSelected, // number of selected dielectron pairs per event
65 kNtracksTotal, // total number of tracks
66 kNtracksSelected, // number of selected tracks
67 kSPDntracklets, // SPD number of tracklets in |eta|<1.0
68 kEventMixingId, // Id of the event mixing category
69 // VZERO event plane related variables
70 kVZEROAemptyChannels, // Number of empty VZERO channels in A side
71 kVZEROCemptyChannels, // Number of empty VZERO channels in C side
72 kVZEROChannelMult, // VZERO multiplicity per channel
73 kVZEROChannelEta = kVZEROChannelMult+64, // pseudo-rapidity of a VZERO channel
74 kVZEROQvecX = kVZEROChannelEta+64, // Q-vector components for harmonics 1-6 and
75 kVZEROQvecY = kVZEROQvecX+6*3, // 6- n-harmonics; 3- A,C and A&C options
76 kVZERORP = kVZEROQvecY+6*3, // VZERO reaction plane from A,C and A&C sides (harmonics 1-6)
77 kVZERORPres = kVZERORP+6*3, // VZERO reaction plane resolution (sqrt(n*(RPa-RPc)))
78 kVZEROXaXc = kVZERORPres+6, // correlations for the components of the Q vector
79 kVZEROXaYa = kVZEROXaXc+6,
80 kVZEROXaYc = kVZEROXaYa+6,
81 kVZEROYaXc = kVZEROXaYc+6,
82 kVZEROYaYc = kVZEROYaXc+6,
83 kVZEROXcYc = kVZEROYaYc+6,
84 kVZEROdeltaRPac = kVZEROXcYc+6, // Psi_VZEROA-Psi_VZEROC
85 kVZEROdeltaRPa = kVZEROdeltaRPac+6, // Psi_VZEROA6-Psi_VZEROA5, Psi_VZEROA7-Psi_VZEROA5, Psi_VZEROA6-Psi_VZEROA5
86 kVZEROdeltaRPc = kVZEROdeltaRPa+6*2, // Psi_VZEROA2-Psi_VZEROA1, Psi_VZEROA3-Psi_VZEROA1, Psi_VZEROA4-Psi_VZEROA1
87 kVZEROflowV2TPC = kVZEROdeltaRPc+6*2, // vzero v2 using TPC event plane
88 // TPC event plane variables
89 kTPCQvecX = kVZEROflowV2TPC+64, // TPC Q-vector components for harmonics 1-6
90 kTPCQvecY = kTPCQvecX+6,
91 kTPCRP = kTPCQvecY+6, // Event plane using TPC
92 kTPCRPres = kTPCRP+6, // Event plane resolution variables sqrt(n*(RPtpc-RPvzeroa)),sqrt(n*(RPtpc-RPvzeroc))
93 // Correlations between TPC and VZERO event planes
94 kRPXtpcXvzeroa = kTPCRPres+6*2,
95 kRPXtpcXvzeroc = kRPXtpcXvzeroa+6,
96 kRPYtpcYvzeroa = kRPXtpcXvzeroc+6,
97 kRPYtpcYvzeroc = kRPYtpcYvzeroa+6,
98 kRPXtpcYvzeroa = kRPYtpcYvzeroc+6,
99 kRPXtpcYvzeroc = kRPXtpcYvzeroa+6,
100 kRPYtpcXvzeroa = kRPXtpcYvzeroc+6,
101 kRPYtpcXvzeroc = kRPYtpcXvzeroa+6,
102 kRPdeltaVZEROAtpc = kRPYtpcXvzeroc+6,
103 kRPdeltaVZEROCtpc = kRPdeltaVZEROAtpc+6,
104 kNEventVars = kRPdeltaVZEROCtpc+6, // number of event variables
105 // Pair variables --------------------------------------
106 kMass=kNEventVars,
107 kCandidateId,
108 kPairType, // 0 ++; 1 +-; 2 --
109 kPairPt,
110 kPairPx,
111 kPairPy,
112 kPairPz,
113 kPairP,
114 kPairRap,
115 kPairEta,
116 kPairTheta,
117 kPairPhi,
118 kPairLxy,
119 kPairOpeningAngle,
120 kPairCosNPhi, // cos (n*phi)
121 kPairSinNPhi = kPairCosNPhi+6, // sin (n*phi)
122 kPairDeltaPhiVZEROFlowVn = kPairSinNPhi+6, // phi - Psi_{VZERO}
123 kPairDeltaPhiTPCFlowVn = kPairDeltaPhiVZEROFlowVn+6*3, // phi - Psi_{TPC}
124 kPairVZEROFlowVn = kPairDeltaPhiTPCFlowVn+6, // vn{Psi_{n,VZERO}}
125 kPairTPCFlowVn = kPairVZEROFlowVn+6*3, // vn{Psi_{n,TPC}}
126 kPairVZEROFlowV3Psi2 = kPairTPCFlowVn+6, // v3{Psi_{2,VZERO}}
127 kPairTPCFlowV3Psi2 = kPairVZEROFlowV3Psi2+3, // v3{Psi_{2,TPC}}
128 // Track variables -------------------------------------
129 kPt=kPairTPCFlowV3Psi2+1,
130 kP,
131 kTheta,
132 kPhi,
133 kEta,
134 kRap,
135 kPtTPC,
136 kPhiTPC,
137 kEtaTPC,
138 kPin,
139 kTrackDeltaPhiVZEROFlowVn,
140 kTrackVZEROFlowVn = kTrackDeltaPhiVZEROFlowVn+6*2,
141 kDcaXY = kTrackVZEROFlowVn+6*2,
142 kDcaZ,
143 kITSncls,
144 kITSsignal,
145 kTPCncls,
146 kTPCcrossedRows,
147 kTPCnclsIter1,
148 kTPCnclsF,
149 kTPCnclsRatio,
150 kTPCsignal,
151 kTPCnSig,
152 kTOFbeta=kTPCnSig+4,
153 kTOFnSig,
154 kTRDntracklets=kTOFnSig+4,
155 kTRDntrackletsPID,
156 kTRDpidProbabilities,
157 kEMCALmatchedEnergy=kTRDpidProbabilities+2,
158 kEMCALmatchedEOverP,
159 // Calorimeter cluster variables --------------------------------------
160 kEMCALclusterEnergy,
161 kEMCALclusterDx,
162 kEMCALclusterDz,
163 kEMCALdetector, // 0 - EMCAL; 1 - PHOS
164 // Tracking flags -----------------------------------------------------
165 kTrackingFlag,
166 // Correlation variables ----------------------------------------------
167 kDeltaPhi,
168 kDeltaTheta,
169 kDeltaEta,
170 kNVars
171 };
172
173
174 Bool_t gUsedVars[kNVars] = {kFALSE};
175
176 // tracking flags as in AliESDtrack.h
177 // NOTE: check consistency with aliroot
178 enum TrackingFlags {
179 kITSin=0,
180 kITSout,
181 kITSrefit,
182 kITSpid,
183 kTPCin,
184 kTPCout,
185 kTPCrefit,
186 kTPCpid,
187 kTRDin,
188 kTRDout,
189 kTRDrefit,
190 kTRDpid,
191 kTOFin,
192 kTOFout,
193 kTOFrefit,
194 kTOFpid,
195 kTOFmismatch,
196 kHMPIDout,
197 kHMPIDpid,
198 kEMCALmatch,
199 kPHOSmatch,
200 kTRDbackup,
201 kTRDStop,
202 kESDpid,
203 kTIME,
204 kGlobalMerge,
205 kITSpureSA,
206 kMultInV0,
207 kMultSec,
208 kTRDnPlanes,
209 kEMCALNoMatch,
210 kNTrackingFlags
211 };
212
213 const Char_t* gkTrackingFlagNames[kNTrackingFlags] = {
214 "kITSin", "kITSout", "kITSrefit", "kITSpid",
215 "kTPCin", "kTPCout", "kTPCrefit", "kTPCpid",
216 "kTRDin", "kTRDout", "kTRDrefit", "kTRDpid",
217 "kTOFin", "kTOFout", "kTOFrefit", "kTOFpid", "kTOFmismatch",
218 "kHMPIDout", "kHMPIDpid",
219 "kEMCALmatch", "kPHOSmatch",
220 "kTRDbackup", "kTRDStop",
221 "kESDpid", "kTIME", "kGlobalMerge",
222 "kITSpureSA",
223 "kMultInV0",
224 "kMultSec",
225 "kTRDnPlanes",
226 "kEMCALNoMatch"
227 };
228
229 // offline triggers as defined in AliVEvent.h
230 // NOTE: Check consistency with updates in aliroot!!!
231 enum EOfflineTriggerTypes {
232 kMB = BIT(0), // Minimum bias trigger, i.e. interaction trigger, offline SPD or V0 selection
233 kINT7 = BIT(1), // V0AND trigger, offline V0 selection
234 kMUON = BIT(2), // Muon trigger, offline SPD or V0 selection
235 kHighMult = BIT(3), // High-multiplicity trigger (threshold defined online), offline SPD or V0 selection
236 kEMC1 = BIT(4), // EMCAL trigger
237 kCINT5 = BIT(5), // Minimum bias trigger without SPD. i.e. interaction trigger, offline V0 selection
238 kCMUS5 = BIT(6), // Muon trigger, offline V0 selection
239 kMUSPB = BIT(6), // idem for PbPb
240 kMUSH7 = BIT(7), // Muon trigger: high pt, single muon, offline V0 selection, CINT7 suite
241 kMUSHPB = BIT(7), // idem for PbPb
242 kMUL7 = BIT(8), // Muon trigger: like sign dimuon, offline V0 selection, CINT7 suite
243 kMuonLikePB = BIT(8), // idem for PbPb
244 kMUU7 = BIT(9), // Muon trigger, unlike sign dimuon, offline V0 selection, CINT7 suite
245 kMuonUnlikePB = BIT(9), // idem for PbPb
246 kEMC7 = BIT(10), // EMCAL trigger, CINT7 suite
247 kMUS7 = BIT(11), // Muon trigger: low pt, single muon, offline V0 selection, CINT7 suite
248 kPHI1 = BIT(12), // PHOS trigger, CINT1 suite
249 kPHI7 = BIT(13), // PHOS trigger, CINT7 suite
250 kPHOSPb = BIT(13), // idem for PbPb
251 kEMCEJE = BIT(14), // EMCAL jet patch trigger
252 kEMCEGA = BIT(15), // EMCAL gamma trigger
253 kCentral = BIT(16), // PbPb central collision trigger
254 kSemiCentral = BIT(17), // PbPb semicentral collision trigger
255 kDG5 = BIT(18), // Double gap diffractive
256 kZED = BIT(19), // ZDC electromagnetic dissociation
257 kUserDefined = BIT(27), // Set when custom trigger classes are set in AliPhysicsSelection, offline SPD or V0 selection
258 // Bits 28 and above are reserved for FLAGS
259 kFastOnly = BIT(30), // The fast cluster fired. This bit is set in to addition another trigger bit, e.g. kMB
260 kAny = 0xffffffff, // to accept any trigger
261 kAnyINT = kMB | kINT7 | kCINT5 // to accept any interaction (aka minimum bias) trigger
262 };
263
264 const Char_t* gkOfflineTriggerNames[64] = {
265 "MB", "INT7", "MUON", "HighMult", "EMC1", "CINT5", "CMUS5/MUSPB", "MUSH7/MUSHPB",
266 "MUL7/MuonLikePB", "MUU7/MuonUnlikePB", "EMC7", "MUS7", "PHI1", "PHI7/PHOSPb", "EMCEJE", "EMCEGA",
267 "Central", "SemiCentral", "DG5", "ZED", "N/A", "N/A", "N/A", "N/A",
268 "N/A", "N/A", "N/A", "UserDefined", "N/A", "N/A", "FastOnly", "N/A",
269 "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A",
270 "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A",
271 "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A",
272 "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A"
273 };
274
275 enum ITSLayerMap {
276 kITSfirst = 1,
277 kITSsecond = 2,
278 kITSthird = 4,
279 kITSfourth = 8,
280 kITSfifth = 16,
281 kITSsixth = 32
282 };
283
284 // radii of VZERO channels centers (in cm)
285 const Double_t gkVZEROChannelRadii[64] = {6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567,
286 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977,
287 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504,
288 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031,
289 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347,
290 10.685, 10.685, 10.685, 10.685, 10.685, 10.685, 10.685, 10.685,
291 18.116, 18.116, 18.116, 18.116, 18.116, 18.116, 18.116, 18.116,
292 31.84, 31.84, 31.84, 31.84, 31.84, 31.84, 31.84, 31.84};
293 const Double_t gkVZEROAz = 340.0; // cm
294 const Double_t gkVZEROCz = 90.0; // cm
295 const Double_t gkVZEROminMult = 0.5; // minimum VZERO channel multiplicity
296
297 // Pointer to the current event
298 AliCorrelationReducedEvent* gCurrentEvent = 0x0;
299 AliCorrelationReducedEventFriend* gCurrentEventFriend = 0x0;
300
301 // Event mixing variables
302 Int_t gEMCategories=0;
303 TString* gEMCategoryNames;
304
305 TObjArray* gHistLists=0x0; // main histogram list for the current running process
306 TDirectoryFile* gHistListsOld=0x0; // main directory for a standard tree analysis output (used for calibration, plotting etc.)
307 TFile* gHistFile = 0x0; // pointer to a TFile opened for reading
308
309 TProfile2D* gVzeroAvMult[64] = {0x0}; // pointer to the array of average VZERO multiplicity per channel
310 TProfile2D* gQvecCentering[AliCorrelationReducedEventFriend::kNdetectors][fgkNMaxHarmonics][2] = {{{0x0}}}; // pointer to the array of Qvec centering histograms
311
312 // pt range for J/psi's
313 const Float_t gkJpsiPtCut[2] = {0.0, 20.0};
314
315 // Function prototypes
316 void WriteOutput(TFile* saveFile);
317 //void DefineHistograms(const Char_t* histClasses);
318 TChain* GetChain(const Char_t* filename, Int_t howMany, Int_t offset, Long64_t& entries, TChain* friendChain=0x0, const Char_t* friendChainFile=0x0);
319 void FillEventInfo(AliCorrelationReducedEvent* event, Float_t* values, AliCorrelationReducedEventFriend* eventF=0x0);
320 void FillEventOfflineTriggers(UShort_t triggerBit, Float_t* values);
321 void FillTrackingFlag(AliCorrelationReducedTrack* track, UShort_t flag, Float_t* values);
322 void FillTrackInfo(AliCorrelationReducedTrack* p, Float_t* values);
323 void FillPairInfo(AliCorrelationReducedPair* p, Float_t* values);
324 void FillPairInfo(AliCorrelationReducedTrack* t1, AliCorrelationReducedTrack* t2, Int_t type, Float_t* values);
325 void FillCorrelationInfo(AliCorrelationReducedPair* p, AliCorrelationReducedTrack* t, Float_t* values);
326 void FillCaloClusterInfo(AliCorrelationReducedCaloCluster* cl, Float_t* values);
327 void FillHistClass(const Char_t* className, Float_t* values);
328 void DoEventMixing(TList* list1, TList* list2, Float_t* values, Int_t mixingType, const Char_t* histClass);
329 void EventMixingPairTracks(TList* pairs, TList* tracks, Float_t* values);
330 void EventMixingResonanceLegs(TList* posLegs, TList* negLegs, Float_t* values, Int_t type, const Char_t* histClass);
331 Double_t DeltaPhi(Double_t phi1, Double_t phi2); // calculate delta phi in the (-pi,+pi) interval
332 Bool_t IsPairSelectedEM(Float_t* values); // pair selection used in the mixed event
333 void AddHistClass(const Char_t* histClass);
334 Int_t ValidateHistogramName(THashList* hList, const Char_t* name);
335 void AddHistogram(const Char_t* histClass,
336 const Char_t* name, const Char_t* title, Bool_t isProfile,
337 Int_t nXbins, Double_t xmin, Double_t xmax, Int_t varX,
338 Int_t nYbins=0, Double_t ymin=0, Double_t ymax=0, Int_t varY=kNothing,
339 Int_t nZbins=0, Double_t zmin=0, Double_t zmax=0, Int_t varZ=kNothing,
340 const Char_t* xLabels="", const Char_t* yLabels="", const Char_t* zLabels="");
341 void AddHistogram(const Char_t* histClass,
342 const Char_t* name, const Char_t* title, Bool_t isProfile,
343 Int_t nXbins, Double_t* xbins, Int_t varX,
344 Int_t nYbins=0, Double_t* ybins=0x0, Int_t varY=kNothing,
345 Int_t nZbins=0, Double_t* zbins=0x0, Int_t varZ=kNothing,
346 const Char_t* xLabels="", const Char_t* yLabels="", const Char_t* zLabels="");
347 void MakeAxisLabels(TAxis* ax, const Char_t* labels);
348 void InitFile(const Char_t* filename); // open an old output filename
349 void CloseFile();
350 TObject* GetHistogram(const Char_t* listname, const Char_t* hname); // get a histogram from an old output
351
352} // end declarations for namespace DstCommonMacros
353
354//__________________________________________________________________
355void DstCommonMacros::FillEventInfo(AliCorrelationReducedEvent* event, Float_t* values, AliCorrelationReducedEventFriend* eventF/*=0x0*/) {
356 //
357 // fill event wise info
358 //
359 values[kRunNo] = event->RunNo();
360 values[kBC] = event->BC();
361 values[kTriggerMask] = event->TriggerMask();
362 values[kIsPhysicsSelection] = (event->IsPhysicsSelection() ? 1.0 : 0.0);
363 values[kVtxX] = event->Vertex(0);
364 values[kVtxY] = event->Vertex(1);
365 values[kVtxZ] = event->Vertex(2);
366 values[kVtxXtpc] = event->VertexTPC(0);
367 values[kVtxYtpc] = event->VertexTPC(1);
368 values[kVtxZtpc] = event->VertexTPC(2);
369 values[kCentVZERO] = event->CentralityVZERO();
370 values[kCentSPD] = event->CentralitySPD();
371 values[kCentTPC] = event->CentralityTPC();
372 values[kCentZDC] = event->CentralityZEMvsZDC();
373 values[kCentQuality] = event->CentralityQuality();
374 values[kNV0total] = event->NV0CandidatesTotal();
375 values[kNV0selected] = event->NV0Candidates();
376 values[kNdielectrons] = event->NDielectrons();
377 values[kNtracksTotal] = event->NTracksTotal();
378 values[kNtracksSelected] = event->NTracks();
379 values[kSPDntracklets] = event->SPDntracklets();
380 values[kVZEROAemptyChannels] = 0;
381 values[kVZEROCemptyChannels] = 0;
382 for(Int_t ich=0;ich<64;++ich) gUsedVars[kVZEROChannelMult+ich] = kTRUE;
383 Float_t theta=0.0;
384 for(Int_t ich=0;ich<64;++ich) {
385 if(gUsedVars[kVZEROChannelMult+ich]) {
386 values[kVZEROChannelMult+ich] = event->MultChannelVZERO(ich);
387 if(!gVzeroAvMult[0] && values[kVZEROChannelMult+ich]<gkVZEROminMult) {
388 gUsedVars[kVZEROChannelMult+ich] = kFALSE; // will not be filled in histograms by the histogram manager
389 if(ich<32) values[kVZEROCemptyChannels] += 1;
390 else values[kVZEROAemptyChannels] += 1;
391 }
392 }
393 if(gUsedVars[kVZEROChannelEta+ich]) {
394 if(ich<32) theta = TMath::ATan(gkVZEROChannelRadii[ich]/(gkVZEROCz-values[kVtxZ]));
395 else theta = TMath::Pi()-TMath::ATan(gkVZEROChannelRadii[ich]/(gkVZEROAz-values[kVtxZ]));
396 values[kVZEROChannelEta+ich] = -1.0*TMath::Log(TMath::Tan(theta/2.0));
397 }
398 }
399
400 if(eventF) {
401 for(Int_t ih=0; ih<6; ++ih) {
402 // VZERO event plane variables
403 values[kVZEROQvecX+2*6+ih] = 0.0;
404 values[kVZEROQvecY+2*6+ih] = 0.0;
405 values[kVZERORP +2*6+ih] = 0.0;
406 for(Int_t iVZEROside=0; iVZEROside<2; ++iVZEROside) {
407 values[kVZEROQvecX+iVZEROside*6+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1);
408 values[kVZEROQvecY+iVZEROside*6+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1);
409 values[kVZERORP +iVZEROside*6+ih] = eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1);
410 values[kVZEROQvecX+2*6 +ih] += values[kVZEROQvecX+iVZEROside*6+ih];
411 values[kVZEROQvecY+2*6 +ih] += values[kVZEROQvecY+iVZEROside*6+ih];
412 // cos(n(EPtpc-EPvzero A/C))
413 values[kTPCRPres+iVZEROside*6+ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kTPC, ih+1), eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1));
414 values[kTPCRPres+iVZEROside*6+ih] = TMath::Cos(values[kTPCRPres+iVZEROside*6+ih]*(ih+1));
415 }
416 values[kVZERORP +2*6+ih] = TMath::ATan2(values[kVZEROQvecY+2*6+ih],values[kVZEROQvecX+2*6+ih]);
417 // cos (n*(psi_A-psi_C))
418 values[kVZERORPres + ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA, ih+1),
419 eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROC, ih+1));
420 values[kVZERORPres + ih] = TMath::Cos(values[kVZERORPres + ih]*(ih+1));
421 // Qx,Qy correlations for VZERO
422 values[kVZEROXaXc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1);
423 values[kVZEROXaYa+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1);
424 values[kVZEROXaYc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1);
425 values[kVZEROYaXc+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1);
426 values[kVZEROYaYc+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1);
427 values[kVZEROXcYc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1);
428 // Psi_A - Psi_C
429 values[kVZEROdeltaRPac+ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA, ih+1),
430 eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROC, ih+1));
431
432 // TPC event plane
433 values[kTPCQvecX+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kTPC, ih+1);
434 values[kTPCQvecY+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kTPC, ih+1);
435 values[kTPCRP +ih] = eventF->EventPlane(AliCorrelationReducedEventFriend::kTPC, ih+1);
436 // TPC VZERO Q-vector correlations
437 values[kRPXtpcXvzeroa+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecX+ih];
438 values[kRPXtpcXvzeroc+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecX+6+ih];
439 values[kRPYtpcYvzeroa+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecY+ih];
440 values[kRPYtpcYvzeroc+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecY+6+ih];
441 values[kRPXtpcYvzeroa+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecY+ih];
442 values[kRPXtpcYvzeroc+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecY+6+ih];
443 values[kRPYtpcXvzeroa+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecX+ih];
444 values[kRPYtpcXvzeroc+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecX+6+ih];
445 // Psi_TPC - Psi_VZERO A/C
446 values[kRPdeltaVZEROAtpc+ih] = DeltaPhi(values[kVZERORP+0*6+ih], values[kTPCRP+ih]);
447 values[kRPdeltaVZEROCtpc+ih] = DeltaPhi(values[kVZERORP+1*6+ih], values[kTPCRP+ih]);
448 } // end loop over harmonics
449
450 // VZERO v2 using TPC event plane
451 Double_t vzeroChannelPhi[8] = {0.3927, 1.1781, 1.9635, 2.7489, -2.7489, -1.9635, -1.1781, -0.3927};
452 for(Int_t ich=0; ich<64; ++ich) {
453 if(gUsedVars[kVZEROflowV2TPC]) {
454 values[kVZEROflowV2TPC+ich] = values[kVZEROChannelMult+ich]*TMath::Cos(2.0*DeltaPhi(vzeroChannelPhi[ich%8],values[kTPCRP+1]));
455 }
456 }
457 } // end if (eventF)
458}
459
460
461//_________________________________________________________________
462void DstCommonMacros::FillTrackingFlag(AliCorrelationReducedTrack* track, UShort_t flag, Float_t* values) {
463 //
464 // fill the tracking flag
465 //
466 values[kTrackingFlag] = -1;
467 if(track->CheckTrackStatus(flag)) values[kTrackingFlag] = flag;
468}
469
470
471//_________________________________________________________________
472void DstCommonMacros::FillEventOfflineTriggers(UShort_t triggerBit, Float_t* values) {
473 //
474 // fill the trigger bit input
475 //
476 if(triggerBit>=64) return;
477 if(!gCurrentEvent) return;
478 ULong64_t trigger = BIT(0);
479 values[kOfflineTrigger] = triggerBit;
480 values[kOfflineTriggerFired] = (gCurrentEvent->TriggerMask()&(trigger<<triggerBit) ? 1.0 : 0.0);
481 values[kOfflineTriggerFired2] = (values[kOfflineTriggerFired]>0.01 ? triggerBit : -1.0);
482}
483
484
485//_________________________________________________________________
486void DstCommonMacros::FillTrackInfo(AliCorrelationReducedTrack* p, Float_t* values) {
487 //
488 // fill track information
489 //
490 values[kPt] = p->Pt();
491 values[kPtTPC] = p->PtTPC();
492 if(gUsedVars[kP]) values[kP] = p->P();
493 if(gUsedVars[kTheta]) values[kTheta] = p->Theta();
494 values[kPhi] = p->Phi();
495 values[kPhiTPC] = p->PhiTPC();
496 values[kEta] = p->Eta();
497 values[kEtaTPC] = p->EtaTPC();
498 values[kPin] = p->Pin();
499 values[kDcaXY] = p->DCAxy();
500 values[kDcaZ] = p->DCAz();
501
502 if(gUsedVars[kITSncls]) values[kITSncls] = p->ITSncls();
503 values[kITSsignal] = p->ITSsignal();
504
505 values[kTPCncls] = p->TPCncls();
506 if(gUsedVars[kTPCnclsRatio])
507 values[kTPCnclsRatio] = (p->TPCFindableNcls()>0 ? Float_t(p->TPCncls())/Float_t(p->TPCFindableNcls()) : 0.0);
508 values[kTPCnclsIter1] = p->TPCnclsIter1();
509 values[kTPCnclsF] = p->TPCFindableNcls();
510 values[kTPCsignal] = p->TPCsignal();
511
512 values[kTOFbeta] = p->TOFbeta();
513 for(Int_t specie=kElectron; specie<=kProton; ++specie) {
514 values[kTPCnSig+specie] = p->TPCnSig(specie);
515 values[kTOFnSig+specie] = p->TOFnSig(specie);
516 }
517 values[kTRDpidProbabilities] = p->TRDpid(0);
518 values[kTRDpidProbabilities+1] = p->TRDpid(1);
519
520 values[kTRDntracklets] = p->TRDntracklets(0);
521 values[kTRDntrackletsPID] = p->TRDntracklets(1);
522
523 if(gUsedVars[kEMCALmatchedEnergy] || gUsedVars[kEMCALmatchedEOverP]) {
524 AliCorrelationReducedCaloCluster* cluster = gCurrentEvent->GetCaloCluster(p->CaloClusterId());
525 values[kEMCALmatchedEnergy] = (cluster ? cluster->Energy() : -999.0);
526 Float_t mom = p->P();
527 values[kEMCALmatchedEnergy] = (TMath::Abs(mom)>1.e-8 && cluster ? values[kEMCALmatchedEOverP]/mom : -999.0);
528 }
529
530 // Fill track flow variables
531 for(Int_t iVZEROside=0; iVZEROside<2; ++iVZEROside) {
532 for(Int_t ih=0; ih<6; ++ih) {
533 if(gUsedVars[kTrackVZEROFlowVn+iVZEROside*6+ih] || gUsedVars[kTrackDeltaPhiVZEROFlowVn+iVZEROside*6+ih]) {
534 values[kTrackVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPhi],values[kVZERORP+iVZEROside*6+ih])*(ih+1));
535 values[kTrackDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPhi],values[kVZERORP+iVZEROside*6+ih]);
536 }
537 }
538 }
539}
540
541
542//_________________________________________________________________
543void DstCommonMacros::FillCaloClusterInfo(AliCorrelationReducedCaloCluster* cl, Float_t* values) {
544 //
545 // Fill calorimeter cluster information
546 //
547 values[kEMCALclusterEnergy] = cl->Energy();
548 values[kEMCALclusterDx] = cl->Dx();
549 values[kEMCALclusterDz] = cl->Dz();
550 values[kEMCALdetector] = (cl->IsEMCAL() ? AliCorrelationReducedCaloCluster::kEMCAL : AliCorrelationReducedCaloCluster::kPHOS);
551}
552
553
554//_________________________________________________________________
555void DstCommonMacros::FillPairInfo(AliCorrelationReducedPair* p, Float_t* values) {
556 //
557 // fill pair information
558 //
559
560 values[kCandidateId] = p->CandidateId();
561 values[kPairType] = p->PairType();
562 if(gUsedVars[kMass]) {
563 values[kMass] = p->Mass();
564 if(p->CandidateId()==AliCorrelationReducedPair::kLambda0ToPPi) values[kMass] = p->Mass(1);
565 if(p->CandidateId()==AliCorrelationReducedPair::kALambda0ToPPi) values[kMass] = p->Mass(2);
566 }
567 values[kPairPt] = p->Pt();
568 if(gUsedVars[kPairP]) values[kPairP] = p->P();
569 if(gUsedVars[kPairPx]) values[kPairPx] = p->Px();
570 if(gUsedVars[kPairPy]) values[kPairPy] = p->Py();
571 if(gUsedVars[kPairPz]) values[kPairPz] = p->Pz();
572 values[kPairEta] = p->Eta();
573 if(gUsedVars[kPairRap]) values[kPairRap] = p->Rapidity();
574 values[kPairPhi] = p->Phi();
575 values[kPairLxy] = p->Lxy();
576 values[kPairOpeningAngle] = p->OpeningAngle();
577 if(gUsedVars[kPairTheta]) values[kPairTheta] = p->Theta();
578
579 // Flow variables
580 // cos(n*phi), sin(n*phi)
581 for(Int_t ih=0; ih<6; ++ih) {
582 if(gUsedVars[kPairCosNPhi+ih]) values[kPairCosNPhi+ih] = TMath::Cos(values[kPairPhi]*(1.0+ih));
583 if(gUsedVars[kPairSinNPhi+ih]) values[kPairSinNPhi+ih] = TMath::Sin(values[kPairPhi]*(1.0+ih));
584 }
585 // VZERO
586 for(Int_t iVZEROside=0; iVZEROside<3; ++iVZEROside) {
587 // v3 using VZERO Psi_2
588 if(gUsedVars[kPairVZEROFlowV3Psi2+iVZEROside])
589 values[kPairVZEROFlowV3Psi2+iVZEROside] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kVZERORP+iVZEROside*6+1]));
590 for(Int_t ih=0; ih<6; ++ih) {
591 // vn using VZERO Psi_n
592 if(gUsedVars[kPairVZEROFlowVn+iVZEROside*6+ih])
593 values[kPairVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih])*(ih+1));
594 // phi - Psi_n
595 if(gUsedVars[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih])
596 values[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih]);
597 }
598 }
599 // TPC
600 // v3 using Psi_2
601 if(gUsedVars[kPairTPCFlowV3Psi2])
602 values[kPairTPCFlowV3Psi2] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kTPCRP+1]));
603 for(Int_t ih=0; ih<6; ++ih) {
604 // vn using Psi_n
605 if(gUsedVars[kPairTPCFlowVn+ih])
606 values[kPairTPCFlowVn+ih] = TMath::Cos(DeltaPhi(values[kPairPhi],values[kTPCRP+ih])*(ih+1));
607 if(gUsedVars[kPairDeltaPhiTPCFlowVn+ih])
608 values[kPairDeltaPhiTPCFlowVn+ih] = DeltaPhi(values[kPairPhi], values[kTPCRP+ih]);
609 }
610}
611
612
613//_________________________________________________________________
614void DstCommonMacros::FillPairInfo(AliCorrelationReducedTrack* t1, AliCorrelationReducedTrack* t2, Int_t type, Float_t* values) {
615 //
616 // fill pair information from 2 tracks
617 //
618 // type - Parameter encoding the resonance type
619 // This is needed for making a mass assumption on the legs
620 //
621 if(gUsedVars[kPairType]) {
622 if(t1->Charge()*t2->Charge()<0) values[kPairType] = 1;
623 else if(t1->Charge()>0) values[kPairType] = 0;
624 else values[kPairType] = 2;
625 }
626 Float_t kMass1 = 0.0;
627 Float_t kMass2 = 0.0;
628 switch (type) {
629 case AliCorrelationReducedPair::kK0sToPiPi :
630 kMass1 = 0.13957; kMass2 = 0.13957;
631 break;
632 case AliCorrelationReducedPair::kPhiToKK :
633 kMass1 = 0.493677; kMass2 = 0.493677;
634 break;
635 case AliCorrelationReducedPair::kLambda0ToPPi :
636 kMass1 = 0.938272; kMass2 = 0.13957;
637 break;
638 case AliCorrelationReducedPair::kALambda0ToPPi :
639 kMass1 = 0.13957; kMass2 = 0.938272;
640 break;
641 case AliCorrelationReducedPair::kJpsiToEE :
642 kMass1 = 0.000511; kMass2 = 0.000511;
643 break;
644 default :
645 break;
646 }
647
648 if(gUsedVars[kMass]) {
649 values[kMass] = kMass1*kMass1+kMass2*kMass2 +
650 2.0*(TMath::Sqrt(kMass1*kMass1+t1->P()*t1->P())*TMath::Sqrt(kMass2*kMass2+t2->P()*t2->P()) -
651 t1->Px()*t2->Px() - t1->Py()*t2->Py() - t1->Pz()*t2->Pz());
652 if(values[kMass]<0.0) {
653 cout << "FillPairInfo(track, track, type, values): Warning: Very small squared mass found. "
654 << " Could be negative due to resolution of Float_t so it will be set to a small positive value." << endl;
655 cout << " mass2: " << values[kMass] << endl;
656 cout << "p1(p,x,y,z): " << t1->P() << ", " << t1->Px() << ", " << t1->Py() << ", " << t1->Pz() << endl;
657 cout << "p2(p,x,y,z): " << t2->P() << ", " << t2->Px() << ", " << t2->Py() << ", " << t2->Pz() << endl;
658 values[kMass] = 0.0;
659 }
660 else
661 values[kMass] = TMath::Sqrt(values[kMass]);
662 }
663
664 if(gUsedVars[kPairPt]) values[kPairPt] = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) +
665 (t1->Py()+t2->Py())*(t1->Py()+t2->Py()));
666 if(gUsedVars[kPairP]) values[kPairP] = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) +
667 (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) +
668 (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz()));
669 if(gUsedVars[kPairEta]) {
670 Float_t p = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) +
671 (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) +
672 (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz()));
673 values[kPairEta] = p-t1->Pz()-t2->Pz();
674 values[kPairEta] = (TMath::Abs(values[kPairEta])>1.0e-8 ? (p+t1->Pz()+t2->Pz())/values[kPairEta] : 0.0);
675 values[kPairEta] = (values[kPairEta]>1.0e-8 ? 0.5*TMath::Log(values[kPairEta]) : -999.);
676 }
677 if(gUsedVars[kPairRap]) {
678 Float_t mass = kMass1*kMass1+kMass2*kMass2 +
679 2.0*(TMath::Sqrt(kMass1*kMass1+t1->P()*t1->P())*TMath::Sqrt(kMass2*kMass2+t2->P()*t2->P()) -
680 t1->Px()*t2->Px() - t1->Py()*t2->Py() - t1->Pz()*t2->Pz());
681 if(mass<0.0) {
682 cout << "Negative squared mass (Float_t resolution). Setting mass to zero" << endl;
683 mass = 0.0;
684 }
685 else mass = TMath::Sqrt(mass);
686 Float_t e = TMath::Sqrt(mass*mass+
687 (t1->Px()+t2->Px())*(t1->Px()+t2->Px()) +
688 (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) +
689 (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz()));
690 values[kPairRap] = 0.5*TMath::Log((e+t1->Pz()+t2->Pz())/(e-t1->Pz()-t2->Pz()));
691 }
692 if(gUsedVars[kPairPhi]) {
693 values[kPairPhi] = TMath::ATan2(t1->Py()+t2->Py(),t1->Px()+t2->Px());
694 if(values[kPairPhi]<0.0) values[kPairPhi] = 2.0*TMath::Pi() + values[kPairPhi];
695 }
696 if(gUsedVars[kPairTheta]) values[kPairTheta] = TMath::ACos((t1->Pz()+t2->Pz())/
697 TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) +
698 (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) +
699 (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz())));
700 // Flow variables
701 // cos(n*phi), sin(n*phi)
702 for(Int_t ih=0; ih<6; ++ih) {
703 if(gUsedVars[kPairCosNPhi+ih]) values[kPairCosNPhi+ih] = TMath::Cos(values[kPairPhi]*(1.0+ih));
704 if(gUsedVars[kPairSinNPhi+ih]) values[kPairSinNPhi+ih] = TMath::Sin(values[kPairPhi]*(1.0+ih));
705 }
706 // VZERO
707 for(Int_t iVZEROside=0; iVZEROside<3; ++iVZEROside) {
708 // v3 using VZERO Psi_2
709 if(gUsedVars[kPairVZEROFlowV3Psi2+iVZEROside])
710 values[kPairVZEROFlowV3Psi2+iVZEROside] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kVZERORP+iVZEROside*6+1]));
711 for(Int_t ih=0; ih<6; ++ih) {
712 // vn using VZERO Psi_n
713 if(gUsedVars[kPairVZEROFlowVn+iVZEROside*6+ih])
714 values[kPairVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih])*(ih+1));
715 // phi - Psi_n
716 if(gUsedVars[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih])
717 values[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih]);
718 }
719 }
720 // TPC
721 // v3 using Psi_2
722 if(gUsedVars[kPairTPCFlowV3Psi2])
723 values[kPairTPCFlowV3Psi2] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kTPCRP+1]));
724 for(Int_t ih=0; ih<6; ++ih) {
725 // vn using Psi_n
726 if(gUsedVars[kPairTPCFlowVn+ih])
727 values[kPairTPCFlowVn+ih] = TMath::Cos(DeltaPhi(values[kPairPhi],values[kTPCRP+ih])*(ih+1));
728 if(gUsedVars[kPairDeltaPhiTPCFlowVn+ih])
729 values[kPairDeltaPhiTPCFlowVn+ih] = DeltaPhi(values[kPairPhi], values[kTPCRP+ih]);
730 }
731}
732
733
734//__________________________________________________________________
735void DstCommonMacros::FillCorrelationInfo(AliCorrelationReducedPair* p, AliCorrelationReducedTrack* t, Float_t* values) {
736 //
737 // fill pair-track correlation information
738 //
739 if(gUsedVars[kDeltaPhi])
740 values[kDeltaPhi] = DeltaPhi(p->Phi(), t->Phi());
741
742 if(gUsedVars[kDeltaTheta])
743 values[kDeltaTheta] = p->Theta() - t->Theta();
744
745 if(gUsedVars[kDeltaEta]) values[kDeltaEta] = p->Eta() - t->Eta();
746}
747
748
749//__________________________________________________________________
750void DstCommonMacros::DoEventMixing(TList* list1, TList* list2, Float_t* values, Int_t type, const Char_t* histClass) {
751 //
752 // Do the event mixing
753 //
754 if(type==0) return;
755
756 cout << "mixing ..." << endl;
757 switch(type) {
758 case -1:
759 EventMixingPairTracks(list1, list2, values); // list1 contains pairs; list2 contains tracks
760 break;
761 case AliCorrelationReducedPair::kJpsiToEE:
762 // type is needed to encode the type of resonance which is needed for the mass assumption of legs when calculating the invariant mass
763 EventMixingResonanceLegs(list1, list2, values, type, histClass); // list1 contains positive legs; list2 contains negative legs
764 break;
765 default:
766 break;
767 }
768
769 for(Int_t i=list1->GetEntries()-1; i>=0; --i) ((TList*)list1->At(i))->RemoveAll();
770 for(Int_t j=list2->GetEntries()-1; j>=0; --j) ((TList*)list2->At(j))->RemoveAll();
771 list1->Clear();
772 list2->Clear();
773}
774
775
776//__________________________________________________________________
777void DstCommonMacros::EventMixingPairTracks(TList* pairLists, TList* trackLists, Float_t* values) {
778 //
779 // Make event mixing for pair - track correlations
780 //
781 cout << "Mixing category " << gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data() << endl;
782 Int_t entries = pairLists->GetEntries(); // we should have the same number of entries in both pair and track lists
783 if(entries<2) return;
784
785 TList* pairs = 0x0;
786 TList* tracks = 0x0;
787 AliCorrelationReducedPair* pair = 0x0;
788 AliCorrelationReducedTrack* track = 0x0;
789
790 TIter nextPairList(pairLists);
791 for(Int_t ie1=0; ie1<entries; ++ie1) {
792 pairs = (TList*)nextPairList();
793
794 TIter nextTrackList(trackLists);
795 for(Int_t ie2=0; ie2<entries; ++ie2) {
796 tracks = (TList*)nextTrackList();
797 if(ie1==ie2) continue;
798
799 TIter nextPair(pairs);
800 for(Int_t ip=0; ip<pairs->GetEntries(); ++ip) {
801 pair = (AliCorrelationReducedPair*)nextPair();
802 FillPairInfo(pair, values);
803
804 TIter nextTrack(tracks);
805 Int_t leadingIdx = -1;
806 Float_t leadingPt = 0.0;
807 for(Int_t it=0; it<tracks->GetEntries(); ++it) {
808 track = (AliCorrelationReducedTrack*)nextTrack();
809 FillTrackInfo(track, values);
810
811 FillCorrelationInfo(pair, track, values);
812 if(pair->PairType()==1)
813 FillHistClass(Form("Correlation_ME_US%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
814 if(pair->PairType()==0) {
815 FillHistClass(Form("Correlation_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values);
816 FillHistClass(Form("Correlation_ME_LSpp%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values);
817 }
818 if(pair->PairType()==2) {
819 FillHistClass(Form("Correlation_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values);
820 FillHistClass(Form("Correlation_ME_LSmm%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values);
821 }
822 if(track->Pt()>leadingPt) {
823 leadingPt = track->Pt();
824 leadingIdx = it;
825 }
826 } // end loop over tracks
827 if(leadingIdx!=-1) {
828 track = (AliCorrelationReducedTrack*)tracks->At(leadingIdx);
829 FillTrackInfo(track, values);
830 FillCorrelationInfo(pair, track, values);
831 if(pair->PairType()==1)
832 FillHistClass(Form("Correlation_LeadingPt_ME_US%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
833 if(pair->PairType()==0) {
834 FillHistClass(Form("Correlation_LeadingPt_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
835 FillHistClass(Form("Correlation_LeadingPt_ME_LSpp%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
836 }
837 if(pair->PairType()==2) {
838 FillHistClass(Form("Correlation_LeadingPt_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
839 FillHistClass(Form("Correlation_LeadingPt_ME_LSmm%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values);
840 }
841 }
842 } // end loop over pairs
843 } // end loop over second event
844 } // end loop over first event
845}
846
847
848//__________________________________________________________________
849void DstCommonMacros::EventMixingResonanceLegs(TList* posLists, TList* negLists, Float_t* values, Int_t type, const Char_t* histClass) {
850 //
851 // Do event mixing with the legs of a resonance to obtain the background
852 //
853 cout << "Mixing event class " << histClass << ", (cent/vtx): " << values[kCentVZERO] << "/" << values[kVtxZ] << endl;
854 Int_t entries = posLists->GetEntries();
855 cout << "mixing positives: " << entries << endl;
856 cout << "mixing negatives: " << negLists->GetEntries() << endl;
857 if(entries<2) return;
858
859 TList* positives1 = 0x0; // list of tracks in the first event
860 //TList* negatives1 = 0x0;
861 //TList* positives2 = 0x0; // list of tracks in the second event
862 TList* negatives2 = 0x0;
863 AliCorrelationReducedTrack* posTrack1 = 0x0;
864 //AliCorrelationReducedTrack* negTrack1 = 0x0;
865 //AliCorrelationReducedTrack* posTrack2 = 0x0;
866 AliCorrelationReducedTrack* negTrack2 = 0x0;
867
868 TIter nextPosList1(posLists);
869 TIter nextNegList1(negLists);
870 for(Int_t ie1=0; ie1<entries; ++ie1) {
871 positives1 = (TList*)nextPosList1();
872 //negatives1 = (TList*)nextNegList1();
873
874 TIter nextPosList2(posLists);
875 TIter nextNegList2(negLists);
876 for(Int_t ie2=0; ie2<entries; ++ie2) {
877 //positives2 = (TList*)nextPosList2();
878 negatives2 = (TList*)nextNegList2();
879 if(ie1==ie2) continue; // no SE mixing
880
881 TIter nextPosLeg1(positives1);
882 while((posTrack1=(AliCorrelationReducedTrack*)nextPosLeg1())) {
883 /*TIter nextPosLeg2(positives2);
884 while((posTrack2=(AliCorrelationReducedTrack*)nextPosLeg2())) {
885 FillPairInfo(posTrack1, posTrack2, type, values); // ++ combinations
886 FillHistClass(Form("%s_LSpp", histClass),values);
887 } // end while over pos legs in event 2
888 */
889 TIter nextNegLeg2(negatives2);
890 while((negTrack2=(AliCorrelationReducedTrack*)nextNegLeg2())) {
891 FillPairInfo(posTrack1, negTrack2, type, values); // +- combinations
892 if(IsPairSelectedEM(values))
893 FillHistClass(histClass, values);
894 } // end while over neg legs in event 2
895 } // end while over pos legs in event 1
896 /*TIter nextNegLeg1(negatives1);
897 while((negTrack1=(AliCorrelationReducedTrack*)nextNegLeg1())) {
898 TIter nextNegLeg2(negatives2);
899 while((negTrack2=(AliCorrelationReducedTrack*)nextNegLeg2())) {
900 FillPairInfo(negTrack1, negTrack2, type, values); // -- combinations
901 FillHistClass(Form("%s_LSmm", histClass),values);
902 } // end while over neg legs in event 2
903 } // end while over neg legs in event 1
904 */
905 } // end for over event 2
906 } // end for over event 1
907}
908
909
910
911
912//________________________________________________________________
913Double_t DstCommonMacros::DeltaPhi(Double_t phi1, Double_t phi2) {
914 //
915 // compute the delta of two angles defined in the (-pi,+pi) interval
916 //
917 Double_t delta = phi1-phi2;
918 if(delta>2.0*TMath::Pi()) delta -= 2.0*TMath::Pi();
919 if(delta<0.0) delta += 2.0*TMath::Pi();
920 /*Double_t delta = phi2;
921 if(phi2<0.0) delta += 2.0*TMath::Pi();
922 delta = phi1-delta;
923 if(delta>TMath::Pi()) delta = delta - 2.0*TMath::Pi();
924 if(delta<-1.*TMath::Pi()) delta = 2.0*TMath::Pi() + delta;
925 */
926 return delta;
927}
928
929
930//__________________________________________________________________
931void DstCommonMacros::FillHistClass(const Char_t* className, Float_t* values) {
932 //
933 // fill a class of histograms
934 //
935 THashList* hList = (THashList*)gHistLists->FindObject(className);
936 if(!hList) {
937 //cout << "Warning in FillHistClass(): Histogram class " << className << " not found" << endl;
938 return;
939 }
940
941 TIter next(hList);
942 TObject* h=0x0;
943 while((h=next())) {
944 UInt_t id = h->GetUniqueID();
945 UInt_t varX = id%kNVars;
946 UInt_t varY = (id/kNVars)%kNVars;
947 UInt_t varZ = (id/kNVars/kNVars)%kNVars;
948 if(id<kNVars) {
949 ((TH1F*)h)->Fill(values[varX]);
950 } else if(id<kNVars*kNVars) {
951 if(((TH1*)h)->GetDimension()==1) ((TProfile*)h)->Fill(values[varX],values[varY]);
952 if(((TH1*)h)->GetDimension()==2) ((TH2F*)h)->Fill(values[varX],values[varY]);
953 } else {
954 if(((TH1*)h)->GetDimension()==2) ((TProfile2D*)h)->Fill(values[varX],values[varY],values[varZ]);
955 if(((TH1*)h)->GetDimension()==3) ((TH3F*)h)->Fill(values[varX],values[varY],values[varZ]);
956 };
957 }
958}
959
960
961//__________________________________________________________________
962void DstCommonMacros::AddHistClass(const Char_t* histClass) {
963 //
964 // Add a new histogram list
965 //
966 if(!gHistLists)
967 gHistLists = new TObjArray();
968 gHistLists->SetOwner();
969 gHistLists->SetName("histos");
970
971 if(gHistLists->FindObject(histClass)) {
972 cout << "Warning in AddHistClass: Cannot add histogram class " << histClass
973 << " because it already exists." << endl;
974 return;
975 }
976 THashList* hList=new THashList;
977 hList->SetOwner(kTRUE);
978 hList->SetName(histClass);
979 gHistLists->Add(hList);
980}
981
982
983//_________________________________________________________________
984Int_t DstCommonMacros::ValidateHistogramName(THashList* hList, const Char_t* name) {
985 //
986 // check whether a histogram with this name already exist, and how many of them
987 //
988 for(Int_t i=0; i<hList->GetEntries(); ++i) {
989 TString hname = hList->At(i)->GetName();
990 TObjArray* arr=hname.Tokenize("#");
991 TString nameRoot = arr->At(0)->GetName();
992 if(nameRoot.CompareTo(name)==0) {
993 cout << "Warning in AddHistogram(): Histogram " << name << " already exists in this list (" << nameRoot.Data() << ")" << endl;
994 return -1;
995 }
996 }
997 Int_t nnames = 0;
998 for(Int_t il=0; il<gHistLists->GetEntries(); ++il) {
999 THashList* tlist = (THashList*)gHistLists->At(il);
1000 for(Int_t ih=0; ih<tlist->GetEntries(); ++ih) {
1001 TString hname = tlist->At(ih)->GetName();
1002 TObjArray* arr=hname.Tokenize("#");
1003 TString nameRoot = arr->At(0)->GetName();
1004 if(nameRoot.CompareTo(name)==0) ++nnames;
1005 }
1006 }
1007 return nnames;
1008}
1009
1010
1011//_________________________________________________________________
1012void DstCommonMacros::AddHistogram(const Char_t* histClass,
1013 const Char_t* name, const Char_t* title, Bool_t isProfile,
1014 Int_t nXbins, Double_t xmin, Double_t xmax, Int_t varX,
1015 Int_t nYbins, Double_t ymin, Double_t ymax, Int_t varY,
1016 Int_t nZbins, Double_t zmin, Double_t zmax, Int_t varZ,
1017 const Char_t* xLabels, const Char_t* yLabels, const Char_t* zLabels) {
1018 //
1019 // add a histogram
1020 //
1021 THashList* hList = (THashList*)gHistLists->FindObject(histClass);
1022 if(hList->FindObject(name)) {
1023 cout << "Warning in AddHistogram(): Histogram " << name << " already exists" << endl;
1024 return;
1025 }
1026 /*Int_t nnames = ValidateHistogramName(hList, name);
1027 if(nnames<0) return;
1028 TString hname = Form("%s#v%d", name, nnames);
1029 cout << "Assigned name " << hname.Data() << endl;*/
1030 TString hname = name;
1031
1032 Int_t dimension = 1;
1033 if(varY!=kNothing) dimension = 2;
1034 if(varZ!=kNothing) dimension = 3;
1035
1036 TString titleStr(title);
1037 TObjArray* arr=titleStr.Tokenize(";");
1038
1039 TH1* h=0x0;
1040 switch(dimension) {
1041 case 1:
1042 h=new TH1F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax);
1043 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX) << endl;
1044 h->SetUniqueID(UInt_t(varX));
1045 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1046 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1047 gUsedVars[varX] = kTRUE;
1048 hList->Add(h);
1049 break;
1050 case 2:
1051 if(isProfile)
1052 h=new TProfile(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax);
1053 else
1054 h=new TH2F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax);
1055 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY) << endl;
1056 h->SetUniqueID(UInt_t(varX+kNVars*varY));
1057 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1058 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1059 if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName());
1060 if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels);
1061 gUsedVars[varX] = kTRUE;
1062 gUsedVars[varY] = kTRUE;
1063 hList->Add(h);
1064 break;
1065 case 3:
1066 if(isProfile)
1067 h=new TProfile2D(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax);
1068 else
1069 h=new TH3F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax,nZbins,zmin,zmax);
1070 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ) << endl;
1071 h->SetUniqueID(UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ));
1072 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1073 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1074 if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName());
1075 if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels);
1076 if(arr->At(3)) h->GetZaxis()->SetTitle(arr->At(3)->GetName());
1077 if(zLabels[0]!='\0') MakeAxisLabels(h->GetZaxis(), zLabels);
1078 gUsedVars[varX] = kTRUE;
1079 gUsedVars[varY] = kTRUE;
1080 gUsedVars[varZ] = kTRUE;
1081 hList->Add(h);
1082 break;
1083 }
1084}
1085
1086//_________________________________________________________________
1087void DstCommonMacros::AddHistogram(const Char_t* histClass,
1088 const Char_t* name, const Char_t* title, Bool_t isProfile,
1089 Int_t nXbins, Double_t* xbins, Int_t varX,
1090 Int_t nYbins, Double_t* ybins, Int_t varY,
1091 Int_t nZbins, Double_t* zbins, Int_t varZ,
1092 const Char_t* xLabels, const Char_t* yLabels, const Char_t* zLabels) {
1093 //
1094 // add a histogram
1095 //
1096 THashList* hList = (THashList*)gHistLists->FindObject(histClass);
1097 if(hList->FindObject(name)) {
1098 cout << "Warning in AddHistogram(): Histogram " << name << " already exists" << endl;
1099 return;
1100 }
1101 //Int_t nnames = ValidateHistogramName(hList, name);
1102 //if(nnames<0) return;
1103 //TString hname = Form("%s#v%d", name, nnames);
1104 TString hname = name;
1105
1106 Int_t dimension = 1;
1107 if(varY!=kNothing) dimension = 2;
1108 if(varZ!=kNothing) dimension = 3;
1109
1110 TString titleStr(title);
1111 TObjArray* arr=titleStr.Tokenize(";");
1112
1113 TH1* h=0x0;
1114 switch(dimension) {
1115 case 1:
1116 h=new TH1F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins);
1117 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX) << endl;
1118 h->SetUniqueID(UInt_t(varX));
1119 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1120 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1121 gUsedVars[varX] = kTRUE;
1122 hList->Add(h);
1123 break;
1124 case 2:
1125 if(isProfile)
1126 h=new TProfile(hname.Data(),arr->At(0)->GetName(),nXbins,xbins);
1127 else
1128 h=new TH2F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins);
1129 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY) << endl;
1130 h->SetUniqueID(UInt_t(varX+kNVars*varY));
1131 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1132 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1133 if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName());
1134 if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels);
1135 gUsedVars[varX] = kTRUE;
1136 gUsedVars[varY] = kTRUE;
1137 hList->Add(h);
1138 break;
1139 case 3:
1140 if(isProfile)
1141 h=new TProfile2D(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins);
1142 else
1143 h=new TH3F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins,nZbins,zbins);
1144 //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ) << endl;
1145 h->SetUniqueID(UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ));
1146 if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName());
1147 if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels);
1148 if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName());
1149 if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels);
1150 if(arr->At(3)) h->GetZaxis()->SetTitle(arr->At(3)->GetName());
1151 if(zLabels[0]!='\0') MakeAxisLabels(h->GetZaxis(), zLabels);
1152 gUsedVars[varX] = kTRUE;
1153 gUsedVars[varY] = kTRUE;
1154 gUsedVars[varZ] = kTRUE;
1155 hList->Add(h);
1156 break;
1157 }
1158}
1159
1160
1161//____________________________________________________________________________________
1162void DstCommonMacros::MakeAxisLabels(TAxis* ax, const Char_t* labels) {
1163 //
1164 // add bin labels to an axis
1165 //
1166 TString labelsStr(labels);
1167 TObjArray* arr=labelsStr.Tokenize(";");
1168 for(Int_t ib=1; ib<=ax->GetNbins(); ++ib) {
1169 if(ib>=arr->GetEntries()+1) break;
1170 ax->SetBinLabel(ib, arr->At(ib-1)->GetName());
1171 }
1172}
1173
1174
1175//____________________________________________________________________________________
1176void DstCommonMacros::InitFile(const Char_t* filename) {
1177 //
1178 // Open an existing ROOT file containing lists of histograms and initialize the global list pointer
1179 //
1180 //if(gHistFile) gHistFile->Close();
1181 TString histfilename="";
1182 if(gHistFile) histfilename = gHistFile->GetName();
1183 if(!histfilename.Contains(filename))
1184 gHistFile = new TFile(filename); // open file only if not already open
1185
1186 if(!gHistFile) {
1187 cout << "GlobalMacros.C::GetHistogram() : File " << filename << " not opened!!" << endl;
1188 return;
1189 }
1190 if(gHistFile->IsZombie()) {
1191 cout << "GlobalMacros.C::GetHistogram() : File " << filename << " not opened!!" << endl;
1192 return;
1193 }
1194 TList* list1 = gHistFile->GetListOfKeys();
1195 TKey* key1 = (TKey*)list1->At(0);
1196 gHistListsOld = (TDirectoryFile*)key1->ReadObj();
1197}
1198
1199
1200//____________________________________________________________________________________
1201void DstCommonMacros::CloseFile() {
1202 //
1203 // Close the opened file
1204 //
1205 gHistListsOld = 0x0;
1206 if(gHistFile && gHistFile->IsOpen()) gHistFile->Close();
1207}
1208
1209
1210//____________________________________________________________________________________
1211TObject* DstCommonMacros::GetHistogram(const Char_t* listname, const Char_t* hname) {
1212 //
1213 // Retrieve a histogram from the list hlist
1214 //
1215 if(!gHistListsOld) {
1216 cout << "GlobalMacros.C::GetHistogram() : The main list was not initialized." << endl;
1217 cout << " A ROOT file must pe initialized first!!" << endl;
1218 return 0x0;
1219 }
1220 TKey* listKey = gHistListsOld->FindKey(listname);
1221 TDirectoryFile* hlist = (TDirectoryFile*)listKey->ReadObj();
1222 TKey* key = hlist->FindKey(hname);
1223 return key->ReadObj();
1224}
1225
1226
1227//_________________________________________________________________
1228TChain* DstCommonMacros::GetChain(const Char_t* filename, Int_t howMany, Int_t offset, Long64_t& entries,
1229 TChain* friendChain/*=0x0*/, const Char_t* friendChainFile/*=0x0*/) {
1230 //
1231 // read an ascii file containing a list of root files with reduced trees
1232 // and build a TChain
1233 //
1234 cout << "Creating the data chain from " << filename << " ..." << flush;
1235 TChain* chain = new TChain("DstTree");
1236 ifstream inBuf;
1237 inBuf.open(filename);
1238 Int_t index = 0;
1239 while(inBuf.good()) {
1240 Char_t str[512];
1241 inBuf.getline(str,512,'\n');
1242
1243 if(index<offset) {++index; continue;}
1244 if(index>=offset+howMany) break;
1245
1246 TString strstr = str;
1247 if(!strstr.Contains(".root")) continue;
1248 cout << endl << "Adding file " << str << endl;
1249 chain->Add(str);
1250 if(friendChain) {
1251 TObjArray* arr = strstr.Tokenize("/");
1252 strstr.ReplaceAll(arr->At(arr->GetEntries()-1)->GetName(),friendChainFile);
1253 friendChain->Add(strstr.Data());
1254 }
1255 ++index;
1256 }
1257 inBuf.close();
1258 entries = chain->GetEntries();
1259 Long64_t entriesFriend = (friendChain ? friendChain->GetEntries() : 0);
1260 cout << "DstCommonMacros::GetChain() Chain entries = " << entries << endl;
1261 if(friendChain)
1262 cout << "DstCommonMacros::GetChain() Friend chain entries = " << entriesFriend << endl;
1263 if(friendChain && (entries!=entriesFriend)) {
1264 cout << "DstCommonMacros::GetChain() The friend chain does not have the same number of entries as the main chain!!!" << endl;
1265 cout << " Check it out and retry !!" << endl;
1266 return 0x0;
1267 }
1268 cout << " done" << endl;
1269 return chain;
1270}
1271
1272
1273//__________________________________________________________________
1274void DstCommonMacros::WriteOutput(TFile* save) {
1275 //
1276 // Write the histogram lists in the output file
1277 //
1278 cout << "Writing the output to " << save->GetName() << " ... " << flush;
1279 //TFile* save=new TFile(filename,"RECREATE");
1280 TDirectory* mainDir = save->mkdir(gHistLists->GetName());
1281 mainDir->cd();
1282 for(Int_t i=0; i<gHistLists->GetEntries(); ++i) {
1283 THashList* list = (THashList*)gHistLists->At(i);
1284 TDirectory* dir = mainDir->mkdir(list->GetName());
1285 dir->cd();
1286 list->Write();
1287 mainDir->cd();
1288 }
1289 save->Close();
1290 cout << "done" << endl;
1291}
1292
1293
1294//__________________________________________________________________
1295Bool_t DstCommonMacros::IsPairSelectedEM(Float_t* values) {
1296 //
1297 // pair selection on the pair resulting from the event mixing
1298 //
1299 if(values[kPairPt]<gkJpsiPtCut[0]) return kFALSE;
1300 if(values[kPairPt]>gkJpsiPtCut[1]) return kFALSE;
1301 return kTRUE;
1302}