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