K0s code update (Matt Steinpreis)
[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
501 int zBin=0;
502 double zStep=2*10/double(kZVertexBins), zstart=-10.;
503
504 /////////////////////////////////////////////////
505
506
507 //Centrality selection
508
509 AliCentrality *centrality = fAOD->GetCentrality();
510 float percent = centrality->GetCentralityPercentile("V0M");
511 int centBin=0;
512 //Printf("Centrality percent = %f", percent);
513
514 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
515 float multV0A=aodV0->GetMTotV0A();
516 float multV0C=aodV0->GetMTotV0C();
517
518 if(percent < 0) {
519 //Printf("No centrality info");
520 return;
521 }
522 if(percent < 0.1 && (multV0A + multV0C < 19500)){
523 //Printf("No centrality info");
524 return;
525 }
526 else if(percent <= 5) centBin=15;
527 else if(percent <= 10) centBin=14;
528 else if(percent <= 15) centBin=13;
529 else if(percent <= 20) centBin=12;
530 else if(percent <= 25) centBin=11;
531 else if(percent <= 30) centBin=10;
532 else if(percent <= 35) centBin=9;
533 else if(percent <= 40) centBin=8;
534 else if(percent <= 45) centBin=7;
535 else if(percent <= 50) centBin=6;
536 else if(percent <= 55) centBin=5;
537 else if(percent <= 60) centBin=4;
538 else if(percent <= 65) centBin=3;
539 else if(percent <= 70) centBin=2;
540 else if(percent <= 75) centBin=1;
541 else if(percent <= 80) centBin=0;
542 else {
543 //Printf("Skipping Peripheral Event");
544 return;
545 }
546 if(percent > 10 && isCentral) return;
547 ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
548
549 //flatten centrality dist.
550 if(percent < 9){
551 if(fFlatCent){
552 if(RejectEventCentFlat(bField,percent)) return;
553 }
554 }
555 ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
556
557 //Vertexing
558 AliAODVertex *primaryVertex;
559 double vertex[3]={0};
560 primaryVertex = fAOD->GetPrimaryVertex();
561 vertex[0]=primaryVertex->GetX();
562 vertex[1]=primaryVertex->GetY();
563 vertex[2]=primaryVertex->GetZ();
564 if(vertex[0]<10e-5 && vertex[1]<10e-5 && vertex[2]<10e-5) return;
565 if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
566
567 for(int i=0; i<kZVertexBins; i++)
568 {
569 if((vertex[2] > zstart+i*zStep) && (vertex[2] < zstart+(i+1)*zStep))
570 {
571 zBin=i;
572 break;
573 }
574 }
575
576 //Event plane
577 int psiBin = 0;
578 AliEventplane *eventplane = fAOD->GetEventplane();
579 if(fPsiBinning && !eventplane) return;
580 double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
581 if(psiEP < -0.25*PI) psiBin = 0;
582 else if(psiEP < 0.0) psiBin = 1;
583 else if(psiEP < 0.25*PI) psiBin = 2;
584 else psiBin = 3;
585 if(!fPsiBinning) psiBin = 0;
586 ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
587
588////////////////////////////////////////////////////////////////
589//Cut Values and constants
590
591 //const bool kMCCase = kFALSE; //switch for MC analysis
592 const int kMaxNumK0 = 300; //maximum number of K0s, array size
593 const float kMinDCAPrimaryPion = 0.4; //minimum dca of pions to primary
594 const float kMaxDCADaughtersK0 = 0.3; //maximum dca of pions to each other - 3D
595 const float kMaxDCAK0 = 0.3; //maximum dca of K0 to primary
596 const float kMaxDLK0 = 30.0; //maximum decay length of K0
597 const float kMinDLK0 = fMinDecayLength; //minimum decay length of K0
598 const float kEtaCut = 0.8; //maximum |pseudorapidity|
599 const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle
600
601 const float kMinSeparation = fMinSep; //minimum daughter (pair) separation
602
603 const float kTOFLow = 0.8; //boundary for TOF usage
604 const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas
605 const float kMaxTPCSigmaPion = 3.0; //TPC # of sigmas
606
607 //const float kMassPion = .13957;
608 const float kMassK0Short = .497614; //true PDG masses
609
610////////////////////////////////////////////////////////////////
611 //v0 tester
612////////////////////////////////////////////////////////////////
613 int v0Count = 0; //number of v0s (entries in array)
614 int k0Count = 0; //number of good K0s
615
616 AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
617
618 //for daughter sharing studies
619 //int idArray[100] = {0};
620 //int idCount = 0;
621
622 //for MC
623 //TClonesArray *mcArray = 0x0;
624 //if(kMCCase){
625 //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
626 //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
627
628 for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
629 {
630 bool goodK0 = kFALSE;
631 bool goodPiPlus = kFALSE;
632 bool goodPiMinus = kFALSE;
633
634 //load v0 track
635 AliAODv0* v0 = fAOD->GetV0(i);
636 if(!v0) continue;
637 if(fOnlineCase){
638 if(!(v0->GetOnFlyStatus())) continue;
639 } //for online
640 else{
641 if((v0->GetOnFlyStatus())) continue; //for offline
642 }
643
644 //for on-the-fly ordering
645 AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
646 short int pos0or1;
647 short int neg0or1;
648 bool orderswitch = kFALSE;
649 if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
650 else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
651
652 //load daughter tracks
653 AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
654 AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
655 if(!prongTrackPos) continue;
656 if(!prongTrackNeg) continue;
657
658 //daughter cuts
659 if(v0->PtProng(pos0or1) < .15) continue;
660 if(v0->PtProng(neg0or1) < .15) continue;
661 if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
662 if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
663
664 //load status for PID
665 statusPos=prongTrackPos->GetStatus();
666 if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
667 prongTrackPos->SetAODEvent(fAOD);
668 statusNeg=prongTrackNeg->GetStatus();
669 if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
670 prongTrackNeg->SetAODEvent(fAOD);
671
672 //TPC PID
673 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
674 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
675
676 //Positive daughter identification TOF
677 float probMis;
678 AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
679 double Ppos = v0->PProng(pos0or1);
680 if(Ppos > kTOFLow) //PiPlus
681 {
682 //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
683 if(AliPIDResponse::kDetPidOk == statusPosTOF)
684 {
685 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
686 if(probMis < 0.01) //avoid TOF-TPC mismatch
687 {
688 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
689 else goodPiPlus = kFALSE;
690 }
691 }
692 }
693 //Negative daughter identification TOF
694 AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
695 double Pneg = v0->PProng(neg0or1);
696 if(Pneg > kTOFLow) //PiMinus
697 {
698 //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
699 if(AliPIDResponse::kDetPidOk == statusNegTOF)
700 {
701 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
702 if(probMis < 0.01) //avoid TOF-TPC mismatch
703 {
704 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
705 else goodPiMinus = kFALSE;
706 }
707 }
708 }
709
710 //K0 cuts
711 if(v0->Eta() > kEtaCut) continue;
712 if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;
713 if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;
714 if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;
715 if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue;
716 if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;
717 if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;
718 if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;
719 double v0Dca = v0->DcaV0ToPrimVertex();
720 if(v0Dca > kMaxDCAK0) continue;
721 if(!goodPiMinus || !goodPiPlus) continue;
722
723 //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
724
725 //for MC
726 //bool MCgood = kFALSE;
727 //if(kMCCase){
728 //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
729 //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
730 //if(mck0dp->GetMother() >= 0){
731 //if(mck0dp->GetMother() == mck0dn->GetMother()){
732 //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
733 //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
734 //if(abs(mck0->GetPdgCode()) == 310){
735 //MCgood = kTRUE;
736 //}
737 //}
738 //}
739 //}
740 //}// if kMCCase
741
742 if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
743 //else continue; //removed, Apr 18
744
745 //Check for shared daughters, using v0 DCA to judge
746 bool v0JudgeNew; //true if new v0 beats old
747 tempK0[v0Count].fSkipShared = kFALSE;
748 double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
749 if(fMeritCase){
750 for(int iID = 0; iID<v0Count; iID++){
751 if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old
752 if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
753 double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
754 v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
755 if(!v0JudgeNew){ //if old beats new...
756 if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing...
757 else{ //but if bad old beats new bad, or good old beats anything, skip new
758 tempK0[v0Count].fSkipShared = kTRUE; //skip new
759 break; //no need to keep checking others
760 }
761 }
762 else{ //if new beats old...
763 if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...
764 else{ //but if bad new beats bad old, or good new beats anything, skip old
765 tempK0[iID].fSkipShared = kTRUE; //skip old
766 if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
767 }
768 }
769 }
770 }
771 }
772 if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0
773 }//if MeritCase
774
775 //load parameters into temporary class instance
776 if(v0Count < kMaxNumK0)
777 {
778 if(goodK0){
779 tempK0[v0Count].fK0 = kTRUE;
780 k0Count++;
781 }
782 else tempK0[v0Count].fK0 = kFALSE;
783
784 //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
785 //else tempK0[v0Count].fSideLeft = kFALSE;
786 //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
787 //else tempK0[v0Count].fSideRight = kFALSE;
788 //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
789
790 tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();
791 tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();
792 tempK0[v0Count].fMomentum[0] = v0->Px();
793 tempK0[v0Count].fMomentum[1] = v0->Py();
794 tempK0[v0Count].fMomentum[2] = v0->Pz();
795 tempK0[v0Count].fPt = v0->Pt();
796 tempK0[v0Count].fMass = v0->MassK0Short();
797 tempK0[v0Count].fV0Dca = v0Dca;
798
799 //for hists
800 tempK0[v0Count].fDDDca = v0->DcaV0Daughters();
801 tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);
802 tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);
803 tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);
804 tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);
805 tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);
806 if(!orderswitch){
807 tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();
808 tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();
809 }
810 else{
811 tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();
812 tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();
813 }
814
815 //for psi studies
816 double v0Phi = v0->Phi(); //between [0,2pi]
817 double v0PhiPsi = v0Phi-psiEP;
818 if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
819 else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
820 else{};
821 tempK0[v0Count].fPhiPsi = v0PhiPsi;
822
823 //for separation
824 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
825 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
826 //for DPC
827 prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
828 prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
829
830 //if(idCount < 50){
831 // if(goodK0){
832 // idArray[idCount*2] = prongTrackPos->GetID();
833 // idArray[idCount*2+1] = prongTrackNeg->GetID();
834 // idCount++;
835 //}}
836
837 v0Count++;
838 }
839
840 }//v0
841
842 if(k0Count<2) return; //only keep events with more than 1 good K0
843
844 //Add Event to buffer - this is for event mixing
845 fEC[zBin][centBin][psiBin]->FIFOShift();
846 (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
847 (fEvt)->fFillStatus = 1;
848 int unskippedCount = 0;
849 for(int i=0;i<v0Count;i++)
850 {
851 if(!tempK0[i].fSkipShared) //don't include skipped v0s (from shared daughters)
852 {
853 ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
854 if(tempK0[i].fK0) //make sure particle is good (mass)
855 {
856 (fEvt)->fK0Particle[unskippedCount] = tempK0[i]; //load good, unskipped K0s
857 unskippedCount++; //count good, unskipped K0s
858 }
859 }
860 }
861 (fEvt)->fNumV0s = unskippedCount;
862 //Printf("Number of v0s: %d", v0Count);
863 //Printf("Number of K0s: %d", k0Count);
864 delete [] tempK0; tempK0 = NULL;
865
866 ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount); // changed 3/25, used to be "k0Count"
867 ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
868
869 //Printf("Reconstruction Finished. Starting pair studies.");
870
871 //////////////////////////////////////////////////////////////////////
872 // Correlations
873 //////////////////////////////////////////////////////////////////////
874
875 float px1, py1, pz1, px2, py2, pz2; //single kaon values
876 float en1, en2; //single kaon values
877 //float pt1, pt2; //single kaon values
878 float pairPx, pairPy, pairPz, pairP0; //pair momentum values
879 float pairPt, pairMt, pairKt; //pair momentum values
880 float pairMInv, pairPDotQ;
881 float qinv, q0, qx, qy, qz; //pair q values
882 //float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values
883 float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF
884 //float qOutLCMS;
885 float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame
886 float betasq, gamma;
887 float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
888 float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
889
890
891 for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
892 {
893 //single particle histograms (done here to avoid "skipped" v0s
894 ((TH1F*)fOutputList->FindObject("fHistDCADaughters")) ->Fill((fEvt)->fK0Particle[i].fDDDca);
895 ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0")) ->Fill((fEvt)->fK0Particle[i].fDecayLength);
896 ((TH1F*)fOutputList->FindObject("fHistDCAK0")) ->Fill((fEvt)->fK0Particle[i].fV0Dca);
897 ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
898 ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus")) ->Fill((fEvt)->fK0Particle[i].fPosDca);
899 ((TH2F*)fOutputList->FindObject("fHistPtK0")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
900 ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
901 ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
902 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fPosPhi);
903 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fNegPhi);
904
905 ((TH1F*)fOutputList->FindObject("fHistPx")) ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
906 ((TH1F*)fOutputList->FindObject("fHistPy")) ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
907 ((TH1F*)fOutputList->FindObject("fHistPz")) ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
908
909 ((TH2F*)fOutputList->FindObject("fHistPhiPsi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
910
911 for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
912 {
913 int startbin=0;
914 if(evnum==0) startbin=i+1;
915
916 for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
917 {
918 if(evnum==0) // Get rid of shared tracks
919 {
920 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
921 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
922 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
923 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
924 }
925
926 px1 = (fEvt)->fK0Particle[i].fMomentum[0];
927 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
928 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
929 //pt1 = (fEvt)->fK0Particle[i].fPt;
930 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
931 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
932 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
933 //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
934 if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
935 double tempvar;
936 tempvar = px1; px1 = px2; px2 = tempvar;
937 tempvar = py1; py1 = py2; py2 = tempvar;
938 tempvar = pz1; pz1 = pz2; pz2 = tempvar;
939 }
940
941 en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
942 en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
943
944 q0 = en1 - en2;
945 qx = px1 - px2;
946 qy = py1 - py2;
947 qz = pz1 - pz2;
948 qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
949
950 pairPx = px1 + px2;
951 pairPy = py1 + py2;
952 pairPz = pz1 + pz2;
953 pairP0 = en1 + en2;
954 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
955 pairKt = pairPt/2.; //used for KT binning
956 pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)
957 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
958 pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF
959
960 //PRF (this section will probably be removed in favor of later boosting section)
961 p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length
962 am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)
963 epm = en1+en2+am12; //"energy plus mass"
964 p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP
965 if(am12 == 0) continue;
966 h1 = (p112/epm - en1)/am12;
967 ppx = px1+pairPx*h1; //px in PRF
968 ppy = py1+pairPy*h1; //py in PRF
969 ppz = pz1+pairPz*h1; //pz in PRF
970 ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*
971 ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
972 ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
973 ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
974 ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
975
976 //relative momentum in out-side-long for LCMS and PRF
977 if(pairMt == 0 || pairPt == 0) continue;
978 qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
979 qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
980 //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
981 qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
982
983 //finding gamma for gamma binning/hists (likely will be removed after tests)
984 p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;
985 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
986 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
987 p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;
988 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
989 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
990 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
991 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
992 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
993 gamma = 1./sqrt(1-betasq);
994 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
995 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
996 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
997 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
998 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
999 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1000 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1001 //getting bin numbers and names for 3D histogram
1002 TString *histname3D = new TString("fHist3DOSL");
1003 int ktBin;
1004 if(pairKt < 0.6) ktBin = 0;
1005 else if(pairKt < 0.8) ktBin = 1;
1006 else if(pairKt < 1.0) ktBin = 2;
1007 else ktBin = 3;
1008 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1009 *histname3D += ktBin;
1010
1011 //Spherical harmonics
1012 //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1013 //thetaSHCos = qLong/qLength;
1014 //thetaSH = acos(thetaSHCos);
1015 //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1016
1017 //Finding average separation of daughters throughout TPC - two-track cut
1018 float posPositions1[9][3] = {{0}};
1019 float negPositions1[9][3] = {{0}};
1020 float posPositions2[9][3] = {{0}};
1021 float negPositions2[9][3] = {{0}};
1022 for(int iPos = 0; iPos < 9; iPos++){
1023 for(int jPos = 0; jPos < 3; jPos++){
1024 posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1025 negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1026 posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1027 negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1028 }
1029 }
1030 float pMean = 0.; //average separation for positive daughters
1031 float nMean = 0.; //average separation for negative daughters
1032 float pDiff;
1033 float nDiff;
1034 float pMin = 9999.; //minimum separation (updates) - pos
1035 float nMin = 9999.; //minimum separation (updates) - neg
1036 double pCount=0; //counter for number of points used - pos
1037 double nCount=0; //counter for number of points used - neg
1038 for(int ss=0;ss<9;ss++){
1039 if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1040 pCount++;
1041 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));
1042 pMean = pMean + pDiff;
1043 if(pDiff < pMin) pMin = pDiff;
1044 }
1045 if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1046 nCount++;
1047 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));
1048 nMean = nMean + nDiff;
1049 if(nDiff < nMin) nMin = nDiff;
1050 }
1051 }
1052 pMean = pMean/pCount;
1053 nMean = nMean/nCount;
1054
1055 if(evnum==0){
1056 ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1057 ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1058 ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1059 ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1060 }
1061 else{
1062 ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1063 ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1064 ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1065 ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1066 }
1067
1068 //Decay plane coincidence
1069 //daughter momenta
1070 float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1071 float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1072 float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1073 float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1074 float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1075 float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1076 float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1077 float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1078 float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1079 float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1080 float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1081 float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1082
1083 float cross1[3];
1084 float cross2[3];
1085 cross1[0] = b1*f1-c1*e1;
1086 cross1[1] = c1*d1-a1*f1;
1087 cross1[2] = a1*e1-b1*d1;
1088 cross2[0] = b2*f2-c2*e2;
1089 cross2[1] = c2*d2-a2*f2;
1090 cross2[2] = a2*e2-b2*d2;
1091 float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1092 float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1093 float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1094
1095 if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1096 else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1097
1098 if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1099 //end separation studies
1100
1101 //Fill Histograms
1102 bool center1K0 = kFALSE; //accepted mass K0
1103 bool center2K0 = kFALSE;
1104 if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1105 if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1106 //bool CL1 = kFALSE;
1107 //bool CL2 = kFALSE;
1108 //bool CR1 = kFALSE;
1109 //bool CR2 = kFALSE;
1110 //if(center1K0 && center2K0){
1111 // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1112 // else CR1 = kTRUE;
1113 // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1114 // else CR2 = kTRUE;
1115 //}
1116
1117 //bool SideLeft1 = kFALSE;
1118 //bool SideLeft2 = kFALSE;
1119 //bool SideRight1 = kFALSE;
1120 //bool SideRight2 = kFALSE;
1121 //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1122 //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1123 //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1124 //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1125
1126 //for psi binning
1127 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1128 float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1129 bool inPlane1 = kFALSE;
1130 bool inPlane2 = kFALSE;
1131 if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1132 if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1133 if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1134 if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
1135
1136
1137 if(evnum==0) //Same Event
1138 {
1139 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1140 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1141 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1142 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1143 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1144 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1145
1146 if(center1K0 && center2K0){
1147 //1D
1148 ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1149 //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1150 ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1151
1152 //eventplane
1153 if(inPlane1 && inPlane2)
1154 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1155 else if(!inPlane1 && !inPlane2)
1156 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
1157
1158 //for mass bin study
1159 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1160 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1161 //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1162
1163 //3D
1164 if(fCase3D){
1165 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1166 histname3D->Append("Signal");
1167 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1168 }
1169 }
1170 /*if(pairKt < 1.0){
1171 if(centBin > 13){
1172 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1173 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1174 else if(centBin > 9){
1175 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1176 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1177 else if(pairKt < 2.0){
1178 if(centBin > 13){
1179 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1180 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1181 else if(centBin > 9){
1182 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1183
1184 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1185
1186 }//centercenter
1187
1188 //side-side correlations
1189 //if(!splitK0sides){
1190 // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1191 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1192 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1193 //}
1194
1195 }//same event
1196
1197 else //Mixed Events
1198 {
1199 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1200 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1201 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1202 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1203
1204 if(center1K0 && center2K0){
1205 //1D
1206 ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1207
1208 //eventplane
1209 if(inPlane1 && inPlane2)
1210 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1211 else if(!inPlane1 && !inPlane2)
1212 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1213
1214 //for mass bin study
1215 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1216 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1217 //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1218
1219 //3D
1220 if(fCase3D){
1221 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1222 histname3D->Replace(12,6,"Bkg");
1223 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1224 }
1225 }
1226 /*if(pairKt < 1.0){
1227 if(centBin > 13){
1228 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1229 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1230 else if(centBin > 9){
1231 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1232 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1233 else if(pairKt < 2.0){
1234 if(centBin > 13){
1235 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1236 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1237 else if(centBin > 9){
1238 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1239 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1240 }
1241
1242 //side-side correlations
1243 //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1244 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1245 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1246
1247 }//Mixed Events
1248
1249 }//past event
1250 }//event buffer
1251 }//current event
1252
1253 // Post output data.
1254 PostData(1, fOutputList);
1255
1256 }
1257//________________________________________________________________________
1258void AliFemtoK0Analysis::Terminate(Option_t *)
1259{
1260 // Called once at the end of the query
1261 cout<<"Done"<<endl;
1262
1263}
1264
1265//_________________________________________________________________________
1266void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1267 // Gets the global position of the track at nine different radii in the TPC
1268 // track is the track you want to propagate
1269 // bfield is the magnetic field of your event
1270 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1271
1272 // Initialize the array to something indicating there was no propagation
1273 for(Int_t i=0;i<9;i++){
1274 for(Int_t j=0;j<3;j++){
1275 globalPositionsAtRadii[i][j]=-9999.;
1276 }
1277 }
1278
1279 // Make a copy of the track to not change parameters of the track
1280 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1281 //printf("\nAfter CopyFromVTrack\n");
1282 //etp.Print();
1283
1284 // The global position of the the track
1285 Double_t xyz[3]={-9999.,-9999.,-9999.};
1286
1287 // Counter for which radius we want
1288 Int_t iR=0;
1289 // The radii at which we get the global positions
1290 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1291 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1292 // The global radius we are at
1293 Float_t globalRadius=0;
1294
1295 // Propagation is done in local x of the track
1296 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1297 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1298 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1299 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1300 // adds to the global radius
1301
1302 // Stop if the propagation was not succesful. This can happen for low pt tracks
1303 // that don't reach outer radii
1304 if(!etp.PropagateTo(x,bfield))break;
1305 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1306 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1307
1308 // Roughly reached the radius we want
1309 if(globalRadius > Rwanted[iR]){
1310
1311 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1312 while (globalRadius>Rwanted[iR]){
1313 x-=.1;
1314 // printf("propagating to x %5.2f\n",x);
1315 if(!etp.PropagateTo(x,bfield))break;
1316 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1317 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1318 }
1319 //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]);
1320 globalPositionsAtRadii[iR][0]=xyz[0];
1321 globalPositionsAtRadii[iR][1]=xyz[1];
1322 globalPositionsAtRadii[iR][2]=xyz[2];
1323 //subtract primary vertex, "zero" track for correct mixed-event comparison
1324 globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1325 globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1326 globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1327
1328 // Indicate we want the next radius
1329 iR+=1;
1330 }
1331 if(iR>=8){
1332 // TPC edge reached
1333 return;
1334 }
1335 }
1336}
1337
1338bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1339 //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1340 //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1341
1342 bool newV0Wins = kFALSE;
1343 double pardiff[3] = {newPars[0]-oldPars[0],
1344 newPars[1]-oldPars[1],
1345 newPars[2]-oldPars[2]};
1346 if(cutChoice > 0 && cutChoice < 4){
1347 if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1348 }
1349 else if(cutChoice == 4){
1350 int newWinCount = 0;
1351 for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1352 if(newWinCount > 1) newV0Wins = kTRUE;
1353 }
1354 else{};
1355 return newV0Wins;
1356}
1357
1358bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1359{ // to flatten centrality distribution
1360 bool RejectEvent = kFALSE;
1361 int weightBinSign;
1362 if(MagField > 0) weightBinSign = 0;
1363 else weightBinSign = 1;
1364 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1365 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1366 int weightBinCent = (int) CentPercent;
1367 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1368
1369 return RejectEvent;
1370}