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