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