]> 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 code update (Matt Steinpreis)
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / K0Analysis / AliFemtoK0Analysis.cxx
CommitLineData
41dfc4d3
DRG		 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////////////////////////////////
17//
18// This class is used to perform femtoscopic analysis on K0s particles,
19// which are reconstructed using the AliAODv0 class.
20//
21// authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)
22//
23// Change log:
24// - TOF mismatch function calls changed (4/18/13)
25// - added minimum decay length cut (rarely used though) (3/28/13)
26// - K0 multiplicity histogram now filled with "unskippedCount" instead
27// of k0Count (which included skipped k0s with shared daughters)
28// (3/25/13)
29// - added hists for 3D mom. in LF and PRF (3/28/13)
30// - changed calling of PIDResponse (should be same actions) (3/28/13)
31// - keep "side" K0s for mass plot (4/18)
32// - tweaked loading and skipping appropriately
33// - use merit test to skip sides (against good and side K0s)
34// - a good K0 cant be skipped by a side
35// - moved TPC propagation (via Hans' method) up to v0 level, which now
36// uses an AliAODTrack(AliVTrack) instead of AliESDtrack (5/31/13)
37// - added primary vertex subtraction in TPC propagation (5/31/13)
38// - removed all instances of AliESDtrack usage (5/31/13)
39// - removed old separation method/histograms (5/31/13)
40// - tidied up LCMS boost (6/10/13)
41// - added new boosting prescription, get q out-side-long for LCMS and PRF (6/24/13)
42// - added histograms and values for LCMS momenta (for simulation)
43// - added random particle order switch in correlations (9/09/13)
44// - added more bins for 3D OSL analysis (9/19/13)
45// - added merit cut choice, pass as argument (10/16/13)
46// - 1-mass, 2-v0dca, 3-dddca, 4-combination (used to be v0dca)
47// - added passable argument for two-track minimum separation (10/16/13)
48// - added boolean to turn off field-sign dependence for train (10/30/13)
49// - changed destructors to minimize lost memory (11/27/13)
50// - added Case3D to switch off all 3D objects (11/27/13)
51// - added centrality flattening routine (and switch) (12/04/13)
52// - added event plane stuff (12/11/13)
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
d467271e 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
d467271e 462 TH3F *fHist3DOSLCutsSignal[3][5][3]; //3 cent bins, 5 parameters, 3 cut values
463 TH3F *fHist3DOSLCutsBkg[3][5][3];
b0a793ac 464 if(fCutCheck){
d467271e 465 for(int i3D=0;i3D<3;i3D++){
b0a793ac 466 for(int j3D=0;j3D<5;j3D++){
d467271e 467 for(int k3D=0;k3D<3;k3D++){
468 TString *histname = new TString("fHist3DOSLCuts");
469 *histname += i3D;
470 *histname += j3D;
471 *histname += k3D;
472 histname->Append("Signal");
473 fHist3DOSLCutsSignal[i3D][j3D][k3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
474 fOutputList->Add(fHist3DOSLCutsSignal[i3D][j3D][k3D]);
475 histname->Replace(17,6,"Bkg");
476 cout << histname->Data() << endl;
477 fHist3DOSLCutsBkg[i3D][j3D][k3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
478 fOutputList->Add(fHist3DOSLCutsBkg[i3D][j3D][k3D]);
479 }
b0a793ac 480 }
481 }
482 }
483
41dfc4d3
DRG		 484 //3D Spherical Harmonics
485 //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
486 //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
487 //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
488 //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
489 //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
490 //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
491 //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
492 //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
493 //fOutputList->Add(fHistSHCentLowKt);
494 //fOutputList->Add(fHistSHCentHighKt);
495 //fOutputList->Add(fHistSHSemiCentLowKt);
496 //fOutputList->Add(fHistSHSemiCentHighKt);
497 //fOutputList->Add(fHistSHCentLowKtBkg);
498 //fOutputList->Add(fHistSHCentHighKtBkg);
499 //fOutputList->Add(fHistSHSemiCentLowKtBkg);
500 //fOutputList->Add(fHistSHSemiCentHighKtBkg);
501
502 //side-side
503 //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
504 //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
505 //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
506 //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
507 //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
508 //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
509 //fOutputList->Add(fHistLeftLeftSignal);
510 //fOutputList->Add(fHistLeftRightSignal);
511 //fOutputList->Add(fHistRightRightSignal);
512 //fOutputList->Add(fHistLeftLeftBkg);
513 //fOutputList->Add(fHistLeftRightBkg);
514 //fOutputList->Add(fHistRightRightBkg);
515
516 //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
517 //fOutputList->Add(fHistSplitK0Sides);
518 //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
519 //fOutputList->Add(fHistSplitK0Centers);
520 //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
521 //fOutputList->Add(fHistQinvSignalNoSplit);
522
523 PostData(1, fOutputList);
524
525}
526
527//________________________________________________________________________
528void AliFemtoK0Analysis::Exec(Option_t *)
529{
530 // Main loop
531 // Called for each event
532 //cout<<"=========== Event # "<<fEventCount+1<<" ==========="<<endl;
533 fEventCount++;
534 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
535 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
536
537 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
538 bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
539 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
540 if(!isSelected) {
541 //cout << "Failed trigger selection." << endl;
542 return;
543 }
544
545 ///////////////////////////////////////////////////////////
546
547 unsigned int statusPos=0;
548 unsigned int statusNeg=0;
549
550 float bField=0;
551 bField = fAOD->GetMagneticField();
552 if(bField == 0) return;
553 if(fSignDep){
554 if(fFieldPos && bField < 0) return;
555 if(!fFieldPos && bField > 0) return;
556 }
557
558
41dfc4d3
DRG		 559 /////////////////////////////////////////////////
560
41dfc4d3
DRG		 561 //Centrality selection
562
563 AliCentrality *centrality = fAOD->GetCentrality();
564 float percent = centrality->GetCentralityPercentile("V0M");
565 int centBin=0;
566 //Printf("Centrality percent = %f", percent);
567
568 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
569 float multV0A=aodV0->GetMTotV0A();
570 float multV0C=aodV0->GetMTotV0C();
571
572 if(percent < 0) {
573 //Printf("No centrality info");
574 return;
575 }
576 if(percent < 0.1 && (multV0A + multV0C < 19500)){
577 //Printf("No centrality info");
578 return;
579 }
580 else if(percent <= 5) centBin=15;
581 else if(percent <= 10) centBin=14;
582 else if(percent <= 15) centBin=13;
583 else if(percent <= 20) centBin=12;
584 else if(percent <= 25) centBin=11;
585 else if(percent <= 30) centBin=10;
586 else if(percent <= 35) centBin=9;
587 else if(percent <= 40) centBin=8;
588 else if(percent <= 45) centBin=7;
589 else if(percent <= 50) centBin=6;
590 else if(percent <= 55) centBin=5;
591 else if(percent <= 60) centBin=4;
592 else if(percent <= 65) centBin=3;
593 else if(percent <= 70) centBin=2;
594 else if(percent <= 75) centBin=1;
595 else if(percent <= 80) centBin=0;
596 else {
597 //Printf("Skipping Peripheral Event");
598 return;
599 }
600 if(percent > 10 && isCentral) return;
d467271e 601 if(fCutCheck && percent > 50) return; //only looking at 0-50% for Cut Check
41dfc4d3
DRG		 602 ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
603
604 //flatten centrality dist.
605 if(percent < 9){
606 if(fFlatCent){
607 if(RejectEventCentFlat(bField,percent)) return;
608 }
609 }
610 ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
611
612 //Vertexing
613 AliAODVertex *primaryVertex;
614 double vertex[3]={0};
615 primaryVertex = fAOD->GetPrimaryVertex();
616 vertex[0]=primaryVertex->GetX();
617 vertex[1]=primaryVertex->GetY();
618 vertex[2]=primaryVertex->GetZ();
619 if(vertex[0]<10e-5 && vertex[1]<10e-5 && vertex[2]<10e-5) return;
620 if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
621
18d901d5
DRG		 622 int zBin=0;
623 double zStep=2*10/double(kZVertexBins), zStart=-10.;
41dfc4d3
DRG		 624 for(int i=0; i<kZVertexBins; i++)
625 {
18d901d5
DRG		 626 if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
627 {
628 zBin=i;
629 break;
630 }
41dfc4d3
DRG		 631 }
632
633 //Event plane
634 int psiBin = 0;
635 AliEventplane *eventplane = fAOD->GetEventplane();
636 if(fPsiBinning && !eventplane) return;
637 double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
41dfc4d3
DRG		 638 ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
639
0e22b4bd 640 double psiStep = PI/double(fNPsiBins);
18d901d5 641 double psiStart = -0.5*PI;
0e22b4bd 642 for(int i=0; i<fNPsiBins; i++)
18d901d5
DRG		 643 {
644 if((psiEP > psiStart+i*psiStep) && (psiEP < psiStart+(i+1)*psiStep))
645 {
646 psiBin = i;
647 break;
648 }
649 }
650 if(!fPsiBinning) psiBin = 0;
651
41dfc4d3
DRG		 652////////////////////////////////////////////////////////////////
653//Cut Values and constants
654
655 //const bool kMCCase = kFALSE; //switch for MC analysis
656 const int kMaxNumK0 = 300; //maximum number of K0s, array size
657 const float kMinDCAPrimaryPion = 0.4; //minimum dca of pions to primary
658 const float kMaxDCADaughtersK0 = 0.3; //maximum dca of pions to each other - 3D
659 const float kMaxDCAK0 = 0.3; //maximum dca of K0 to primary
660 const float kMaxDLK0 = 30.0; //maximum decay length of K0
661 const float kMinDLK0 = fMinDecayLength; //minimum decay length of K0
662 const float kEtaCut = 0.8; //maximum |pseudorapidity|
663 const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle
664
665 const float kMinSeparation = fMinSep; //minimum daughter (pair) separation
666
667 const float kTOFLow = 0.8; //boundary for TOF usage
668 const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas
669 const float kMaxTPCSigmaPion = 3.0; //TPC # of sigmas
670
671 //const float kMassPion = .13957;
672 const float kMassK0Short = .497614; //true PDG masses
673
b0a793ac 674 //for cut checks
d467271e 675 double kCheckMassLow [3] = {0.49,0.48,0.45};
676 double kCheckMassHigh [3] = {0.505,0.515,0.550};
677 double kCheckDCAK0 [3] = {0.1,0.3,1.0};
678 double kCheckDCAPi [3] = {1.0,0.4,0.1};
679 double kCheckDCAPiPi [3] = {0.1,0.3,1.0};
680 double kCheckAvgSep [3] = {10.0,5.0,0.0};
b0a793ac 681
41dfc4d3
DRG		 682////////////////////////////////////////////////////////////////
683 //v0 tester
684////////////////////////////////////////////////////////////////
685 int v0Count = 0; //number of v0s (entries in array)
686 int k0Count = 0; //number of good K0s
687
688 AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
689
690 //for daughter sharing studies
691 //int idArray[100] = {0};
692 //int idCount = 0;
693
694 //for MC
695 //TClonesArray *mcArray = 0x0;
696 //if(kMCCase){
697 //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
698 //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
699
700 for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
701 {
702 bool goodK0 = kFALSE;
703 bool goodPiPlus = kFALSE;
704 bool goodPiMinus = kFALSE;
705
706 //load v0 track
707 AliAODv0* v0 = fAOD->GetV0(i);
708 if(!v0) continue;
709 if(fOnlineCase){
710 if(!(v0->GetOnFlyStatus())) continue;
711 } //for online
712 else{
713 if((v0->GetOnFlyStatus())) continue; //for offline
714 }
715
716 //for on-the-fly ordering
717 AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
718 short int pos0or1;
719 short int neg0or1;
720 bool orderswitch = kFALSE;
721 if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
722 else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
723
724 //load daughter tracks
725 AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
726 AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
727 if(!prongTrackPos) continue;
728 if(!prongTrackNeg) continue;
729
730 //daughter cuts
731 if(v0->PtProng(pos0or1) < .15) continue;
732 if(v0->PtProng(neg0or1) < .15) continue;
733 if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
734 if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
735
736 //load status for PID
737 statusPos=prongTrackPos->GetStatus();
738 if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
739 prongTrackPos->SetAODEvent(fAOD);
740 statusNeg=prongTrackNeg->GetStatus();
741 if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
742 prongTrackNeg->SetAODEvent(fAOD);
743
744 //TPC PID
745 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
746 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
747
748 //Positive daughter identification TOF
749 float probMis;
750 AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
751 double Ppos = v0->PProng(pos0or1);
752 if(Ppos > kTOFLow) //PiPlus
753 {
754 //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
755 if(AliPIDResponse::kDetPidOk == statusPosTOF)
756 {
757 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
758 if(probMis < 0.01) //avoid TOF-TPC mismatch
759 {
760 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
761 else goodPiPlus = kFALSE;
762 }
763 }
764 }
765 //Negative daughter identification TOF
766 AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
767 double Pneg = v0->PProng(neg0or1);
768 if(Pneg > kTOFLow) //PiMinus
769 {
770 //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
771 if(AliPIDResponse::kDetPidOk == statusNegTOF)
772 {
773 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
774 if(probMis < 0.01) //avoid TOF-TPC mismatch
775 {
776 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
777 else goodPiMinus = kFALSE;
778 }
779 }
780 }
781
782 //K0 cuts
b0a793ac 783 if(!goodPiMinus || !goodPiPlus) continue;
41dfc4d3
DRG		 784 if(v0->Eta() > kEtaCut) continue;
785 if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;
786 if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;
41dfc4d3
DRG		 787 if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;
788 if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;
41dfc4d3 789
b0a793ac 790 double v0Dca = v0->DcaV0ToPrimVertex();
791 if(!fCutCheck){
792 if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;
793 if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue;
794 if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;
795 if(v0Dca > kMaxDCAK0) continue;
796 }
797 else{
d467271e 798 if(v0->DcaNegToPrimVertex() < kCheckDCAPi[2]) continue;
799 if(v0->DcaPosToPrimVertex() < kCheckDCAPi[2]) continue;
800 if(v0->DcaV0Daughters() > kCheckDCAPiPi[2]) continue;
801 if(v0Dca > kCheckDCAK0[2]) continue;
b0a793ac 802 }
803
41dfc4d3
DRG		 804 //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
805
806 //for MC
807 //bool MCgood = kFALSE;
808 //if(kMCCase){
809 //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
810 //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
811 //if(mck0dp->GetMother() >= 0){
812 //if(mck0dp->GetMother() == mck0dn->GetMother()){
813 //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
814 //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
815 //if(abs(mck0->GetPdgCode()) == 310){
816 //MCgood = kTRUE;
817 //}
818 //}
819 //}
820 //}
821 //}// if kMCCase
822
b0a793ac 823 if(!fCutCheck){
824 if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
825 }
826 else{
d467271e 827 if(v0->MassK0Short() > kCheckMassLow[2] && v0->MassK0Short() < kCheckMassHigh[2]) goodK0 = kTRUE;
b0a793ac 828 }
829
41dfc4d3
DRG		 830 //Check for shared daughters, using v0 DCA to judge
831 bool v0JudgeNew; //true if new v0 beats old
832 tempK0[v0Count].fSkipShared = kFALSE;
833 double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
834 if(fMeritCase){
835 for(int iID = 0; iID<v0Count; iID++){
836 if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old
837 if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
838 double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
839 v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
840 if(!v0JudgeNew){ //if old beats new...
841 if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing...
842 else{ //but if bad old beats new bad, or good old beats anything, skip new
843 tempK0[v0Count].fSkipShared = kTRUE; //skip new
844 break; //no need to keep checking others
845 }
846 }
847 else{ //if new beats old...
848 if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...
849 else{ //but if bad new beats bad old, or good new beats anything, skip old
850 tempK0[iID].fSkipShared = kTRUE; //skip old
851 if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
852 }
853 }
854 }
855 }
856 }
857 if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0
b0a793ac 858 }//if MeritCase
859
860 //for cut check
861 if(fCutCheck){
d467271e 862 for(int iSet=0;iSet<4;iSet++){ //number of cut pars (not counting AvgSep)
863 for(int jSet=0;jSet<3;jSet++){ //number of cut values
864 tempK0[v0Count].fCutPass[iSet][jSet] = kFALSE;
865 }
866 }
867 for(int jCut = 0;jCut<3;jCut++){
868 if(v0->MassK0Short() > kCheckMassLow[jCut] && v0->MassK0Short() < kCheckMassHigh[jCut])
b0a793ac 869 tempK0[v0Count].fCutPass[0][jCut] = kTRUE;
870 if(v0Dca < kCheckDCAK0[jCut]) tempK0[v0Count].fCutPass[1][jCut] = kTRUE;
871 if(v0->DcaPosToPrimVertex() > kCheckDCAPi[jCut] && v0->DcaNegToPrimVertex() > kCheckDCAPi[jCut])
872 tempK0[v0Count].fCutPass[2][jCut] = kTRUE;
873 if(v0->DcaV0Daughters() < kCheckDCAPiPi[jCut]) tempK0[v0Count].fCutPass[3][jCut] = kTRUE;
874 }
875 }
876
41dfc4d3
DRG		 877 //load parameters into temporary class instance
878 if(v0Count < kMaxNumK0)
879 {
880 if(goodK0){
881 tempK0[v0Count].fK0 = kTRUE;
882 k0Count++;
883 }
884 else tempK0[v0Count].fK0 = kFALSE;
885
886 //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
887 //else tempK0[v0Count].fSideLeft = kFALSE;
888 //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
889 //else tempK0[v0Count].fSideRight = kFALSE;
890 //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
891
892 tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();
893 tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();
894 tempK0[v0Count].fMomentum[0] = v0->Px();
895 tempK0[v0Count].fMomentum[1] = v0->Py();
896 tempK0[v0Count].fMomentum[2] = v0->Pz();
897 tempK0[v0Count].fPt = v0->Pt();
898 tempK0[v0Count].fMass = v0->MassK0Short();
899 tempK0[v0Count].fV0Dca = v0Dca;
900
901 //for hists
902 tempK0[v0Count].fDDDca = v0->DcaV0Daughters();
903 tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);
904 tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);
905 tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);
906 tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);
907 tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);
908 if(!orderswitch){
909 tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();
910 tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();
911 }
912 else{
913 tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();
914 tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();
915 }
916
917 //for psi studies
918 double v0Phi = v0->Phi(); //between [0,2pi]
919 double v0PhiPsi = v0Phi-psiEP;
920 if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
921 else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
922 else{};
95a1cff1 923 tempK0[v0Count].fPhi = v0Phi;
41dfc4d3
DRG		 924 tempK0[v0Count].fPhiPsi = v0PhiPsi;
925
926 //for separation
927 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
928 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
929 //for DPC
930 prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
931 prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
932
933 //if(idCount < 50){
934 // if(goodK0){
935 // idArray[idCount*2] = prongTrackPos->GetID();
936 // idArray[idCount*2+1] = prongTrackNeg->GetID();
937 // idCount++;
938 //}}
939
940 v0Count++;
941 }
942
943 }//v0
41dfc4d3
DRG		 944 if(k0Count<2) return; //only keep events with more than 1 good K0
945
946 //Add Event to buffer - this is for event mixing
947 fEC[zBin][centBin][psiBin]->FIFOShift();
948 (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
949 (fEvt)->fFillStatus = 1;
950 int unskippedCount = 0;
951 for(int i=0;i<v0Count;i++)
952 {
953 if(!tempK0[i].fSkipShared) //don't include skipped v0s (from shared daughters)
954 {
955 ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
b0a793ac 956
d467271e 957 //if(fCutCheck){
958 // for(int iCut=1;iCut<4;iCut++){
959 // for(int jCut=0;jCut<5;jCut++){
960 // TString *histname = new TString("fHistMassCuts");
961 // *histname += iCut;
962 // *histname += jCut;
963 // if(tempK0[i].fCutPass[iCut][jCut]) ((TH1F*)fOutputList->FindObject(histname->Data()))->Fill(tempK0[i].fMass);
964 // }
965 // }
966 //}
b0a793ac 967
41dfc4d3
DRG		 968 if(tempK0[i].fK0) //make sure particle is good (mass)
969 {
970 (fEvt)->fK0Particle[unskippedCount] = tempK0[i]; //load good, unskipped K0s
971 unskippedCount++; //count good, unskipped K0s
972 }
973 }
974 }
975 (fEvt)->fNumV0s = unskippedCount;
976 //Printf("Number of v0s: %d", v0Count);
977 //Printf("Number of K0s: %d", k0Count);
978 delete [] tempK0; tempK0 = NULL;
979
980 ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount); // changed 3/25, used to be "k0Count"
981 ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
982
983 //Printf("Reconstruction Finished. Starting pair studies.");
984
985 //////////////////////////////////////////////////////////////////////
986 // Correlations
987 //////////////////////////////////////////////////////////////////////
988
989 float px1, py1, pz1, px2, py2, pz2; //single kaon values
990 float en1, en2; //single kaon values
991 //float pt1, pt2; //single kaon values
992 float pairPx, pairPy, pairPz, pairP0; //pair momentum values
993 float pairPt, pairMt, pairKt; //pair momentum values
994 float pairMInv, pairPDotQ;
995 float qinv, q0, qx, qy, qz; //pair q values
996 //float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values
997 float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF
998 //float qOutLCMS;
999 float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame
1000 float betasq, gamma;
1001 float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
1002 float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
1003
1004
1005 for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
1006 {
1007 //single particle histograms (done here to avoid "skipped" v0s
1008 ((TH1F*)fOutputList->FindObject("fHistDCADaughters")) ->Fill((fEvt)->fK0Particle[i].fDDDca);
1009 ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0")) ->Fill((fEvt)->fK0Particle[i].fDecayLength);
1010 ((TH1F*)fOutputList->FindObject("fHistDCAK0")) ->Fill((fEvt)->fK0Particle[i].fV0Dca);
1011 ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
1012 ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus")) ->Fill((fEvt)->fK0Particle[i].fPosDca);
1013 ((TH2F*)fOutputList->FindObject("fHistPtK0")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
1014 ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
1015 ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
1016 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fPosPhi);
1017 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fNegPhi);
1018
1019 ((TH1F*)fOutputList->FindObject("fHistPx")) ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
1020 ((TH1F*)fOutputList->FindObject("fHistPy")) ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
1021 ((TH1F*)fOutputList->FindObject("fHistPz")) ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
1022
95a1cff1 1023 ((TH2F*)fOutputList->FindObject("fHistPhi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhi);
41dfc4d3
DRG		 1024 ((TH2F*)fOutputList->FindObject("fHistPhiPsi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
1025
1026 for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
1027 {
1028 int startbin=0;
1029 if(evnum==0) startbin=i+1;
1030
1031 for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
1032 {
1033 if(evnum==0) // Get rid of shared tracks
1034 {
1035 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
1036 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
1037 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
1038 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
1039 }
1040
1041 px1 = (fEvt)->fK0Particle[i].fMomentum[0];
1042 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
1043 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
1044 //pt1 = (fEvt)->fK0Particle[i].fPt;
1045 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
1046 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
1047 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
1048 //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
1049 if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
1050 double tempvar;
1051 tempvar = px1; px1 = px2; px2 = tempvar;
1052 tempvar = py1; py1 = py2; py2 = tempvar;
1053 tempvar = pz1; pz1 = pz2; pz2 = tempvar;
1054 }
1055
1056 en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
1057 en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
1058
1059 q0 = en1 - en2;
1060 qx = px1 - px2;
1061 qy = py1 - py2;
1062 qz = pz1 - pz2;
1063 qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
1064
1065 pairPx = px1 + px2;
1066 pairPy = py1 + py2;
1067 pairPz = pz1 + pz2;
1068 pairP0 = en1 + en2;
1069 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
1070 pairKt = pairPt/2.; //used for KT binning
1071 pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)
1072 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
1073 pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF
1074
1075 //PRF (this section will probably be removed in favor of later boosting section)
1076 p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length
1077 am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)
1078 epm = en1+en2+am12; //"energy plus mass"
1079 p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP
1080 if(am12 == 0) continue;
1081 h1 = (p112/epm - en1)/am12;
1082 ppx = px1+pairPx*h1; //px in PRF
1083 ppy = py1+pairPy*h1; //py in PRF
1084 ppz = pz1+pairPz*h1; //pz in PRF
1085 ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*
1086 ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
1087 ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
1088 ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
1089 ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
1090
1091 //relative momentum in out-side-long for LCMS and PRF
1092 if(pairMt == 0 || pairPt == 0) continue;
1093 qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
1094 qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
1095 //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
1096 qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
1097
1098 //finding gamma for gamma binning/hists (likely will be removed after tests)
1099 p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;
1100 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
1101 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
1102 p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;
1103 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
1104 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
1105 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
1106 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
1107 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
1108 gamma = 1./sqrt(1-betasq);
1109 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
1110 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
1111 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
1112 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
1113 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
1114 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1115 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1116 //getting bin numbers and names for 3D histogram
1117 TString *histname3D = new TString("fHist3DOSL");
1118 int ktBin;
1119 if(pairKt < 0.6) ktBin = 0;
1120 else if(pairKt < 0.8) ktBin = 1;
1121 else if(pairKt < 1.0) ktBin = 2;
1122 else ktBin = 3;
1123 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1124 *histname3D += ktBin;
1125
1126 //Spherical harmonics
1127 //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1128 //thetaSHCos = qLong/qLength;
1129 //thetaSH = acos(thetaSHCos);
1130 //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1131
1132 //Finding average separation of daughters throughout TPC - two-track cut
1133 float posPositions1[9][3] = {{0}};
1134 float negPositions1[9][3] = {{0}};
1135 float posPositions2[9][3] = {{0}};
1136 float negPositions2[9][3] = {{0}};
1137 for(int iPos = 0; iPos < 9; iPos++){
1138 for(int jPos = 0; jPos < 3; jPos++){
1139 posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1140 negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1141 posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1142 negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1143 }
1144 }
1145 float pMean = 0.; //average separation for positive daughters
1146 float nMean = 0.; //average separation for negative daughters
1147 float pDiff;
1148 float nDiff;
1149 float pMin = 9999.; //minimum separation (updates) - pos
1150 float nMin = 9999.; //minimum separation (updates) - neg
1151 double pCount=0; //counter for number of points used - pos
1152 double nCount=0; //counter for number of points used - neg
1153 for(int ss=0;ss<9;ss++){
1154 if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1155 pCount++;
1156 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));
1157 pMean = pMean + pDiff;
1158 if(pDiff < pMin) pMin = pDiff;
1159 }
1160 if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1161 nCount++;
1162 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));
1163 nMean = nMean + nDiff;
1164 if(nDiff < nMin) nMin = nDiff;
1165 }
1166 }
1167 pMean = pMean/pCount;
1168 nMean = nMean/nCount;
1169
1170 if(evnum==0){
1171 ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1172 ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1173 ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1174 ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1175 }
1176 else{
1177 ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1178 ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1179 ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1180 ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1181 }
1182
1183 //Decay plane coincidence
1184 //daughter momenta
1185 float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1186 float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1187 float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1188 float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1189 float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1190 float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1191 float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1192 float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1193 float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1194 float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1195 float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1196 float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1197
1198 float cross1[3];
1199 float cross2[3];
1200 cross1[0] = b1*f1-c1*e1;
1201 cross1[1] = c1*d1-a1*f1;
1202 cross1[2] = a1*e1-b1*d1;
1203 cross2[0] = b2*f2-c2*e2;
1204 cross2[1] = c2*d2-a2*f2;
1205 cross2[2] = a2*e2-b2*d2;
1206 float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1207 float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1208 float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1209
1210 if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1211 else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1212
d467271e 1213 bool SepPass[3] = {0};
b0a793ac 1214 if(!fCutCheck){
1215 if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1216 }
1217 else{
d467271e 1218 if(pMean < kCheckAvgSep[2] || nMean < kCheckAvgSep[2]) continue;
1219 for(int jCut=0;jCut<3;jCut++){
851c1532 1220 if(pMean > kCheckAvgSep[jCut] && nMean > kCheckAvgSep[jCut]) SepPass[jCut] = kTRUE;
b0a793ac 1221 }
1222 }
1223
41dfc4d3
DRG		 1224 //end separation studies
1225
1226 //Fill Histograms
1227 bool center1K0 = kFALSE; //accepted mass K0
1228 bool center2K0 = kFALSE;
1229 if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1230 if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1231 //bool CL1 = kFALSE;
1232 //bool CL2 = kFALSE;
1233 //bool CR1 = kFALSE;
1234 //bool CR2 = kFALSE;
1235 //if(center1K0 && center2K0){
1236 // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1237 // else CR1 = kTRUE;
1238 // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1239 // else CR2 = kTRUE;
1240 //}
1241
1242 //bool SideLeft1 = kFALSE;
1243 //bool SideLeft2 = kFALSE;
1244 //bool SideRight1 = kFALSE;
1245 //bool SideRight2 = kFALSE;
1246 //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1247 //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1248 //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1249 //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1250
1251 //for psi binning
1252 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1253 float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1254 bool inPlane1 = kFALSE;
1255 bool inPlane2 = kFALSE;
1256 if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1257 if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1258 if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1259 if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
95a1cff1 1260
1261 //for dphi, dphipsi
1262 float dPhi = fabs((fEvt)->fK0Particle[i].fPhi - (fEvt+evnum)->fK0Particle[j].fPhi);
1263 if(dPhi > PI) dPhi = 2*PI-dPhi;
1264 float dPhiPsi = fabs((fEvt)->fK0Particle[i].fPhiPsi - (fEvt+evnum)->fK0Particle[j].fPhiPsi);
1265 if(dPhiPsi > PI) dPhiPsi = 2*PI-dPhiPsi;
41dfc4d3 1266
d467271e 1267 int CutCentBin; //for cut check
1268 if(centBin > 13) CutCentBin = 0;
1269 else if(centBin > 9) CutCentBin = 1;
1270 else if(centBin > 5) CutCentBin = 2;
1271 else{};
41dfc4d3
DRG		 1272 if(evnum==0) //Same Event
1273 {
1274 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1275 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1276 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1277 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1278 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1279 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1280
1281 if(center1K0 && center2K0){
1282 //1D
1283 ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1284 //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1285 ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1286
1287 //eventplane
1288 if(inPlane1 && inPlane2)
1289 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1290 else if(!inPlane1 && !inPlane2)
1291 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
95a1cff1 1292
1293 //dPhi,dPhiPsi
1294 ((TH3F*)fOutputList->FindObject("fHistDPhi"))->Fill(centBin+1,pairKt,dPhi);
1295 ((TH3F*)fOutputList->FindObject("fHistDPhiPsi"))->Fill(centBin+1,pairKt,dPhiPsi);
1296
41dfc4d3
DRG		 1297
1298 //for mass bin study
1299 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1300 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1301 //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1302
1303 //3D
1304 if(fCase3D){
1305 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1306 histname3D->Append("Signal");
1307 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1308 }
b0a793ac 1309 }
1310
1311 //for cut check (3D, but fCase3D set to false)
1312 if(fCutCheck){
1313 for(int iCut=0;iCut<4;iCut++){//different cuts (4 + AvgSep)
1314 bool Skip = kFALSE;
d467271e 1315 for(int iCut2=0;iCut2<4;iCut2++){//for setting other cuts to usual value
b0a793ac 1316 if(iCut2 != iCut){
d467271e 1317 if(!(fEvt)->fK0Particle[i].fCutPass[iCut2][1] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut2][1]) Skip = kTRUE;
b0a793ac 1318 }
1319 }
d467271e 1320 if(!SepPass[1]) Skip = kTRUE; //set avg sep cut to usual value
b0a793ac 1321 if(Skip) continue;
d467271e 1322 for(int jCut=0;jCut<3;jCut++){//different cut values
b0a793ac 1323 TString *histname = new TString("fHist3DOSLCuts");
d467271e 1324 *histname += CutCentBin;
b0a793ac 1325 *histname += iCut;
1326 *histname += jCut;
1327 histname->Append("Signal");
1328 if((fEvt)->fK0Particle[i].fCutPass[iCut][jCut] && (fEvt+evnum)->fK0Particle[j].fCutPass[iCut][jCut])
1329 ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1330 }//jcut
1331 }//icut
1332
d467271e 1333 //for avg sep cutcheck
b0a793ac 1334 bool asSkip = kFALSE;
a7ac9250 1335 for(int iCut=0;iCut<4;iCut++){ //other parameters
d467271e 1336 if(!(fEvt)->fK0Particle[i].fCutPass[iCut][1] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut][1]) asSkip=kTRUE; //set other cuts to usual values
b0a793ac 1337 }
1338 if(asSkip) continue;
d467271e 1339 for(int jCut=0;jCut<3;jCut++){
b0a793ac 1340 TString *histname = new TString("fHist3DOSLCuts");
d467271e 1341 *histname += CutCentBin;
b0a793ac 1342 *histname += 4; //4 for AvgSep
1343 *histname += jCut;
1344 histname->Append("Signal");
1345 if(SepPass[jCut]) ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1346 }
1347 }//cutCheck
1348
1349
41dfc4d3
DRG		 1350 /*if(pairKt < 1.0){
1351 if(centBin > 13){
1352 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1353 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1354 else if(centBin > 9){
1355 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1356 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1357 else if(pairKt < 2.0){
1358 if(centBin > 13){
1359 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1360 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1361 else if(centBin > 9){
1362 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1363
1364 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1365
1366 }//centercenter
1367
1368 //side-side correlations
1369 //if(!splitK0sides){
1370 // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1371 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1372 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1373 //}
1374
1375 }//same event
1376
1377 else //Mixed Events
1378 {
1379 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1380 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1381 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1382 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1383
1384 if(center1K0 && center2K0){
1385 //1D
1386 ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1387
1388 //eventplane
1389 if(inPlane1 && inPlane2)
1390 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1391 else if(!inPlane1 && !inPlane2)
1392 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1393
95a1cff1 1394 ((TH3F*)fOutputList->FindObject("fHistDPhiBkg"))->Fill(centBin+1,pairKt,dPhi);
1395 ((TH3F*)fOutputList->FindObject("fHistDPhiPsiBkg"))->Fill(centBin+1,pairKt,dPhiPsi);
1396
41dfc4d3
DRG		 1397 //for mass bin study
1398 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1399 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1400 //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1401
1402 //3D
1403 if(fCase3D){
1404 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1405 histname3D->Replace(12,6,"Bkg");
1406 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1407 }
1408 }
b0a793ac 1409
d467271e 1410 //for cut check (3D, but fCase3D set to false)
b0a793ac 1411 if(fCutCheck){
1412 for(int iCut=0;iCut<4;iCut++){//different cuts (4 + AvgSep)
1413 bool Skip = kFALSE;
d467271e 1414 for(int iCut2=0;iCut2<4;iCut2++){//for setting other cuts to usual value
b0a793ac 1415 if(iCut2 != iCut){
d467271e 1416 if(!(fEvt)->fK0Particle[i].fCutPass[iCut2][1] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut2][1]) Skip = kTRUE;
b0a793ac 1417 }
1418 }
d467271e 1419 if(!SepPass[1]) Skip = kTRUE; //set avg sep cut to usual value
b0a793ac 1420 if(Skip) continue;
d467271e 1421 for(int jCut=0;jCut<3;jCut++){//different cut values
b0a793ac 1422 TString *histname = new TString("fHist3DOSLCuts");
d467271e 1423 *histname += CutCentBin;
b0a793ac 1424 *histname += iCut;
1425 *histname += jCut;
1426 histname->Append("Bkg");
1427 if((fEvt)->fK0Particle[i].fCutPass[iCut][jCut] && (fEvt+evnum)->fK0Particle[j].fCutPass[iCut][jCut])
1428 ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1429 }//jcut
1430 }//icut
1431
d467271e 1432 //for avg sep cutcheck
b0a793ac 1433 bool asSkip = kFALSE;
d467271e 1434 for(int iCut=0;iCut<4;iCut++){ //other parameters
1435 if(!(fEvt)->fK0Particle[i].fCutPass[iCut][1] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut][1]) asSkip=kTRUE; //set other cuts to usual values
b0a793ac 1436 }
1437 if(asSkip) continue;
d467271e 1438 for(int jCut=0;jCut<3;jCut++){
b0a793ac 1439 TString *histname = new TString("fHist3DOSLCuts");
d467271e 1440 *histname += CutCentBin;
b0a793ac 1441 *histname += 4; //4 for AvgSep
1442 *histname += jCut;
1443 histname->Append("Bkg");
1444 if(SepPass[jCut]) ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1445 }
d467271e 1446 }//cutCheck
b0a793ac 1447
1448
41dfc4d3
DRG		 1449 /*if(pairKt < 1.0){
1450 if(centBin > 13){
1451 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1452 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1453 else if(centBin > 9){
1454 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1455 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1456 else if(pairKt < 2.0){
1457 if(centBin > 13){
1458 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1459 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1460 else if(centBin > 9){
1461 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1462 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1463 }
1464
1465 //side-side correlations
1466 //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1467 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1468 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1469
1470 }//Mixed Events
1471
1472 }//past event
1473 }//event buffer
1474 }//current event
1475
1476 // Post output data.
1477 PostData(1, fOutputList);
1478
1479 }
1480//________________________________________________________________________
1481void AliFemtoK0Analysis::Terminate(Option_t *)
1482{
1483 // Called once at the end of the query
1484 cout<<"Done"<<endl;
1485
1486}
1487
1488//_________________________________________________________________________
1489void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1490 // Gets the global position of the track at nine different radii in the TPC
1491 // track is the track you want to propagate
1492 // bfield is the magnetic field of your event
1493 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1494
1495 // Initialize the array to something indicating there was no propagation
1496 for(Int_t i=0;i<9;i++){
1497 for(Int_t j=0;j<3;j++){
1498 globalPositionsAtRadii[i][j]=-9999.;
1499 }
1500 }
1501
1502 // Make a copy of the track to not change parameters of the track
1503 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1504 //printf("\nAfter CopyFromVTrack\n");
1505 //etp.Print();
1506
1507 // The global position of the the track
1508 Double_t xyz[3]={-9999.,-9999.,-9999.};
1509
1510 // Counter for which radius we want
1511 Int_t iR=0;
1512 // The radii at which we get the global positions
1513 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1514 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1515 // The global radius we are at
1516 Float_t globalRadius=0;
1517
1518 // Propagation is done in local x of the track
1519 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1520 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1521 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1522 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1523 // adds to the global radius
1524
1525 // Stop if the propagation was not succesful. This can happen for low pt tracks
1526 // that don't reach outer radii
1527 if(!etp.PropagateTo(x,bfield))break;
1528 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1529 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1530
1531 // Roughly reached the radius we want
1532 if(globalRadius > Rwanted[iR]){
1533
1534 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1535 while (globalRadius>Rwanted[iR]){
1536 x-=.1;
1537 // printf("propagating to x %5.2f\n",x);
1538 if(!etp.PropagateTo(x,bfield))break;
1539 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1540 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1541 }
1542 //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]);
1543 globalPositionsAtRadii[iR][0]=xyz[0];
1544 globalPositionsAtRadii[iR][1]=xyz[1];
1545 globalPositionsAtRadii[iR][2]=xyz[2];
1546 //subtract primary vertex, "zero" track for correct mixed-event comparison
1547 globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1548 globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1549 globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1550
1551 // Indicate we want the next radius
1552 iR+=1;
1553 }
1554 if(iR>=8){
1555 // TPC edge reached
1556 return;
1557 }
1558 }
1559}
1560
1561bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1562 //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1563 //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1564
1565 bool newV0Wins = kFALSE;
1566 double pardiff[3] = {newPars[0]-oldPars[0],
1567 newPars[1]-oldPars[1],
1568 newPars[2]-oldPars[2]};
1569 if(cutChoice > 0 && cutChoice < 4){
1570 if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1571 }
1572 else if(cutChoice == 4){
1573 int newWinCount = 0;
1574 for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1575 if(newWinCount > 1) newV0Wins = kTRUE;
1576 }
1577 else{};
1578 return newV0Wins;
1579}
1580
1581bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1582{ // to flatten centrality distribution
1583 bool RejectEvent = kFALSE;
1584 int weightBinSign;
1585 if(MagField > 0) weightBinSign = 0;
1586 else weightBinSign = 1;
1587 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1588 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1589 int weightBinCent = (int) CentPercent;
1590 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1591
1592 return RejectEvent;
1593}