1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
16 ////////////////////////////////////////////////////////////////////////////
\r
18 // This class is used to perform femtoscopic analysis on K0s particles,
\r
19 // which are reconstructed using the AliAODv0 class.
\r
21 // authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)
\r
24 // - TOF mismatch function calls changed (4/18/13)
\r
25 // - added minimum decay length cut (rarely used though) (3/28/13)
\r
26 // - K0 multiplicity histogram now filled with "unskippedCount" instead
\r
27 // of k0Count (which included skipped k0s with shared daughters)
\r
29 // - added hists for 3D mom. in LF and PRF (3/28/13)
\r
30 // - changed calling of PIDResponse (should be same actions) (3/28/13)
\r
31 // - keep "side" K0s for mass plot (4/18)
\r
32 // - tweaked loading and skipping appropriately
\r
33 // - use merit test to skip sides (against good and side K0s)
\r
34 // - a good K0 cant be skipped by a side
\r
35 // - moved TPC propagation (via Hans' method) up to v0 level, which now
\r
36 // uses an AliAODTrack(AliVTrack) instead of AliESDtrack (5/31/13)
\r
37 // - added primary vertex subtraction in TPC propagation (5/31/13)
\r
38 // - removed all instances of AliESDtrack usage (5/31/13)
\r
39 // - removed old separation method/histograms (5/31/13)
\r
40 // - tidied up LCMS boost (6/10/13)
\r
41 // - added new boosting prescription, get q out-side-long for LCMS and PRF (6/24/13)
\r
42 // - added histograms and values for LCMS momenta (for simulation)
\r
43 // - added random particle order switch in correlations (9/09/13)
\r
44 // - added more bins for 3D OSL analysis (9/19/13)
\r
45 // - added merit cut choice, pass as argument (10/16/13)
\r
46 // - 1-mass, 2-v0dca, 3-dddca, 4-combination (used to be v0dca)
\r
47 // - added passable argument for two-track minimum separation (10/16/13)
\r
48 // - added boolean to turn off field-sign dependence for train (10/30/13)
\r
49 // - changed destructors to minimize lost memory (11/27/13)
\r
50 // - added Case3D to switch off all 3D objects (11/27/13)
\r
51 ////////////////////////////////////////////////////////////////////////////////
\r
61 #include "TObject.h"
\r
62 #include "TObjString.h"
\r
69 #include "TProfile.h"
\r
70 #include "TCanvas.h"
\r
71 #include "TRandom3.h"
\r
73 #include "AliAnalysisTask.h"
\r
74 #include "AliAnalysisManager.h"
\r
76 #include "AliAODEvent.h"
\r
77 #include "AliAODInputHandler.h"
\r
78 #include "AliAODMCParticle.h"
\r
79 #include "AliAODv0.h"
\r
80 #include "AliAODRecoDecay.h"
\r
81 #include "AliCentrality.h"
\r
83 #include "AliFemtoK0Analysis.h"
\r
85 #define PI 3.1415927
\r
88 // Author: Matt Steinpreis, adapted from Dhevan Gangadharan
\r
90 ClassImp(AliFemtoK0Analysis)
\r
92 //________________________________________________________________________
\r
93 AliFemtoK0Analysis::AliFemtoK0Analysis():
\r
94 AliAnalysisTaskSE(),
\r
100 fMinDecayLength(0.0),
\r
101 fMeritCutChoice(0),
\r
107 fRandomNumber(0x0),
\r
114 //________________________________________________________________________
\r
115 AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool SignDep, bool FieldPositive, bool OnlineCase, bool MeritCase, bool Case3D, float MinDL, int MeritCutChoice, float MinSep, bool FlatCent)
\r
116 : AliAnalysisTaskSE(name),
\r
118 fFieldPos(FieldPositive),
\r
119 fOnlineCase(OnlineCase),
\r
120 fMeritCase(MeritCase),
\r
122 fMinDecayLength(MinDL),
\r
123 fMeritCutChoice(MeritCutChoice),
\r
125 fFlatCent(FlatCent),
\r
129 fRandomNumber(0x0),
\r
136 fSignDep = SignDep;
\r
137 fFieldPos = FieldPositive;
\r
138 fOnlineCase = OnlineCase;
\r
139 fMeritCase = MeritCase;
\r
141 fMinDecayLength = MinDL;
\r
142 fMeritCutChoice = MeritCutChoice;
\r
144 fFlatCent = FlatCent;
\r
146 // Define output slots here
\r
148 DefineOutput(1, TList::Class());
\r
151 //________________________________________________________________________
\r
152 AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
\r
153 : AliAnalysisTaskSE(obj.fName),
\r
154 fSignDep(obj.fSignDep),
\r
155 fFieldPos(obj.fFieldPos),
\r
156 fOnlineCase(obj.fOnlineCase),
\r
157 fMeritCase(obj.fMeritCase),
\r
158 fCase3D(obj.fCase3D),
\r
159 fMinDecayLength(obj.fMinDecayLength),
\r
160 fMeritCutChoice(obj.fMeritCutChoice),
\r
161 fMinSep(obj.fMinSep),
\r
162 fFlatCent(obj.fFlatCent),
\r
163 fEventCount(obj.fEventCount),
\r
166 fRandomNumber(obj.fRandomNumber),
\r
169 fOutputList(obj.fOutputList),
\r
170 fPidAOD(obj.fPidAOD)
\r
173 //________________________________________________________________________
\r
174 AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
\r
176 //Assignment operator
\r
177 if (this == &obj) return *this;
\r
179 fSignDep = obj.fSignDep;
\r
180 fFieldPos = obj.fFieldPos;
\r
181 fOnlineCase = obj.fOnlineCase;
\r
182 fMeritCase = obj.fMeritCase;
\r
183 fCase3D = obj.fCase3D;
\r
184 fMinDecayLength= obj.fMinDecayLength;
\r
185 fMeritCutChoice= obj.fMeritCutChoice;
\r
186 fMinSep = obj.fMinSep;
\r
187 fFlatCent = obj.fFlatCent;
\r
188 fEventCount = obj.fEventCount;
\r
191 fRandomNumber = obj.fRandomNumber;
\r
194 fOutputList = obj.fOutputList;
\r
195 fPidAOD = obj.fPidAOD;
\r
199 //________________________________________________________________________
\r
200 AliFemtoK0Analysis::~AliFemtoK0Analysis()
\r
203 for(unsigned short i=0; i<kZVertexBins; i++)
\r
205 for(unsigned short j=0; j<kCentBins; j++)
\r
207 fEC[i][j]->~AliFemtoK0EventCollection();
\r
210 delete[] fEC[i]; fEC[i] = NULL;
\r
212 delete[] fEC; fEC = NULL;
\r
214 if(fEC){ delete fEC; fEC = NULL;}
\r
215 if(fRandomNumber){ delete fRandomNumber; fRandomNumber = NULL;}
\r
216 if(fAOD){ delete fAOD; fAOD = NULL;}
\r
217 if(fOutputList){ delete fOutputList; fOutputList = NULL;}
\r
218 if(fPidAOD){ delete fPidAOD; fPidAOD = NULL;}
\r
220 //________________________________________________________________________
\r
221 void AliFemtoK0Analysis::MyInit()
\r
224 // One can set global variables here
\r
227 fEC = new AliFemtoK0EventCollection **[kZVertexBins];
\r
228 for(unsigned short i=0; i<kZVertexBins; i++)
\r
230 fEC[i] = new AliFemtoK0EventCollection *[kCentBins];
\r
232 for(unsigned short j=0; j<kCentBins; j++)
\r
234 fEC[i][j] = new AliFemtoK0EventCollection(kEventsToMix+1, kMultLimit);
\r
238 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
\r
239 fPidAOD = aodH->GetPIDResponse();
\r
241 fRandomNumber = new TRandom3(); //for 3D, random sign switching
\r
242 fRandomNumber->SetSeed(0);
\r
245 //________________________________________________________________________
\r
246 void AliFemtoK0Analysis::UserCreateOutputObjects()
\r
248 // Create histograms
\r
251 MyInit();// Initialize my settings
\r
253 fOutputList = new TList();
\r
254 fOutputList->SetOwner();
\r
256 TH1F *fHistCent = new TH1F("fHistCent","",100,0,100);
\r
257 fOutputList->Add(fHistCent);
\r
258 TH1F *fHistCentFlat = new TH1F("fHistCentFlat","",100,0,100);
\r
259 fOutputList->Add(fHistCentFlat);
\r
260 TH1F *fHistCentUsed = new TH1F("fHistCentUsed","",100,0,100);
\r
261 fOutputList->Add(fHistCentUsed);
\r
264 TH1F* fHistDCAPiPlus = new TH1F("fHistDCAPiPlus","",100,0,10);
\r
265 fOutputList->Add(fHistDCAPiPlus);
\r
266 TH1F* fHistDCAPiMinus = new TH1F("fHistDCAPiMinus","",100,0,10);
\r
267 fOutputList->Add(fHistDCAPiMinus);
\r
268 TH1F* fHistDCADaughters = new TH1F("fHistDCADaughters", "DCA of pions to each other", 50, 0., 0.5);
\r
269 fOutputList->Add(fHistDCADaughters);
\r
270 TH2F* fHistK0PiPlusPt = new TH2F("fHistK0PiPlusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
\r
271 fOutputList->Add(fHistK0PiPlusPt);
\r
272 TH2F* fHistK0PiMinusPt = new TH2F("fHistK0PiMinusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
\r
273 fOutputList->Add(fHistK0PiMinusPt);
\r
274 TH1F* fHistDaughterPhi = new TH1F("fHistDaughterPhi","",180,-PI,PI);
\r
275 fOutputList->Add(fHistDaughterPhi);
\r
278 TH1F* fHistMultK0 = new TH1F("fHistMultK0", "K0 multiplicity", 51, -0.5, 51-0.5);
\r
279 fOutputList->Add(fHistMultK0);
\r
280 TH2F* fHistPtK0 = new TH2F("fHistPtK0", "K0 pt distribution",kCentBins,.5,kCentBins+.5, 100, 0., 10.);
\r
281 fOutputList->Add(fHistPtK0);
\r
282 TH1F* fHistDecayLengthK0 = new TH1F("fHistDecayLengthK0", "K0 decay length", 100, 0., 100.);
\r
283 fOutputList->Add(fHistDecayLengthK0);
\r
284 TH1F* fHistDCAK0 = new TH1F("fHistDCAK0", "DCA of K0 to primary vertex", 40, 0., 0.4);
\r
285 fOutputList->Add(fHistDCAK0);
\r
286 TH2F* fHistKtK0 = new TH2F("fHistKtK0", "Kt distribution of K0 pairs", kCentBins, .5, kCentBins+.5, 300, 0., 3.);
\r
287 fOutputList->Add(fHistKtK0);
\r
289 TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);
\r
290 TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);
\r
291 TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);
\r
292 TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);
\r
293 TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);
\r
294 TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);
\r
295 TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);
\r
296 fOutputList->Add(fHistPx);
\r
297 fOutputList->Add(fHistPy);
\r
298 fOutputList->Add(fHistPz);
\r
299 fOutputList->Add(fHistPxCM);
\r
300 fOutputList->Add(fHistPyCM);
\r
301 fOutputList->Add(fHistPzCM);
\r
302 fOutputList->Add(fHistKsCM);
\r
304 TH1F* fHistPOutLCMS = new TH1F("fHistPOutLCMS","",200,0,2);
\r
305 TH1F* fHistPSideLCMS = new TH1F("fHistPSideLCMS","",200,0,2);
\r
306 TH1F* fHistPLongLCMS = new TH1F("fHistPLongLCMS","",200,0,2);
\r
307 fOutputList->Add(fHistPOutLCMS);
\r
308 fOutputList->Add(fHistPSideLCMS);
\r
309 fOutputList->Add(fHistPLongLCMS);
\r
311 //pair gamma (LCMS to PRF, OSL)
\r
312 TH2F* fHistGamma = new TH2F("fHistGamma","Gamma from LCMS to PRF",500,1,5,100,0,1);
\r
313 fOutputList->Add(fHistGamma);
\r
315 //invariant mass distributions
\r
316 TH3F* fHistMass = new TH3F("fHistMass","",kCentBins,.5,kCentBins+.5,50,0.,5.,400,.3,.7);
\r
317 fOutputList->Add(fHistMass);
\r
318 //TH3F* fHistMassPtCFK0 = new TH3F("fHistMassPtCFK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
\r
319 //fOutputList->Add(fHistMassPtCFK0);
\r
320 //TH3F* fHistMassPtCFBkgK0 = new TH3F("fHistMassPtCFBkgK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
\r
321 //fOutputList->Add(fHistMassPtCFBkgK0);
\r
322 //TH3F* fHistMassQKt = new TH3F("fHistMassQKt","",100,0,1,200,0,2,200,.4,.6);
\r
323 //fOutputList->Add(fHistMassQKt);
\r
324 //TH3F* fHistMassKtK0 = new TH3F("fHistMassKtK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
\r
325 //fOutputList->Add(fHistMassKtK0);
\r
326 //TH3F* fHistMassKtBkgK0 = new TH3F("fHistMassKtBkgK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
\r
327 //fOutputList->Add(fHistMassKtBkgK0);
\r
329 //separation studies
\r
330 TH1F* fHistSepNumPos = new TH1F("fHistSepNumPos","",200,0,20);
\r
331 fOutputList->Add(fHistSepNumPos);
\r
332 TH1F* fHistSepDenPos = new TH1F("fHistSepDenPos","",200,0,20);
\r
333 fOutputList->Add(fHistSepDenPos);
\r
334 TH1F* fHistSepNumNeg = new TH1F("fHistSepNumNeg","",200,0,20);
\r
335 fOutputList->Add(fHistSepNumNeg);
\r
336 TH1F* fHistSepDenNeg = new TH1F("fHistSepDenNeg","",200,0,20);
\r
337 fOutputList->Add(fHistSepDenNeg);
\r
339 TH2F* fHistSepNumPos2 = new TH2F("fHistSepNumPos2","",100,0,20,100,0,20);
\r
340 TH2F* fHistSepDenPos2 = new TH2F("fHistSepDenPos2","",100,0,20,100,0,20);
\r
341 TH2F* fHistSepNumNeg2 = new TH2F("fHistSepNumNeg2","",100,0,20,100,0,20);
\r
342 TH2F* fHistSepDenNeg2 = new TH2F("fHistSepDenNeg2","",100,0,20,100,0,20);
\r
343 fOutputList->Add(fHistSepNumPos2);
\r
344 fOutputList->Add(fHistSepDenPos2);
\r
345 fOutputList->Add(fHistSepNumNeg2);
\r
346 fOutputList->Add(fHistSepDenNeg2);
\r
349 TH2F* fHistSepDPC = new TH2F("fHistSepDPC","",200,-1,1,50,0,10);
\r
350 TH2F* fHistSepDPCBkg = new TH2F("fHistSepDPCBkg","",200,-1,1,50,0,10);
\r
351 fOutputList->Add(fHistSepDPC);
\r
352 fOutputList->Add(fHistSepDPCBkg);
\r
354 /////////Signal Distributions///////////////////
\r
357 TH3F* fHistQinvSignal = new TH3F("fHistQinvSignal","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
358 fOutputList->Add(fHistQinvSignal);
\r
359 TH3F* fHistQinvBkg = new TH3F("fHistQinvBkg","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
\r
360 fOutputList->Add(fHistQinvBkg);
\r
362 //mass bins within peak
\r
363 //TH3F* fHistCLCLSignal = new TH3F("fHistCLCLSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
364 //TH3F* fHistCLCLBkg = new TH3F("fHistCLCLBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
365 //TH3F* fHistCLCRSignal = new TH3F("fHistCLCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
366 //TH3F* fHistCLCRBkg = new TH3F("fHistCLCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
367 //TH3F* fHistCRCRSignal = new TH3F("fHistCRCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
368 //TH3F* fHistCRCRBkg = new TH3F("fHistCRCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
369 //fOutputList->Add(fHistCLCLSignal);
\r
370 //fOutputList->Add(fHistCLCLBkg);
\r
371 //fOutputList->Add(fHistCLCRSignal);
\r
372 //fOutputList->Add(fHistCLCRBkg);
\r
373 //fOutputList->Add(fHistCRCRSignal);
\r
374 //fOutputList->Add(fHistCRCRBkg);
\r
377 TH3F *fHist3DOSLSignal[10][4];
\r
378 TH3F *fHist3DOSLBkg[10][4];
\r
381 for(int i3D=0;i3D<10;i3D++){
\r
382 for(int j3D=0;j3D<4;j3D++){
\r
383 TString *histname = new TString("fHist3DOSL");
\r
386 histname->Append("Signal");
\r
387 fHist3DOSLSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
\r
388 fOutputList->Add(fHist3DOSLSignal[i3D][j3D]);
\r
389 histname->Replace(12,6,"Bkg");
\r
390 fHist3DOSLBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
\r
391 fOutputList->Add(fHist3DOSLBkg[i3D][j3D]);
\r
396 //3D Spherical Harmonics
\r
397 //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
398 //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
399 //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
400 //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
401 //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
402 //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
403 //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
404 //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
\r
405 //fOutputList->Add(fHistSHCentLowKt);
\r
406 //fOutputList->Add(fHistSHCentHighKt);
\r
407 //fOutputList->Add(fHistSHSemiCentLowKt);
\r
408 //fOutputList->Add(fHistSHSemiCentHighKt);
\r
409 //fOutputList->Add(fHistSHCentLowKtBkg);
\r
410 //fOutputList->Add(fHistSHCentHighKtBkg);
\r
411 //fOutputList->Add(fHistSHSemiCentLowKtBkg);
\r
412 //fOutputList->Add(fHistSHSemiCentHighKtBkg);
\r
415 //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
416 //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
417 //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
418 //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
419 //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
420 //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
421 //fOutputList->Add(fHistLeftLeftSignal);
\r
422 //fOutputList->Add(fHistLeftRightSignal);
\r
423 //fOutputList->Add(fHistRightRightSignal);
\r
424 //fOutputList->Add(fHistLeftLeftBkg);
\r
425 //fOutputList->Add(fHistLeftRightBkg);
\r
426 //fOutputList->Add(fHistRightRightBkg);
\r
428 //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
429 //fOutputList->Add(fHistSplitK0Sides);
\r
430 //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
431 //fOutputList->Add(fHistSplitK0Centers);
\r
432 //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
\r
433 //fOutputList->Add(fHistQinvSignalNoSplit);
\r
435 PostData(1, fOutputList);
\r
439 //________________________________________________________________________
\r
440 void AliFemtoK0Analysis::Exec(Option_t *)
\r
443 // Called for each event
\r
444 //cout<<"=========== Event # "<<fEventCount+1<<" ==========="<<endl;
\r
446 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
\r
447 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
\r
449 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
\r
450 bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
\r
451 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
\r
453 //cout << "Failed trigger selection." << endl;
\r
457 ///////////////////////////////////////////////////////////
\r
459 unsigned int statusPos=0;
\r
460 unsigned int statusNeg=0;
\r
463 bField = fAOD->GetMagneticField();
\r
464 if(bField == 0) return;
\r
466 if(fFieldPos && bField < 0) return;
\r
467 if(!fFieldPos && bField > 0) return;
\r
472 double zStep=2*10/double(kZVertexBins), zstart=-10.;
\r
474 /////////////////////////////////////////////////
\r
477 //Centrality selection
\r
479 AliCentrality *centrality = fAOD->GetCentrality();
\r
480 float percent = centrality->GetCentralityPercentile("V0M");
\r
482 //Printf("Centrality percent = %f", percent);
\r
484 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
\r
485 float multV0A=aodV0->GetMTotV0A();
\r
486 float multV0C=aodV0->GetMTotV0C();
\r
489 //Printf("No centrality info");
\r
492 if(percent < 0.1 && (multV0A + multV0C < 19500)){
\r
493 //Printf("No centrality info");
\r
496 else if(percent <= 5) centBin=15;
\r
497 else if(percent <= 10) centBin=14;
\r
498 else if(percent <= 15) centBin=13;
\r
499 else if(percent <= 20) centBin=12;
\r
500 else if(percent <= 25) centBin=11;
\r
501 else if(percent <= 30) centBin=10;
\r
502 else if(percent <= 35) centBin=9;
\r
503 else if(percent <= 40) centBin=8;
\r
504 else if(percent <= 45) centBin=7;
\r
505 else if(percent <= 50) centBin=6;
\r
506 else if(percent <= 55) centBin=5;
\r
507 else if(percent <= 60) centBin=4;
\r
508 else if(percent <= 65) centBin=3;
\r
509 else if(percent <= 70) centBin=2;
\r
510 else if(percent <= 75) centBin=1;
\r
511 else if(percent <= 80) centBin=0;
\r
513 //Printf("Skipping Peripheral Event");
\r
516 if(percent > 10 && isCentral) return;
\r
517 ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
\r
519 //flatten centrality dist.
\r
522 if(RejectEventCentFlat(bField,percent)) return;
\r
525 ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
\r
528 AliAODVertex *primaryVertex;
\r
529 double vertex[3]={0};
\r
530 primaryVertex = fAOD->GetPrimaryVertex();
\r
531 vertex[0]=primaryVertex->GetX();
\r
532 vertex[1]=primaryVertex->GetY();
\r
533 vertex[2]=primaryVertex->GetZ();
\r
534 if(vertex[0]<10e-5 && vertex[1]<10e-5 && vertex[2]<10e-5) return;
\r
535 if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
\r
537 for(int i=0; i<kZVertexBins; i++)
\r
539 if((vertex[2] > zstart+i*zStep) && (vertex[2] < zstart+(i+1)*zStep))
\r
546 ////////////////////////////////////////////////////////////////
\r
547 //Cut Values and constants
\r
549 //const bool kMCCase = kFALSE; //switch for MC analysis
\r
550 const int kMaxNumK0 = 300; //maximum number of K0s, array size
\r
551 const float kMinDCAPrimaryPion = 0.4; //minimum dca of pions to primary
\r
552 const float kMaxDCADaughtersK0 = 0.3; //maximum dca of pions to each other - 3D
\r
553 const float kMaxDCAK0 = 0.3; //maximum dca of K0 to primary
\r
554 const float kMaxDLK0 = 30.0; //maximum decay length of K0
\r
555 const float kMinDLK0 = fMinDecayLength; //minimum decay length of K0
\r
556 const float kEtaCut = 0.8; //maximum |pseudorapidity|
\r
557 const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle
\r
559 const float kMinSeparation = fMinSep; //minimum daughter (pair) separation
\r
561 const float kTOFLow = 0.8; //boundary for TOF usage
\r
562 const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas
\r
563 const float kMaxTPCSigmaPion = 3.0; //TPC # of sigmas
\r
565 //const float kMassPion = .13957;
\r
566 const float kMassK0Short = .497614; //true PDG masses
\r
568 ////////////////////////////////////////////////////////////////
\r
570 ////////////////////////////////////////////////////////////////
\r
571 int v0Count = 0; //number of v0s (entries in array)
\r
572 int k0Count = 0; //number of good K0s
\r
574 AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
\r
576 //for daughter sharing studies
\r
577 //int idArray[100] = {0};
\r
581 //TClonesArray *mcArray = 0x0;
\r
583 //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
\r
584 //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
\r
586 for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
\r
588 bool goodK0 = kFALSE;
\r
589 bool goodPiPlus = kFALSE;
\r
590 bool goodPiMinus = kFALSE;
\r
593 AliAODv0* v0 = fAOD->GetV0(i);
\r
596 if(!(v0->GetOnFlyStatus())) continue;
\r
599 if((v0->GetOnFlyStatus())) continue; //for offline
\r
602 //for on-the-fly ordering
\r
603 AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
\r
606 bool orderswitch = kFALSE;
\r
607 if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
\r
608 else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
\r
610 //load daughter tracks
\r
611 AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
\r
612 AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
\r
613 if(!prongTrackPos) continue;
\r
614 if(!prongTrackNeg) continue;
\r
617 if(v0->PtProng(pos0or1) < .15) continue;
\r
618 if(v0->PtProng(neg0or1) < .15) continue;
\r
619 if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
\r
620 if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
\r
622 //load status for PID
\r
623 statusPos=prongTrackPos->GetStatus();
\r
624 if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
\r
625 prongTrackPos->SetAODEvent(fAOD);
\r
626 statusNeg=prongTrackNeg->GetStatus();
\r
627 if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
\r
628 prongTrackNeg->SetAODEvent(fAOD);
\r
631 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
\r
632 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
\r
634 //Positive daughter identification TOF
\r
636 AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
\r
637 double Ppos = v0->PProng(pos0or1);
\r
638 if(Ppos > kTOFLow) //PiPlus
\r
640 //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
\r
641 if(AliPIDResponse::kDetPidOk == statusPosTOF)
\r
643 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
\r
644 if(probMis < 0.01) //avoid TOF-TPC mismatch
\r
646 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
\r
647 else goodPiPlus = kFALSE;
\r
651 //Negative daughter identification TOF
\r
652 AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
\r
653 double Pneg = v0->PProng(neg0or1);
\r
654 if(Pneg > kTOFLow) //PiMinus
\r
656 //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
\r
657 if(AliPIDResponse::kDetPidOk == statusNegTOF)
\r
659 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
\r
660 if(probMis < 0.01) //avoid TOF-TPC mismatch
\r
662 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
\r
663 else goodPiMinus = kFALSE;
\r
669 if(v0->Eta() > kEtaCut) continue;
\r
670 if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;
\r
671 if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;
\r
672 if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;
\r
673 if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue;
\r
674 if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;
\r
675 if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;
\r
676 if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;
\r
677 double v0Dca = v0->DcaV0ToPrimVertex();
\r
678 if(v0Dca > kMaxDCAK0) continue;
\r
679 if(!goodPiMinus || !goodPiPlus) continue;
\r
681 //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
\r
684 //bool MCgood = kFALSE;
\r
686 //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
\r
687 //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
\r
688 //if(mck0dp->GetMother() >= 0){
\r
689 //if(mck0dp->GetMother() == mck0dn->GetMother()){
\r
690 //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
\r
691 //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
\r
692 //if(abs(mck0->GetPdgCode()) == 310){
\r
700 if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
\r
701 //else continue; //removed, Apr 18
\r
703 //Check for shared daughters, using v0 DCA to judge
\r
704 bool v0JudgeNew; //true if new v0 beats old
\r
705 tempK0[v0Count].fSkipShared = kFALSE;
\r
706 double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
\r
708 for(int iID = 0; iID<v0Count; iID++){
\r
709 if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old
\r
710 if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
\r
711 double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
\r
712 v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
\r
713 if(!v0JudgeNew){ //if old beats new...
\r
714 if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing...
\r
715 else{ //but if bad old beats new bad, or good old beats anything, skip new
\r
716 tempK0[v0Count].fSkipShared = kTRUE; //skip new
\r
717 break; //no need to keep checking others
\r
720 else{ //if new beats old...
\r
721 if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...
\r
722 else{ //but if bad new beats bad old, or good new beats anything, skip old
\r
723 tempK0[iID].fSkipShared = kTRUE; //skip old
\r
724 if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
\r
730 if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0
\r
733 //load parameters into temporary class instance
\r
734 if(v0Count < kMaxNumK0)
\r
737 tempK0[v0Count].fK0 = kTRUE;
\r
740 else tempK0[v0Count].fK0 = kFALSE;
\r
742 //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
\r
743 //else tempK0[v0Count].fSideLeft = kFALSE;
\r
744 //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
\r
745 //else tempK0[v0Count].fSideRight = kFALSE;
\r
746 //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
\r
748 tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();
\r
749 tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();
\r
750 tempK0[v0Count].fMomentum[0] = v0->Px();
\r
751 tempK0[v0Count].fMomentum[1] = v0->Py();
\r
752 tempK0[v0Count].fMomentum[2] = v0->Pz();
\r
753 tempK0[v0Count].fPt = v0->Pt();
\r
754 tempK0[v0Count].fMass = v0->MassK0Short();
\r
755 tempK0[v0Count].fV0Dca = v0Dca;
\r
758 tempK0[v0Count].fDDDca = v0->DcaV0Daughters();
\r
759 tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);
\r
760 tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);
\r
761 tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);
\r
762 tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);
\r
763 tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);
\r
765 tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();
\r
766 tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();
\r
769 tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();
\r
770 tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();
\r
774 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
\r
775 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
\r
777 prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
\r
778 prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
\r
780 //if(idCount < 50){
\r
782 // idArray[idCount*2] = prongTrackPos->GetID();
\r
783 // idArray[idCount*2+1] = prongTrackNeg->GetID();
\r
792 if(k0Count<2) return; //only keep events with more than 1 good K0
\r
794 //Add Event to buffer - this is for event mixing
\r
795 fEC[zBin][centBin]->FIFOShift();
\r
796 (fEvt) = fEC[zBin][centBin]->fEvt;
\r
797 (fEvt)->fFillStatus = 1;
\r
798 int unskippedCount = 0;
\r
799 for(int i=0;i<v0Count;i++)
\r
801 if(!tempK0[i].fSkipShared) //don't include skipped v0s (from shared daughters)
\r
803 ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
\r
804 if(tempK0[i].fK0) //make sure particle is good (mass)
\r
806 (fEvt)->fK0Particle[unskippedCount] = tempK0[i]; //load good, unskipped K0s
\r
807 unskippedCount++; //count good, unskipped K0s
\r
811 (fEvt)->fNumV0s = unskippedCount;
\r
812 //Printf("Number of v0s: %d", v0Count);
\r
813 //Printf("Number of K0s: %d", k0Count);
\r
814 delete [] tempK0; tempK0 = NULL;
\r
816 ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount); // changed 3/25, used to be "k0Count"
\r
817 ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
\r
819 //Printf("Reconstruction Finished. Starting pair studies.");
\r
821 //////////////////////////////////////////////////////////////////////
\r
823 //////////////////////////////////////////////////////////////////////
\r
825 float px1, py1, pz1, px2, py2, pz2; //single kaon values
\r
826 float en1, en2; //single kaon values
\r
827 //float pt1, pt2; //single kaon values
\r
828 float pairPx, pairPy, pairPz, pairP0; //pair momentum values
\r
829 float pairPt, pairMt, pairKt; //pair momentum values
\r
830 float pairMInv, pairPDotQ;
\r
831 float qinv, q0, qx, qy, qz; //pair q values
\r
832 //float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values
\r
833 float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF
\r
835 float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame
\r
836 float betasq, gamma;
\r
837 float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
\r
838 float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
\r
841 for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
\r
843 //single particle histograms (done here to avoid "skipped" v0s
\r
844 ((TH1F*)fOutputList->FindObject("fHistDCADaughters")) ->Fill((fEvt)->fK0Particle[i].fDDDca);
\r
845 ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0")) ->Fill((fEvt)->fK0Particle[i].fDecayLength);
\r
846 ((TH1F*)fOutputList->FindObject("fHistDCAK0")) ->Fill((fEvt)->fK0Particle[i].fV0Dca);
\r
847 ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
\r
848 ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus")) ->Fill((fEvt)->fK0Particle[i].fPosDca);
\r
849 ((TH2F*)fOutputList->FindObject("fHistPtK0")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
\r
850 ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
\r
851 ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
\r
852 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fPosPhi);
\r
853 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fNegPhi);
\r
855 ((TH1F*)fOutputList->FindObject("fHistPx")) ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
\r
856 ((TH1F*)fOutputList->FindObject("fHistPy")) ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
\r
857 ((TH1F*)fOutputList->FindObject("fHistPz")) ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
\r
859 for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
\r
862 if(evnum==0) startbin=i+1;
\r
864 for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
\r
866 if(evnum==0) // Get rid of shared tracks
\r
868 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
\r
869 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
\r
870 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
\r
871 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
\r
874 px1 = (fEvt)->fK0Particle[i].fMomentum[0];
\r
875 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
\r
876 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
\r
877 //pt1 = (fEvt)->fK0Particle[i].fPt;
\r
878 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
\r
879 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
\r
880 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
\r
881 //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
\r
882 if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
\r
884 tempvar = px1; px1 = px2; px2 = tempvar;
\r
885 tempvar = py1; py1 = py2; py2 = tempvar;
\r
886 tempvar = pz1; pz1 = pz2; pz2 = tempvar;
\r
889 en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
\r
890 en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
\r
896 qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
\r
898 pairPx = px1 + px2;
\r
899 pairPy = py1 + py2;
\r
900 pairPz = pz1 + pz2;
\r
901 pairP0 = en1 + en2;
\r
902 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
\r
903 pairKt = pairPt/2.; //used for KT binning
\r
904 pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)
\r
905 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
\r
906 pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF
\r
908 //PRF (this section will probably be removed in favor of later boosting section)
\r
909 p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length
\r
910 am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)
\r
911 epm = en1+en2+am12; //"energy plus mass"
\r
912 p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP
\r
913 if(am12 == 0) continue;
\r
914 h1 = (p112/epm - en1)/am12;
\r
915 ppx = px1+pairPx*h1; //px in PRF
\r
916 ppy = py1+pairPy*h1; //py in PRF
\r
917 ppz = pz1+pairPz*h1; //pz in PRF
\r
918 ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*
\r
919 ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
\r
920 ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
\r
921 ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
\r
922 ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
\r
924 //relative momentum in out-side-long for LCMS and PRF
\r
925 if(pairMt == 0 || pairPt == 0) continue;
\r
926 qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
\r
927 qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
\r
928 //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
\r
929 qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
\r
931 //finding gamma for gamma binning/hists (likely will be removed after tests)
\r
932 p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;
\r
933 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
\r
934 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
\r
935 p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;
\r
936 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
\r
937 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
\r
938 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
\r
939 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
\r
940 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
\r
941 gamma = 1./sqrt(1-betasq);
\r
942 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
\r
943 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
\r
944 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
\r
945 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
\r
946 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
\r
947 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
\r
948 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
\r
949 //getting bin numbers and names for 3D histogram
\r
950 TString *histname3D = new TString("fHist3DOSL");
\r
952 if(pairKt < 0.6) ktBin = 0;
\r
953 else if(pairKt < 0.8) ktBin = 1;
\r
954 else if(pairKt < 1.0) ktBin = 2;
\r
956 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
\r
957 *histname3D += ktBin;
\r
959 //Spherical harmonics
\r
960 //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
\r
961 //thetaSHCos = qLong/qLength;
\r
962 //thetaSH = acos(thetaSHCos);
\r
963 //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
\r
965 //Finding average separation of daughters throughout TPC - two-track cut
\r
966 float posPositions1[9][3] = {{0}};
\r
967 float negPositions1[9][3] = {{0}};
\r
968 float posPositions2[9][3] = {{0}};
\r
969 float negPositions2[9][3] = {{0}};
\r
970 for(int iPos = 0; iPos < 9; iPos++){
\r
971 for(int jPos = 0; jPos < 3; jPos++){
\r
972 posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
\r
973 negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
\r
974 posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
\r
975 negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
\r
978 float pMean = 0.; //average separation for positive daughters
\r
979 float nMean = 0.; //average separation for negative daughters
\r
982 float pMin = 9999.; //minimum separation (updates) - pos
\r
983 float nMin = 9999.; //minimum separation (updates) - neg
\r
984 double pCount=0; //counter for number of points used - pos
\r
985 double nCount=0; //counter for number of points used - neg
\r
986 for(int ss=0;ss<9;ss++){
\r
987 if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
\r
989 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));
\r
990 pMean = pMean + pDiff;
\r
991 if(pDiff < pMin) pMin = pDiff;
\r
993 if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
\r
995 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));
\r
996 nMean = nMean + nDiff;
\r
997 if(nDiff < nMin) nMin = nDiff;
\r
1000 pMean = pMean/pCount;
\r
1001 nMean = nMean/nCount;
\r
1004 ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
\r
1005 ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
\r
1006 ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
\r
1007 ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
\r
1010 ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
\r
1011 ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
\r
1012 ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
\r
1013 ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
\r
1016 //Decay plane coincidence
\r
1017 //daughter momenta
\r
1018 float a1 = (fEvt)->fK0Particle[i].fPPos[0];
\r
1019 float b1 = (fEvt)->fK0Particle[i].fPPos[1];
\r
1020 float c1 = (fEvt)->fK0Particle[i].fPPos[2];
\r
1021 float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
\r
1022 float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
\r
1023 float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
\r
1024 float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
\r
1025 float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
\r
1026 float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
\r
1027 float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
\r
1028 float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
\r
1029 float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
\r
1033 cross1[0] = b1*f1-c1*e1;
\r
1034 cross1[1] = c1*d1-a1*f1;
\r
1035 cross1[2] = a1*e1-b1*d1;
\r
1036 cross2[0] = b2*f2-c2*e2;
\r
1037 cross2[1] = c2*d2-a2*f2;
\r
1038 cross2[2] = a2*e2-b2*d2;
\r
1039 float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
\r
1040 float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
\r
1041 float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
\r
1043 if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
\r
1044 else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
\r
1046 if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
\r
1047 //end separation studies
\r
1050 bool center1K0 = kFALSE; //accepted mass K0
\r
1051 bool center2K0 = kFALSE;
\r
1052 if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
\r
1053 if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
\r
1054 //bool CL1 = kFALSE;
\r
1055 //bool CL2 = kFALSE;
\r
1056 //bool CR1 = kFALSE;
\r
1057 //bool CR2 = kFALSE;
\r
1058 //if(center1K0 && center2K0){
\r
1059 // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
\r
1060 // else CR1 = kTRUE;
\r
1061 // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
\r
1062 // else CR2 = kTRUE;
\r
1065 //bool SideLeft1 = kFALSE;
\r
1066 //bool SideLeft2 = kFALSE;
\r
1067 //bool SideRight1 = kFALSE;
\r
1068 //bool SideRight2 = kFALSE;
\r
1069 //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
\r
1070 //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
\r
1071 //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
\r
1072 //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
\r
1075 if(evnum==0) //Same Event
\r
1077 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
\r
1078 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
\r
1079 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
\r
1080 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
\r
1081 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
\r
1082 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
\r
1084 if(center1K0 && center2K0){
\r
1086 ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1087 //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
\r
1088 ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
\r
1090 //for mass bin study
\r
1091 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1092 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1093 //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1097 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
\r
1098 histname3D->Append("Signal");
\r
1099 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
\r
1102 /*if(pairKt < 1.0){
\r
1104 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
\r
1105 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
\r
1106 else if(centBin > 9){
\r
1107 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
\r
1108 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
\r
1109 else if(pairKt < 2.0){
\r
1111 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
\r
1112 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
\r
1113 else if(centBin > 9){
\r
1114 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
\r
1116 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
\r
1120 //side-side correlations
\r
1121 //if(!splitK0sides){
\r
1122 // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1123 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1124 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
\r
1129 else //Mixed Events
\r
1131 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
\r
1132 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
\r
1133 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
\r
1134 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
\r
1136 if(center1K0 && center2K0){
\r
1138 ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1140 //for mass bin study
\r
1141 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1142 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1143 //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1147 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
\r
1148 histname3D->Replace(12,6,"Bkg");
\r
1149 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
\r
1152 /*if(pairKt < 1.0){
\r
1154 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
\r
1155 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
\r
1156 else if(centBin > 9){
\r
1157 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
\r
1158 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
\r
1159 else if(pairKt < 2.0){
\r
1161 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
\r
1162 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
\r
1163 else if(centBin > 9){
\r
1164 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
\r
1165 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
\r
1168 //side-side correlations
\r
1169 //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1170 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1171 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
\r
1179 // Post output data.
\r
1180 PostData(1, fOutputList);
\r
1183 //________________________________________________________________________
\r
1184 void AliFemtoK0Analysis::Terminate(Option_t *)
\r
1186 // Called once at the end of the query
\r
1187 cout<<"Done"<<endl;
\r
1191 //_________________________________________________________________________
\r
1192 void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
\r
1193 // Gets the global position of the track at nine different radii in the TPC
\r
1194 // track is the track you want to propagate
\r
1195 // bfield is the magnetic field of your event
\r
1196 // globalPositionsAtRadii is the array of global positions in the radii and xyz
\r
1198 // Initialize the array to something indicating there was no propagation
\r
1199 for(Int_t i=0;i<9;i++){
\r
1200 for(Int_t j=0;j<3;j++){
\r
1201 globalPositionsAtRadii[i][j]=-9999.;
\r
1205 // Make a copy of the track to not change parameters of the track
\r
1206 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
\r
1207 //printf("\nAfter CopyFromVTrack\n");
\r
1210 // The global position of the the track
\r
1211 Double_t xyz[3]={-9999.,-9999.,-9999.};
\r
1213 // Counter for which radius we want
\r
1215 // The radii at which we get the global positions
\r
1216 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
\r
1217 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
\r
1218 // The global radius we are at
\r
1219 Float_t globalRadius=0;
\r
1221 // Propagation is done in local x of the track
\r
1222 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
\r
1223 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
\r
1224 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
\r
1225 // If the track's momentum is smaller than infinite, it will develop a y-component, which
\r
1226 // adds to the global radius
\r
1228 // Stop if the propagation was not succesful. This can happen for low pt tracks
\r
1229 // that don't reach outer radii
\r
1230 if(!etp.PropagateTo(x,bfield))break;
\r
1231 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
\r
1232 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
\r
1234 // Roughly reached the radius we want
\r
1235 if(globalRadius > Rwanted[iR]){
\r
1237 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
\r
1238 while (globalRadius>Rwanted[iR]){
\r
1240 // printf("propagating to x %5.2f\n",x);
\r
1241 if(!etp.PropagateTo(x,bfield))break;
\r
1242 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
\r
1243 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
\r
1245 //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]);
\r
1246 globalPositionsAtRadii[iR][0]=xyz[0];
\r
1247 globalPositionsAtRadii[iR][1]=xyz[1];
\r
1248 globalPositionsAtRadii[iR][2]=xyz[2];
\r
1249 //subtract primary vertex, "zero" track for correct mixed-event comparison
\r
1250 globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
\r
1251 globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
\r
1252 globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
\r
1254 // Indicate we want the next radius
\r
1258 // TPC edge reached
\r
1264 bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
\r
1265 //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
\r
1266 //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
\r
1268 bool newV0Wins = kFALSE;
\r
1269 double pardiff[3] = {newPars[0]-oldPars[0],
\r
1270 newPars[1]-oldPars[1],
\r
1271 newPars[2]-oldPars[2]};
\r
1272 if(cutChoice > 0 && cutChoice < 4){
\r
1273 if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
\r
1275 else if(cutChoice == 4){
\r
1276 int newWinCount = 0;
\r
1277 for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
\r
1278 if(newWinCount > 1) newV0Wins = kTRUE;
\r
1284 bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
\r
1285 { // to flatten centrality distribution
\r
1286 bool RejectEvent = kFALSE;
\r
1287 int weightBinSign;
\r
1288 if(MagField > 0) weightBinSign = 0;
\r
1289 else weightBinSign = 1;
\r
1290 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
\r
1291 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
\r
1292 int weightBinCent = (int) CentPercent;
\r
1293 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
\r
1295 return RejectEvent;
\r