]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
K0s analysis update (Matt Steinpries)
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / K0Analysis / AliFemtoK0Analysis.cxx
CommitLineData
41dfc4d3
DRG		 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 is 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// This class is used to perform femtoscopic analysis on K0s particles,
19// which are reconstructed using the AliAODv0 class.
20//
21// authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)
22//
23// Change log:
24// - TOF mismatch function calls changed (4/18/13)
25// - added minimum decay length cut (rarely used though) (3/28/13)
26// - K0 multiplicity histogram now filled with "unskippedCount" instead
27// of k0Count (which included skipped k0s with shared daughters)
28// (3/25/13)
29// - added hists for 3D mom. in LF and PRF (3/28/13)
30// - changed calling of PIDResponse (should be same actions) (3/28/13)
31// - keep "side" K0s for mass plot (4/18)
32// - tweaked loading and skipping appropriately
33// - use merit test to skip sides (against good and side K0s)
34// - a good K0 cant be skipped by a side
35// - moved TPC propagation (via Hans' method) up to v0 level, which now
36// uses an AliAODTrack(AliVTrack) instead of AliESDtrack (5/31/13)
37// - added primary vertex subtraction in TPC propagation (5/31/13)
38// - removed all instances of AliESDtrack usage (5/31/13)
39// - removed old separation method/histograms (5/31/13)
40// - tidied up LCMS boost (6/10/13)
41// - added new boosting prescription, get q out-side-long for LCMS and PRF (6/24/13)
42// - added histograms and values for LCMS momenta (for simulation)
43// - added random particle order switch in correlations (9/09/13)
44// - added more bins for 3D OSL analysis (9/19/13)
45// - added merit cut choice, pass as argument (10/16/13)
46// - 1-mass, 2-v0dca, 3-dddca, 4-combination (used to be v0dca)
47// - added passable argument for two-track minimum separation (10/16/13)
48// - added boolean to turn off field-sign dependence for train (10/30/13)
49// - changed destructors to minimize lost memory (11/27/13)
50// - added Case3D to switch off all 3D objects (11/27/13)
51// - added centrality flattening routine (and switch) (12/04/13)
52// - added event plane stuff (12/11/13)
53////////////////////////////////////////////////////////////////////////////////
54
55
56
57#include <iostream>
58#include <math.h>
59#include "TMath.h"
60#include "TChain.h"
61#include "TFile.h"
62#include "TKey.h"
63#include "TObject.h"
64#include "TObjString.h"
65#include "TList.h"
66#include "TTree.h"
67#include "TH1F.h"
68#include "TH1D.h"
69#include "TH2D.h"
70#include "TH3D.h"
71#include "TProfile.h"
72#include "TCanvas.h"
73#include "TRandom3.h"
74
75#include "AliAnalysisTask.h"
76#include "AliAnalysisManager.h"
77
78#include "AliAODEvent.h"
79#include "AliAODInputHandler.h"
80#include "AliAODMCParticle.h"
81#include "AliAODv0.h"
82#include "AliAODRecoDecay.h"
83#include "AliCentrality.h"
84
85#include "AliFemtoK0Analysis.h"
86
87#define PI 3.1415927
88
89
90// Author: Matt Steinpreis, adapted from Dhevan Gangadharan
91
92ClassImp(AliFemtoK0Analysis)
93
94//________________________________________________________________________
95AliFemtoK0Analysis::AliFemtoK0Analysis():
96AliAnalysisTaskSE(),
97 fSignDep(kFALSE),
98 fFieldPos(kTRUE),
99 fOnlineCase(kTRUE),
100 fMeritCase(kTRUE),
101 fCase3D(kFALSE),
102 fMinDecayLength(0.0),
103 fMeritCutChoice(0),
104 fMinSep(0.0),
105 fFlatCent(kFALSE),
106 fPsiBinning(kFALSE),
107 fEventCount(0),
108 fEC(0x0),
109 fEvt(0X0),
110 fRandomNumber(0x0),
111 fName(0x0),
112 fAOD(0x0),
113 fOutputList(0x0),
114 fPidAOD(0x0)
115{
116}
117//________________________________________________________________________
118AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool SignDep, bool FieldPositive, bool OnlineCase, bool MeritCase, bool Case3D, float MinDL, int MeritCutChoice, float MinSep, bool FlatCent, bool PsiBinning)
119: AliAnalysisTaskSE(name),
120 fSignDep(SignDep),
121 fFieldPos(FieldPositive),
122 fOnlineCase(OnlineCase),
123 fMeritCase(MeritCase),
124 fCase3D(Case3D),
125 fMinDecayLength(MinDL),
126 fMeritCutChoice(MeritCutChoice),
127 fMinSep(MinSep),
128 fFlatCent(FlatCent),
129 fPsiBinning(PsiBinning),
130 fEventCount(0),
131 fEC(0x0),
132 fEvt(0X0),
133 fRandomNumber(0x0),
134 fName(name),
135 fAOD(0x0),
136 fOutputList(0x0),
137 fPidAOD(0x0)
138{
139 //main constructor
140 fSignDep = SignDep;
141 fFieldPos = FieldPositive;
142 fOnlineCase = OnlineCase;
143 fMeritCase = MeritCase;
144 fCase3D = Case3D;
145 fMinDecayLength = MinDL;
146 fMeritCutChoice = MeritCutChoice;
147 fMinSep = MinSep;
148 fFlatCent = FlatCent;
149 fPsiBinning = PsiBinning;
150
151 // Define output slots here
152 // Output slot #1
153 DefineOutput(1, TList::Class());
154
155}
156//________________________________________________________________________
157AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
158: AliAnalysisTaskSE(obj.fName),
159 fSignDep(obj.fSignDep),
160 fFieldPos(obj.fFieldPos),
161 fOnlineCase(obj.fOnlineCase),
162 fMeritCase(obj.fMeritCase),
163 fCase3D(obj.fCase3D),
164 fMinDecayLength(obj.fMinDecayLength),
165 fMeritCutChoice(obj.fMeritCutChoice),
166 fMinSep(obj.fMinSep),
167 fFlatCent(obj.fFlatCent),
168 fPsiBinning(obj.fPsiBinning),
169 fEventCount(obj.fEventCount),
170 fEC(obj.fEC),
171 fEvt(obj.fEvt),
172 fRandomNumber(obj.fRandomNumber),
173 fName(obj.fName),
174 fAOD(obj.fAOD),
175 fOutputList(obj.fOutputList),
176 fPidAOD(obj.fPidAOD)
177{
178}
179//________________________________________________________________________
180AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
181{
182 //Assignment operator
183 if (this == &obj) return *this;
184
185 fSignDep = obj.fSignDep;
186 fFieldPos = obj.fFieldPos;
187 fOnlineCase = obj.fOnlineCase;
188 fMeritCase = obj.fMeritCase;
189 fCase3D = obj.fCase3D;
190 fMinDecayLength= obj.fMinDecayLength;
191 fMeritCutChoice= obj.fMeritCutChoice;
192 fMinSep = obj.fMinSep;
193 fFlatCent = obj.fFlatCent;
194 fPsiBinning = obj.fPsiBinning;
195 fEventCount = obj.fEventCount;
196 fEC = obj.fEC;
197 fEvt = obj.fEvt;
198 fRandomNumber = obj.fRandomNumber;
199 fName = obj.fName;
200 fAOD = obj.fAOD;
201 fOutputList = obj.fOutputList;
202 fPidAOD = obj.fPidAOD;
203
204 return (*this);
205}
206//________________________________________________________________________
207AliFemtoK0Analysis::~AliFemtoK0Analysis()
208{
209 // Destructor
210 for(unsigned short i=0; i<kZVertexBins; i++)
211 {
212 for(unsigned short j=0; j<kCentBins; j++)
213 {
214 for(unsigned short k=0; k<kPsiBins; k++)
215 {
216 fEC[i][j][k]->~AliFemtoK0EventCollection();
217 fEC[i][j][k] = NULL;
218 }
219 delete [] fEC[i][j]; fEC[i][j] = NULL;
220 }
221 delete[] fEC[i]; fEC[i] = NULL;
222 }
223 delete[] fEC; fEC = NULL;
224
225 if(fEC){ delete fEC; fEC = NULL;}
226 if(fRandomNumber){ delete fRandomNumber; fRandomNumber = NULL;}
227 if(fAOD){ delete fAOD; fAOD = NULL;}
228 if(fOutputList){ delete fOutputList; fOutputList = NULL;}
229 if(fPidAOD){ delete fPidAOD; fPidAOD = NULL;}
230}
231//________________________________________________________________________
232void AliFemtoK0Analysis::MyInit()
233{
234
235 // One can set global variables here
236 fEventCount = 0;
237
238 fEC = new AliFemtoK0EventCollection ***[kZVertexBins];
239 for(unsigned short i=0; i<kZVertexBins; i++)
240 {
241 fEC[i] = new AliFemtoK0EventCollection **[kCentBins];
242
243 for(unsigned short j=0; j<kCentBins; j++)
244 {
245 fEC[i][j] = new AliFemtoK0EventCollection *[kPsiBins];
246
247 for(unsigned short k=0; k<kPsiBins; k++)
248 {
249 fEC[i][j][k] = new AliFemtoK0EventCollection(kEventsToMix+1, kMultLimit);
250 }
251 }
252 }
253
254 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
255 fPidAOD = aodH->GetPIDResponse();
256
257 fRandomNumber = new TRandom3(); //for 3D, random sign switching
258 fRandomNumber->SetSeed(0);
259
260}
261//________________________________________________________________________
262void AliFemtoK0Analysis::UserCreateOutputObjects()
263{
264 // Create histograms
265 // Called once
266
267 MyInit();// Initialize my settings
268
269 fOutputList = new TList();
270 fOutputList->SetOwner();
271
272 TH1F *fHistCent = new TH1F("fHistCent","",100,0,100);
273 fOutputList->Add(fHistCent);
274 TH1F *fHistCentFlat = new TH1F("fHistCentFlat","",100,0,100);
275 fOutputList->Add(fHistCentFlat);
276 TH1F *fHistCentUsed = new TH1F("fHistCentUsed","",100,0,100);
277 fOutputList->Add(fHistCentUsed);
278
279 //pion parameters
280 TH1F* fHistDCAPiPlus = new TH1F("fHistDCAPiPlus","",100,0,10);
281 fOutputList->Add(fHistDCAPiPlus);
282 TH1F* fHistDCAPiMinus = new TH1F("fHistDCAPiMinus","",100,0,10);
283 fOutputList->Add(fHistDCAPiMinus);
284 TH1F* fHistDCADaughters = new TH1F("fHistDCADaughters", "DCA of pions to each other", 50, 0., 0.5);
285 fOutputList->Add(fHistDCADaughters);
286 TH2F* fHistK0PiPlusPt = new TH2F("fHistK0PiPlusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
287 fOutputList->Add(fHistK0PiPlusPt);
288 TH2F* fHistK0PiMinusPt = new TH2F("fHistK0PiMinusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
289 fOutputList->Add(fHistK0PiMinusPt);
290 TH1F* fHistDaughterPhi = new TH1F("fHistDaughterPhi","",180,-PI,PI);
291 fOutputList->Add(fHistDaughterPhi);
292
293 //K0 parameters
294 TH1F* fHistMultK0 = new TH1F("fHistMultK0", "K0 multiplicity", 51, -0.5, 51-0.5);
295 fOutputList->Add(fHistMultK0);
296 TH2F* fHistPtK0 = new TH2F("fHistPtK0", "K0 pt distribution",kCentBins,.5,kCentBins+.5, 100, 0., 10.);
297 fOutputList->Add(fHistPtK0);
298 TH1F* fHistDecayLengthK0 = new TH1F("fHistDecayLengthK0", "K0 decay length", 100, 0., 100.);
299 fOutputList->Add(fHistDecayLengthK0);
300 TH1F* fHistDCAK0 = new TH1F("fHistDCAK0", "DCA of K0 to primary vertex", 40, 0., 0.4);
301 fOutputList->Add(fHistDCAK0);
302 TH2F* fHistKtK0 = new TH2F("fHistKtK0", "Kt distribution of K0 pairs", kCentBins, .5, kCentBins+.5, 300, 0., 3.);
303 fOutputList->Add(fHistKtK0);
304
305 TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);
306 TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);
307 TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);
308 TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);
309 TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);
310 TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);
311 TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);
312 fOutputList->Add(fHistPx);
313 fOutputList->Add(fHistPy);
314 fOutputList->Add(fHistPz);
315 fOutputList->Add(fHistPxCM);
316 fOutputList->Add(fHistPyCM);
317 fOutputList->Add(fHistPzCM);
318 fOutputList->Add(fHistKsCM);
319
320 TH1F* fHistPOutLCMS = new TH1F("fHistPOutLCMS","",200,0,2);
321 TH1F* fHistPSideLCMS = new TH1F("fHistPSideLCMS","",200,0,2);
322 TH1F* fHistPLongLCMS = new TH1F("fHistPLongLCMS","",200,0,2);
323 fOutputList->Add(fHistPOutLCMS);
324 fOutputList->Add(fHistPSideLCMS);
325 fOutputList->Add(fHistPLongLCMS);
326
327 //pair gamma (LCMS to PRF, OSL)
328 TH2F* fHistGamma = new TH2F("fHistGamma","Gamma from LCMS to PRF",500,1,5,100,0,1);
329 fOutputList->Add(fHistGamma);
330
331 //invariant mass distributions
332 TH3F* fHistMass = new TH3F("fHistMass","",kCentBins,.5,kCentBins+.5,50,0.,5.,400,.3,.7);
333 fOutputList->Add(fHistMass);
334 //TH3F* fHistMassPtCFK0 = new TH3F("fHistMassPtCFK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
335 //fOutputList->Add(fHistMassPtCFK0);
336 //TH3F* fHistMassPtCFBkgK0 = new TH3F("fHistMassPtCFBkgK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
337 //fOutputList->Add(fHistMassPtCFBkgK0);
338 //TH3F* fHistMassQKt = new TH3F("fHistMassQKt","",100,0,1,200,0,2,200,.4,.6);
339 //fOutputList->Add(fHistMassQKt);
340 //TH3F* fHistMassKtK0 = new TH3F("fHistMassKtK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
341 //fOutputList->Add(fHistMassKtK0);
342 //TH3F* fHistMassKtBkgK0 = new TH3F("fHistMassKtBkgK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
343 //fOutputList->Add(fHistMassKtBkgK0);
344
345 //separation studies
346 TH1F* fHistSepNumPos = new TH1F("fHistSepNumPos","",200,0,20);
347 fOutputList->Add(fHistSepNumPos);
348 TH1F* fHistSepDenPos = new TH1F("fHistSepDenPos","",200,0,20);
349 fOutputList->Add(fHistSepDenPos);
350 TH1F* fHistSepNumNeg = new TH1F("fHistSepNumNeg","",200,0,20);
351 fOutputList->Add(fHistSepNumNeg);
352 TH1F* fHistSepDenNeg = new TH1F("fHistSepDenNeg","",200,0,20);
353 fOutputList->Add(fHistSepDenNeg);
354
355 TH2F* fHistSepNumPos2 = new TH2F("fHistSepNumPos2","",100,0,20,100,0,20);
356 TH2F* fHistSepDenPos2 = new TH2F("fHistSepDenPos2","",100,0,20,100,0,20);
357 TH2F* fHistSepNumNeg2 = new TH2F("fHistSepNumNeg2","",100,0,20,100,0,20);
358 TH2F* fHistSepDenNeg2 = new TH2F("fHistSepDenNeg2","",100,0,20,100,0,20);
359 fOutputList->Add(fHistSepNumPos2);
360 fOutputList->Add(fHistSepDenPos2);
361 fOutputList->Add(fHistSepNumNeg2);
362 fOutputList->Add(fHistSepDenNeg2);
363
364 TH2F* fHistSepDPC = new TH2F("fHistSepDPC","",200,-1,1,50,0,10);
365 TH2F* fHistSepDPCBkg = new TH2F("fHistSepDPCBkg","",200,-1,1,50,0,10);
366 fOutputList->Add(fHistSepDPC);
367 fOutputList->Add(fHistSepDPCBkg);
368
369 TH1F *fHistPsi = new TH1F("fHistPsi","",90,-PI/2,PI/2);
370 fOutputList->Add(fHistPsi);
371 TH2F *fHistPhiPsi = new TH2F("fHistPhiPsi","",kCentBins,.5,kCentBins+.5,180,0,2*PI);
372 fOutputList->Add(fHistPhiPsi);
373
374/////////Pair Distributions///////////////////
375
376 //1D Q invariant
377 TH3F* fHistQinvSignal = new TH3F("fHistQinvSignal","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
378 fOutputList->Add(fHistQinvSignal);
379 TH3F* fHistQinvBkg = new TH3F("fHistQinvBkg","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
380 fOutputList->Add(fHistQinvBkg);
381
382 //event plane
383 TH3F* fHistQinvSignalEPIn = new TH3F("fHistQinvSignalEPIn","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
384 fOutputList->Add(fHistQinvSignalEPIn);
385 TH3F* fHistQinvBkgEPIn = new TH3F("fHistQinvBkgEPIn","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
386 fOutputList->Add(fHistQinvBkgEPIn);
387 TH3F* fHistQinvSignalEPOut = new TH3F("fHistQinvSignalEPOut","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
388 fOutputList->Add(fHistQinvSignalEPOut);
389 TH3F* fHistQinvBkgEPOut = new TH3F("fHistQinvBkgEPOut","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
390 fOutputList->Add(fHistQinvBkgEPOut);
391
392 //mass bins within peak
393 //TH3F* fHistCLCLSignal = new TH3F("fHistCLCLSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
394 //TH3F* fHistCLCLBkg = new TH3F("fHistCLCLBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
395 //TH3F* fHistCLCRSignal = new TH3F("fHistCLCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
396 //TH3F* fHistCLCRBkg = new TH3F("fHistCLCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
397 //TH3F* fHistCRCRSignal = new TH3F("fHistCRCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
398 //TH3F* fHistCRCRBkg = new TH3F("fHistCRCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
399 //fOutputList->Add(fHistCLCLSignal);
400 //fOutputList->Add(fHistCLCLBkg);
401 //fOutputList->Add(fHistCLCRSignal);
402 //fOutputList->Add(fHistCLCRBkg);
403 //fOutputList->Add(fHistCRCRSignal);
404 //fOutputList->Add(fHistCRCRBkg);
405
406 //3D out-side-long
407 TH3F *fHist3DOSLSignal[10][4];
408 TH3F *fHist3DOSLBkg[10][4];
409
410 if(fCase3D){
411 for(int i3D=0;i3D<10;i3D++){
412 for(int j3D=0;j3D<4;j3D++){
413 TString *histname = new TString("fHist3DOSL");
414 *histname += i3D;
415 *histname += j3D;
416 histname->Append("Signal");
417 fHist3DOSLSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
418 fOutputList->Add(fHist3DOSLSignal[i3D][j3D]);
419 histname->Replace(12,6,"Bkg");
420 fHist3DOSLBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
421 fOutputList->Add(fHist3DOSLBkg[i3D][j3D]);
422 }
423 }
424 }
425
426 //3D Spherical Harmonics
427 //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
428 //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
429 //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
430 //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
431 //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
432 //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
433 //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
434 //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
435 //fOutputList->Add(fHistSHCentLowKt);
436 //fOutputList->Add(fHistSHCentHighKt);
437 //fOutputList->Add(fHistSHSemiCentLowKt);
438 //fOutputList->Add(fHistSHSemiCentHighKt);
439 //fOutputList->Add(fHistSHCentLowKtBkg);
440 //fOutputList->Add(fHistSHCentHighKtBkg);
441 //fOutputList->Add(fHistSHSemiCentLowKtBkg);
442 //fOutputList->Add(fHistSHSemiCentHighKtBkg);
443
444 //side-side
445 //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
446 //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
447 //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
448 //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
449 //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
450 //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
451 //fOutputList->Add(fHistLeftLeftSignal);
452 //fOutputList->Add(fHistLeftRightSignal);
453 //fOutputList->Add(fHistRightRightSignal);
454 //fOutputList->Add(fHistLeftLeftBkg);
455 //fOutputList->Add(fHistLeftRightBkg);
456 //fOutputList->Add(fHistRightRightBkg);
457
458 //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
459 //fOutputList->Add(fHistSplitK0Sides);
460 //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
461 //fOutputList->Add(fHistSplitK0Centers);
462 //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
463 //fOutputList->Add(fHistQinvSignalNoSplit);
464
465 PostData(1, fOutputList);
466
467}
468
469//________________________________________________________________________
470void AliFemtoK0Analysis::Exec(Option_t *)
471{
472 // Main loop
473 // Called for each event
474 //cout<<"=========== Event # "<<fEventCount+1<<" ==========="<<endl;
475 fEventCount++;
476 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
477 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
478
479 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
480 bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
481 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
482 if(!isSelected) {
483 //cout << "Failed trigger selection." << endl;
484 return;
485 }
486
487 ///////////////////////////////////////////////////////////
488
489 unsigned int statusPos=0;
490 unsigned int statusNeg=0;
491
492 float bField=0;
493 bField = fAOD->GetMagneticField();
494 if(bField == 0) return;
495 if(fSignDep){
496 if(fFieldPos && bField < 0) return;
497 if(!fFieldPos && bField > 0) return;
498 }
499
500
41dfc4d3
DRG		 501 /////////////////////////////////////////////////
502
503
504 //Centrality selection
505
506 AliCentrality *centrality = fAOD->GetCentrality();
507 float percent = centrality->GetCentralityPercentile("V0M");
508 int centBin=0;
509 //Printf("Centrality percent = %f", percent);
510
511 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
512 float multV0A=aodV0->GetMTotV0A();
513 float multV0C=aodV0->GetMTotV0C();
514
515 if(percent < 0) {
516 //Printf("No centrality info");
517 return;
518 }
519 if(percent < 0.1 && (multV0A + multV0C < 19500)){
520 //Printf("No centrality info");
521 return;
522 }
523 else if(percent <= 5) centBin=15;
524 else if(percent <= 10) centBin=14;
525 else if(percent <= 15) centBin=13;
526 else if(percent <= 20) centBin=12;
527 else if(percent <= 25) centBin=11;
528 else if(percent <= 30) centBin=10;
529 else if(percent <= 35) centBin=9;
530 else if(percent <= 40) centBin=8;
531 else if(percent <= 45) centBin=7;
532 else if(percent <= 50) centBin=6;
533 else if(percent <= 55) centBin=5;
534 else if(percent <= 60) centBin=4;
535 else if(percent <= 65) centBin=3;
536 else if(percent <= 70) centBin=2;
537 else if(percent <= 75) centBin=1;
538 else if(percent <= 80) centBin=0;
539 else {
540 //Printf("Skipping Peripheral Event");
541 return;
542 }
543 if(percent > 10 && isCentral) return;
544 ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
545
546 //flatten centrality dist.
547 if(percent < 9){
548 if(fFlatCent){
549 if(RejectEventCentFlat(bField,percent)) return;
550 }
551 }
552 ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
553
554 //Vertexing
555 AliAODVertex *primaryVertex;
556 double vertex[3]={0};
557 primaryVertex = fAOD->GetPrimaryVertex();
558 vertex[0]=primaryVertex->GetX();
559 vertex[1]=primaryVertex->GetY();
560 vertex[2]=primaryVertex->GetZ();
561 if(vertex[0]<10e-5 && vertex[1]<10e-5 && vertex[2]<10e-5) return;
562 if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
563
18d901d5
DRG		 564 int zBin=0;
565 double zStep=2*10/double(kZVertexBins), zStart=-10.;
41dfc4d3
DRG		 566 for(int i=0; i<kZVertexBins; i++)
567 {
18d901d5
DRG		 568 if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
569 {
570 zBin=i;
571 break;
572 }
41dfc4d3
DRG		 573 }
574
575 //Event plane
576 int psiBin = 0;
577 AliEventplane *eventplane = fAOD->GetEventplane();
578 if(fPsiBinning && !eventplane) return;
579 double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
41dfc4d3
DRG		 580 ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
581
18d901d5
DRG		 582 double psiStep = PI/double(kPsiBins);
583 double psiStart = -0.5*PI;
584 for(int i=0; i<kPsiBins; i++)
585 {
586 if((psiEP > psiStart+i*psiStep) && (psiEP < psiStart+(i+1)*psiStep))
587 {
588 psiBin = i;
589 break;
590 }
591 }
592 if(!fPsiBinning) psiBin = 0;
593
41dfc4d3
DRG		 594////////////////////////////////////////////////////////////////
595//Cut Values and constants
596
597 //const bool kMCCase = kFALSE; //switch for MC analysis
598 const int kMaxNumK0 = 300; //maximum number of K0s, array size
599 const float kMinDCAPrimaryPion = 0.4; //minimum dca of pions to primary
600 const float kMaxDCADaughtersK0 = 0.3; //maximum dca of pions to each other - 3D
601 const float kMaxDCAK0 = 0.3; //maximum dca of K0 to primary
602 const float kMaxDLK0 = 30.0; //maximum decay length of K0
603 const float kMinDLK0 = fMinDecayLength; //minimum decay length of K0
604 const float kEtaCut = 0.8; //maximum |pseudorapidity|
605 const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle
606
607 const float kMinSeparation = fMinSep; //minimum daughter (pair) separation
608
609 const float kTOFLow = 0.8; //boundary for TOF usage
610 const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas
611 const float kMaxTPCSigmaPion = 3.0; //TPC # of sigmas
612
613 //const float kMassPion = .13957;
614 const float kMassK0Short = .497614; //true PDG masses
615
616////////////////////////////////////////////////////////////////
617 //v0 tester
618////////////////////////////////////////////////////////////////
619 int v0Count = 0; //number of v0s (entries in array)
620 int k0Count = 0; //number of good K0s
621
622 AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
623
624 //for daughter sharing studies
625 //int idArray[100] = {0};
626 //int idCount = 0;
627
628 //for MC
629 //TClonesArray *mcArray = 0x0;
630 //if(kMCCase){
631 //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
632 //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
633
634 for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
635 {
636 bool goodK0 = kFALSE;
637 bool goodPiPlus = kFALSE;
638 bool goodPiMinus = kFALSE;
639
640 //load v0 track
641 AliAODv0* v0 = fAOD->GetV0(i);
642 if(!v0) continue;
643 if(fOnlineCase){
644 if(!(v0->GetOnFlyStatus())) continue;
645 } //for online
646 else{
647 if((v0->GetOnFlyStatus())) continue; //for offline
648 }
649
650 //for on-the-fly ordering
651 AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
652 short int pos0or1;
653 short int neg0or1;
654 bool orderswitch = kFALSE;
655 if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
656 else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
657
658 //load daughter tracks
659 AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
660 AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
661 if(!prongTrackPos) continue;
662 if(!prongTrackNeg) continue;
663
664 //daughter cuts
665 if(v0->PtProng(pos0or1) < .15) continue;
666 if(v0->PtProng(neg0or1) < .15) continue;
667 if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
668 if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
669
670 //load status for PID
671 statusPos=prongTrackPos->GetStatus();
672 if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
673 prongTrackPos->SetAODEvent(fAOD);
674 statusNeg=prongTrackNeg->GetStatus();
675 if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
676 prongTrackNeg->SetAODEvent(fAOD);
677
678 //TPC PID
679 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
680 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
681
682 //Positive daughter identification TOF
683 float probMis;
684 AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
685 double Ppos = v0->PProng(pos0or1);
686 if(Ppos > kTOFLow) //PiPlus
687 {
688 //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
689 if(AliPIDResponse::kDetPidOk == statusPosTOF)
690 {
691 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
692 if(probMis < 0.01) //avoid TOF-TPC mismatch
693 {
694 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
695 else goodPiPlus = kFALSE;
696 }
697 }
698 }
699 //Negative daughter identification TOF
700 AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
701 double Pneg = v0->PProng(neg0or1);
702 if(Pneg > kTOFLow) //PiMinus
703 {
704 //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
705 if(AliPIDResponse::kDetPidOk == statusNegTOF)
706 {
707 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
708 if(probMis < 0.01) //avoid TOF-TPC mismatch
709 {
710 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
711 else goodPiMinus = kFALSE;
712 }
713 }
714 }
715
716 //K0 cuts
717 if(v0->Eta() > kEtaCut) continue;
718 if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;
719 if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;
720 if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;
721 if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue;
722 if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;
723 if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;
724 if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;
725 double v0Dca = v0->DcaV0ToPrimVertex();
726 if(v0Dca > kMaxDCAK0) continue;
727 if(!goodPiMinus || !goodPiPlus) continue;
728
729 //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
730
731 //for MC
732 //bool MCgood = kFALSE;
733 //if(kMCCase){
734 //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
735 //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
736 //if(mck0dp->GetMother() >= 0){
737 //if(mck0dp->GetMother() == mck0dn->GetMother()){
738 //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
739 //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
740 //if(abs(mck0->GetPdgCode()) == 310){
741 //MCgood = kTRUE;
742 //}
743 //}
744 //}
745 //}
746 //}// if kMCCase
747
748 if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
749 //else continue; //removed, Apr 18
750
751 //Check for shared daughters, using v0 DCA to judge
752 bool v0JudgeNew; //true if new v0 beats old
753 tempK0[v0Count].fSkipShared = kFALSE;
754 double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
755 if(fMeritCase){
756 for(int iID = 0; iID<v0Count; iID++){
757 if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old
758 if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
759 double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
760 v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
761 if(!v0JudgeNew){ //if old beats new...
762 if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing...
763 else{ //but if bad old beats new bad, or good old beats anything, skip new
764 tempK0[v0Count].fSkipShared = kTRUE; //skip new
765 break; //no need to keep checking others
766 }
767 }
768 else{ //if new beats old...
769 if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...
770 else{ //but if bad new beats bad old, or good new beats anything, skip old
771 tempK0[iID].fSkipShared = kTRUE; //skip old
772 if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
773 }
774 }
775 }
776 }
777 }
778 if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0
779 }//if MeritCase
780
781 //load parameters into temporary class instance
782 if(v0Count < kMaxNumK0)
783 {
784 if(goodK0){
785 tempK0[v0Count].fK0 = kTRUE;
786 k0Count++;
787 }
788 else tempK0[v0Count].fK0 = kFALSE;
789
790 //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
791 //else tempK0[v0Count].fSideLeft = kFALSE;
792 //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
793 //else tempK0[v0Count].fSideRight = kFALSE;
794 //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
795
796 tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();
797 tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();
798 tempK0[v0Count].fMomentum[0] = v0->Px();
799 tempK0[v0Count].fMomentum[1] = v0->Py();
800 tempK0[v0Count].fMomentum[2] = v0->Pz();
801 tempK0[v0Count].fPt = v0->Pt();
802 tempK0[v0Count].fMass = v0->MassK0Short();
803 tempK0[v0Count].fV0Dca = v0Dca;
804
805 //for hists
806 tempK0[v0Count].fDDDca = v0->DcaV0Daughters();
807 tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);
808 tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);
809 tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);
810 tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);
811 tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);
812 if(!orderswitch){
813 tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();
814 tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();
815 }
816 else{
817 tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();
818 tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();
819 }
820
821 //for psi studies
822 double v0Phi = v0->Phi(); //between [0,2pi]
823 double v0PhiPsi = v0Phi-psiEP;
824 if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
825 else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
826 else{};
827 tempK0[v0Count].fPhiPsi = v0PhiPsi;
828
829 //for separation
830 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
831 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
832 //for DPC
833 prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
834 prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
835
836 //if(idCount < 50){
837 // if(goodK0){
838 // idArray[idCount*2] = prongTrackPos->GetID();
839 // idArray[idCount*2+1] = prongTrackNeg->GetID();
840 // idCount++;
841 //}}
842
843 v0Count++;
844 }
845
846 }//v0
847
848 if(k0Count<2) return; //only keep events with more than 1 good K0
849
850 //Add Event to buffer - this is for event mixing
851 fEC[zBin][centBin][psiBin]->FIFOShift();
852 (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
853 (fEvt)->fFillStatus = 1;
854 int unskippedCount = 0;
855 for(int i=0;i<v0Count;i++)
856 {
857 if(!tempK0[i].fSkipShared) //don't include skipped v0s (from shared daughters)
858 {
859 ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
860 if(tempK0[i].fK0) //make sure particle is good (mass)
861 {
862 (fEvt)->fK0Particle[unskippedCount] = tempK0[i]; //load good, unskipped K0s
863 unskippedCount++; //count good, unskipped K0s
864 }
865 }
866 }
867 (fEvt)->fNumV0s = unskippedCount;
868 //Printf("Number of v0s: %d", v0Count);
869 //Printf("Number of K0s: %d", k0Count);
870 delete [] tempK0; tempK0 = NULL;
871
872 ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount); // changed 3/25, used to be "k0Count"
873 ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
874
875 //Printf("Reconstruction Finished. Starting pair studies.");
876
877 //////////////////////////////////////////////////////////////////////
878 // Correlations
879 //////////////////////////////////////////////////////////////////////
880
881 float px1, py1, pz1, px2, py2, pz2; //single kaon values
882 float en1, en2; //single kaon values
883 //float pt1, pt2; //single kaon values
884 float pairPx, pairPy, pairPz, pairP0; //pair momentum values
885 float pairPt, pairMt, pairKt; //pair momentum values
886 float pairMInv, pairPDotQ;
887 float qinv, q0, qx, qy, qz; //pair q values
888 //float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values
889 float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF
890 //float qOutLCMS;
891 float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame
892 float betasq, gamma;
893 float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
894 float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
895
896
897 for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
898 {
899 //single particle histograms (done here to avoid "skipped" v0s
900 ((TH1F*)fOutputList->FindObject("fHistDCADaughters")) ->Fill((fEvt)->fK0Particle[i].fDDDca);
901 ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0")) ->Fill((fEvt)->fK0Particle[i].fDecayLength);
902 ((TH1F*)fOutputList->FindObject("fHistDCAK0")) ->Fill((fEvt)->fK0Particle[i].fV0Dca);
903 ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
904 ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus")) ->Fill((fEvt)->fK0Particle[i].fPosDca);
905 ((TH2F*)fOutputList->FindObject("fHistPtK0")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
906 ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
907 ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
908 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fPosPhi);
909 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fNegPhi);
910
911 ((TH1F*)fOutputList->FindObject("fHistPx")) ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
912 ((TH1F*)fOutputList->FindObject("fHistPy")) ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
913 ((TH1F*)fOutputList->FindObject("fHistPz")) ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
914
915 ((TH2F*)fOutputList->FindObject("fHistPhiPsi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
916
917 for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
918 {
919 int startbin=0;
920 if(evnum==0) startbin=i+1;
921
922 for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
923 {
924 if(evnum==0) // Get rid of shared tracks
925 {
926 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
927 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
928 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
929 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
930 }
931
932 px1 = (fEvt)->fK0Particle[i].fMomentum[0];
933 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
934 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
935 //pt1 = (fEvt)->fK0Particle[i].fPt;
936 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
937 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
938 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
939 //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
940 if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
941 double tempvar;
942 tempvar = px1; px1 = px2; px2 = tempvar;
943 tempvar = py1; py1 = py2; py2 = tempvar;
944 tempvar = pz1; pz1 = pz2; pz2 = tempvar;
945 }
946
947 en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
948 en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
949
950 q0 = en1 - en2;
951 qx = px1 - px2;
952 qy = py1 - py2;
953 qz = pz1 - pz2;
954 qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
955
956 pairPx = px1 + px2;
957 pairPy = py1 + py2;
958 pairPz = pz1 + pz2;
959 pairP0 = en1 + en2;
960 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
961 pairKt = pairPt/2.; //used for KT binning
962 pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)
963 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
964 pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF
965
966 //PRF (this section will probably be removed in favor of later boosting section)
967 p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length
968 am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)
969 epm = en1+en2+am12; //"energy plus mass"
970 p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP
971 if(am12 == 0) continue;
972 h1 = (p112/epm - en1)/am12;
973 ppx = px1+pairPx*h1; //px in PRF
974 ppy = py1+pairPy*h1; //py in PRF
975 ppz = pz1+pairPz*h1; //pz in PRF
976 ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*
977 ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
978 ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
979 ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
980 ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
981
982 //relative momentum in out-side-long for LCMS and PRF
983 if(pairMt == 0 || pairPt == 0) continue;
984 qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
985 qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
986 //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
987 qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
988
989 //finding gamma for gamma binning/hists (likely will be removed after tests)
990 p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;
991 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
992 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
993 p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;
994 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
995 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
996 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
997 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
998 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
999 gamma = 1./sqrt(1-betasq);
1000 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
1001 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
1002 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
1003 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
1004 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
1005 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1006 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1007 //getting bin numbers and names for 3D histogram
1008 TString *histname3D = new TString("fHist3DOSL");
1009 int ktBin;
1010 if(pairKt < 0.6) ktBin = 0;
1011 else if(pairKt < 0.8) ktBin = 1;
1012 else if(pairKt < 1.0) ktBin = 2;
1013 else ktBin = 3;
1014 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1015 *histname3D += ktBin;
1016
1017 //Spherical harmonics
1018 //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1019 //thetaSHCos = qLong/qLength;
1020 //thetaSH = acos(thetaSHCos);
1021 //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1022
1023 //Finding average separation of daughters throughout TPC - two-track cut
1024 float posPositions1[9][3] = {{0}};
1025 float negPositions1[9][3] = {{0}};
1026 float posPositions2[9][3] = {{0}};
1027 float negPositions2[9][3] = {{0}};
1028 for(int iPos = 0; iPos < 9; iPos++){
1029 for(int jPos = 0; jPos < 3; jPos++){
1030 posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1031 negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1032 posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1033 negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1034 }
1035 }
1036 float pMean = 0.; //average separation for positive daughters
1037 float nMean = 0.; //average separation for negative daughters
1038 float pDiff;
1039 float nDiff;
1040 float pMin = 9999.; //minimum separation (updates) - pos
1041 float nMin = 9999.; //minimum separation (updates) - neg
1042 double pCount=0; //counter for number of points used - pos
1043 double nCount=0; //counter for number of points used - neg
1044 for(int ss=0;ss<9;ss++){
1045 if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1046 pCount++;
1047 pDiff = sqrt(pow(posPositions1[ss][0]-posPositions2[ss][0],2)+pow(posPositions1[ss][1]-posPositions2[ss][1],2)+pow(posPositions1[ss][2]-posPositions2[ss][2],2));
1048 pMean = pMean + pDiff;
1049 if(pDiff < pMin) pMin = pDiff;
1050 }
1051 if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1052 nCount++;
1053 nDiff = sqrt(pow(negPositions1[ss][0]-negPositions2[ss][0],2)+pow(negPositions1[ss][1]-negPositions2[ss][1],2)+pow(negPositions1[ss][2]-negPositions2[ss][2],2));
1054 nMean = nMean + nDiff;
1055 if(nDiff < nMin) nMin = nDiff;
1056 }
1057 }
1058 pMean = pMean/pCount;
1059 nMean = nMean/nCount;
1060
1061 if(evnum==0){
1062 ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1063 ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1064 ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1065 ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1066 }
1067 else{
1068 ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1069 ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1070 ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1071 ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1072 }
1073
1074 //Decay plane coincidence
1075 //daughter momenta
1076 float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1077 float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1078 float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1079 float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1080 float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1081 float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1082 float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1083 float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1084 float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1085 float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1086 float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1087 float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1088
1089 float cross1[3];
1090 float cross2[3];
1091 cross1[0] = b1*f1-c1*e1;
1092 cross1[1] = c1*d1-a1*f1;
1093 cross1[2] = a1*e1-b1*d1;
1094 cross2[0] = b2*f2-c2*e2;
1095 cross2[1] = c2*d2-a2*f2;
1096 cross2[2] = a2*e2-b2*d2;
1097 float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1098 float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1099 float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1100
1101 if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1102 else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1103
1104 if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1105 //end separation studies
1106
1107 //Fill Histograms
1108 bool center1K0 = kFALSE; //accepted mass K0
1109 bool center2K0 = kFALSE;
1110 if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1111 if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1112 //bool CL1 = kFALSE;
1113 //bool CL2 = kFALSE;
1114 //bool CR1 = kFALSE;
1115 //bool CR2 = kFALSE;
1116 //if(center1K0 && center2K0){
1117 // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1118 // else CR1 = kTRUE;
1119 // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1120 // else CR2 = kTRUE;
1121 //}
1122
1123 //bool SideLeft1 = kFALSE;
1124 //bool SideLeft2 = kFALSE;
1125 //bool SideRight1 = kFALSE;
1126 //bool SideRight2 = kFALSE;
1127 //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1128 //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1129 //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1130 //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1131
1132 //for psi binning
1133 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1134 float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1135 bool inPlane1 = kFALSE;
1136 bool inPlane2 = kFALSE;
1137 if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1138 if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1139 if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1140 if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
1141
1142
1143 if(evnum==0) //Same Event
1144 {
1145 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1146 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1147 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1148 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1149 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1150 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1151
1152 if(center1K0 && center2K0){
1153 //1D
1154 ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1155 //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1156 ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1157
1158 //eventplane
1159 if(inPlane1 && inPlane2)
1160 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1161 else if(!inPlane1 && !inPlane2)
1162 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
1163
1164 //for mass bin study
1165 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1166 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1167 //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1168
1169 //3D
1170 if(fCase3D){
1171 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1172 histname3D->Append("Signal");
1173 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1174 }
1175 }
1176 /*if(pairKt < 1.0){
1177 if(centBin > 13){
1178 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1179 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1180 else if(centBin > 9){
1181 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1182 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1183 else if(pairKt < 2.0){
1184 if(centBin > 13){
1185 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1186 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1187 else if(centBin > 9){
1188 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1189
1190 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1191
1192 }//centercenter
1193
1194 //side-side correlations
1195 //if(!splitK0sides){
1196 // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1197 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1198 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1199 //}
1200
1201 }//same event
1202
1203 else //Mixed Events
1204 {
1205 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1206 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1207 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1208 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1209
1210 if(center1K0 && center2K0){
1211 //1D
1212 ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1213
1214 //eventplane
1215 if(inPlane1 && inPlane2)
1216 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1217 else if(!inPlane1 && !inPlane2)
1218 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1219
1220 //for mass bin study
1221 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1222 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1223 //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1224
1225 //3D
1226 if(fCase3D){
1227 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1228 histname3D->Replace(12,6,"Bkg");
1229 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1230 }
1231 }
1232 /*if(pairKt < 1.0){
1233 if(centBin > 13){
1234 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1235 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1236 else if(centBin > 9){
1237 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1238 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1239 else if(pairKt < 2.0){
1240 if(centBin > 13){
1241 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1242 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1243 else if(centBin > 9){
1244 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1245 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1246 }
1247
1248 //side-side correlations
1249 //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1250 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1251 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1252
1253 }//Mixed Events
1254
1255 }//past event
1256 }//event buffer
1257 }//current event
1258
1259 // Post output data.
1260 PostData(1, fOutputList);
1261
1262 }
1263//________________________________________________________________________
1264void AliFemtoK0Analysis::Terminate(Option_t *)
1265{
1266 // Called once at the end of the query
1267 cout<<"Done"<<endl;
1268
1269}
1270
1271//_________________________________________________________________________
1272void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1273 // Gets the global position of the track at nine different radii in the TPC
1274 // track is the track you want to propagate
1275 // bfield is the magnetic field of your event
1276 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1277
1278 // Initialize the array to something indicating there was no propagation
1279 for(Int_t i=0;i<9;i++){
1280 for(Int_t j=0;j<3;j++){
1281 globalPositionsAtRadii[i][j]=-9999.;
1282 }
1283 }
1284
1285 // Make a copy of the track to not change parameters of the track
1286 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1287 //printf("\nAfter CopyFromVTrack\n");
1288 //etp.Print();
1289
1290 // The global position of the the track
1291 Double_t xyz[3]={-9999.,-9999.,-9999.};
1292
1293 // Counter for which radius we want
1294 Int_t iR=0;
1295 // The radii at which we get the global positions
1296 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1297 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1298 // The global radius we are at
1299 Float_t globalRadius=0;
1300
1301 // Propagation is done in local x of the track
1302 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1303 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1304 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1305 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1306 // adds to the global radius
1307
1308 // Stop if the propagation was not succesful. This can happen for low pt tracks
1309 // that don't reach outer radii
1310 if(!etp.PropagateTo(x,bfield))break;
1311 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1312 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1313
1314 // Roughly reached the radius we want
1315 if(globalRadius > Rwanted[iR]){
1316
1317 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1318 while (globalRadius>Rwanted[iR]){
1319 x-=.1;
1320 // printf("propagating to x %5.2f\n",x);
1321 if(!etp.PropagateTo(x,bfield))break;
1322 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1323 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1324 }
1325 //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
1326 globalPositionsAtRadii[iR][0]=xyz[0];
1327 globalPositionsAtRadii[iR][1]=xyz[1];
1328 globalPositionsAtRadii[iR][2]=xyz[2];
1329 //subtract primary vertex, "zero" track for correct mixed-event comparison
1330 globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1331 globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1332 globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1333
1334 // Indicate we want the next radius
1335 iR+=1;
1336 }
1337 if(iR>=8){
1338 // TPC edge reached
1339 return;
1340 }
1341 }
1342}
1343
1344bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1345 //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1346 //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1347
1348 bool newV0Wins = kFALSE;
1349 double pardiff[3] = {newPars[0]-oldPars[0],
1350 newPars[1]-oldPars[1],
1351 newPars[2]-oldPars[2]};
1352 if(cutChoice > 0 && cutChoice < 4){
1353 if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1354 }
1355 else if(cutChoice == 4){
1356 int newWinCount = 0;
1357 for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1358 if(newWinCount > 1) newV0Wins = kTRUE;
1359 }
1360 else{};
1361 return newV0Wins;
1362}
1363
1364bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1365{ // to flatten centrality distribution
1366 bool RejectEvent = kFALSE;
1367 int weightBinSign;
1368 if(MagField > 0) weightBinSign = 0;
1369 else weightBinSign = 1;
1370 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1371 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1372 int weightBinCent = (int) CentPercent;
1373 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1374
1375 return RejectEvent;
1376}