1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ////////////////////////////////////////////////////////////////////////////
18 // This class is used to perform femtoscopic analysis on K0s particles,
19 // which are reconstructed using the AliAODv0 class.
21 // authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)
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)
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)
53 // - added NPsiBins argument (1/8/14)
54 ////////////////////////////////////////////////////////////////////////////////
65 #include "TObjString.h"
76 #include "AliAnalysisTask.h"
77 #include "AliAnalysisManager.h"
79 #include "AliAODEvent.h"
80 #include "AliAODInputHandler.h"
81 #include "AliAODMCParticle.h"
83 #include "AliAODRecoDecay.h"
84 #include "AliCentrality.h"
86 #include "AliFemtoK0Analysis.h"
91 // Author: Matt Steinpreis, adapted from Dhevan Gangadharan
93 ClassImp(AliFemtoK0Analysis)
95 //________________________________________________________________________
96 AliFemtoK0Analysis::AliFemtoK0Analysis():
103 fMinDecayLength(0.0),
119 //________________________________________________________________________
120 AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool SignDep, bool FieldPositive, bool OnlineCase, bool MeritCase, bool Case3D, float MinDL, int MeritCutChoice, float MinSep, bool FlatCent, bool PsiBinning, int NPsiBins)
121 : AliAnalysisTaskSE(name),
123 fFieldPos(FieldPositive),
124 fOnlineCase(OnlineCase),
125 fMeritCase(MeritCase),
127 fMinDecayLength(MinDL),
128 fMeritCutChoice(MeritCutChoice),
131 fPsiBinning(PsiBinning),
144 fFieldPos = FieldPositive;
145 fOnlineCase = OnlineCase;
146 fMeritCase = MeritCase;
148 fMinDecayLength = MinDL;
149 fMeritCutChoice = MeritCutChoice;
151 fFlatCent = FlatCent;
152 fPsiBinning = PsiBinning;
153 fNPsiBins = NPsiBins;
155 // Define output slots here
157 DefineOutput(1, TList::Class());
160 //________________________________________________________________________
161 AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
162 : AliAnalysisTaskSE(obj.fName),
163 fSignDep(obj.fSignDep),
164 fFieldPos(obj.fFieldPos),
165 fOnlineCase(obj.fOnlineCase),
166 fMeritCase(obj.fMeritCase),
167 fCase3D(obj.fCase3D),
168 fMinDecayLength(obj.fMinDecayLength),
169 fMeritCutChoice(obj.fMeritCutChoice),
170 fMinSep(obj.fMinSep),
171 fFlatCent(obj.fFlatCent),
172 fPsiBinning(obj.fPsiBinning),
173 fNPsiBins(obj.fNPsiBins),
174 fEventCount(obj.fEventCount),
177 fRandomNumber(obj.fRandomNumber),
180 fOutputList(obj.fOutputList),
184 //________________________________________________________________________
185 AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
187 //Assignment operator
188 if (this == &obj) return *this;
190 fSignDep = obj.fSignDep;
191 fFieldPos = obj.fFieldPos;
192 fOnlineCase = obj.fOnlineCase;
193 fMeritCase = obj.fMeritCase;
194 fCase3D = obj.fCase3D;
195 fMinDecayLength= obj.fMinDecayLength;
196 fMeritCutChoice= obj.fMeritCutChoice;
197 fMinSep = obj.fMinSep;
198 fFlatCent = obj.fFlatCent;
199 fPsiBinning = obj.fPsiBinning;
200 fNPsiBins = obj.fNPsiBins;
201 fEventCount = obj.fEventCount;
204 fRandomNumber = obj.fRandomNumber;
207 fOutputList = obj.fOutputList;
208 fPidAOD = obj.fPidAOD;
212 //________________________________________________________________________
213 AliFemtoK0Analysis::~AliFemtoK0Analysis()
216 for(unsigned short i=0; i<kZVertexBins; i++)
218 for(unsigned short j=0; j<kCentBins; j++)
220 for(unsigned short k=0; k<fNPsiBins; k++)
222 fEC[i][j][k]->~AliFemtoK0EventCollection();
225 delete [] fEC[i][j]; fEC[i][j] = NULL;
227 delete[] fEC[i]; fEC[i] = NULL;
229 delete[] fEC; fEC = NULL;
231 if(fEC){ delete fEC; fEC = NULL;}
232 if(fRandomNumber){ delete fRandomNumber; fRandomNumber = NULL;}
233 if(fAOD){ delete fAOD; fAOD = NULL;}
234 if(fOutputList){ delete fOutputList; fOutputList = NULL;}
235 if(fPidAOD){ delete fPidAOD; fPidAOD = NULL;}
237 //________________________________________________________________________
238 void AliFemtoK0Analysis::MyInit()
241 // One can set global variables here
244 fEC = new AliFemtoK0EventCollection ***[kZVertexBins];
245 for(unsigned short i=0; i<kZVertexBins; i++)
247 fEC[i] = new AliFemtoK0EventCollection **[kCentBins];
249 for(unsigned short j=0; j<kCentBins; j++)
251 fEC[i][j] = new AliFemtoK0EventCollection *[fNPsiBins];
253 for(unsigned short k=0; k<fNPsiBins; k++)
255 fEC[i][j][k] = new AliFemtoK0EventCollection(kEventsToMix+1, kMultLimit);
260 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
261 fPidAOD = aodH->GetPIDResponse();
263 fRandomNumber = new TRandom3(); //for 3D, random sign switching
264 fRandomNumber->SetSeed(0);
267 //________________________________________________________________________
268 void AliFemtoK0Analysis::UserCreateOutputObjects()
273 MyInit();// Initialize my settings
275 fOutputList = new TList();
276 fOutputList->SetOwner();
278 TH1F *fHistCent = new TH1F("fHistCent","",100,0,100);
279 fOutputList->Add(fHistCent);
280 TH1F *fHistCentFlat = new TH1F("fHistCentFlat","",100,0,100);
281 fOutputList->Add(fHistCentFlat);
282 TH1F *fHistCentUsed = new TH1F("fHistCentUsed","",100,0,100);
283 fOutputList->Add(fHistCentUsed);
286 TH1F* fHistDCAPiPlus = new TH1F("fHistDCAPiPlus","",100,0,10);
287 fOutputList->Add(fHistDCAPiPlus);
288 TH1F* fHistDCAPiMinus = new TH1F("fHistDCAPiMinus","",100,0,10);
289 fOutputList->Add(fHistDCAPiMinus);
290 TH1F* fHistDCADaughters = new TH1F("fHistDCADaughters", "DCA of pions to each other", 50, 0., 0.5);
291 fOutputList->Add(fHistDCADaughters);
292 TH2F* fHistK0PiPlusPt = new TH2F("fHistK0PiPlusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
293 fOutputList->Add(fHistK0PiPlusPt);
294 TH2F* fHistK0PiMinusPt = new TH2F("fHistK0PiMinusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
295 fOutputList->Add(fHistK0PiMinusPt);
296 TH1F* fHistDaughterPhi = new TH1F("fHistDaughterPhi","",180,-PI,PI);
297 fOutputList->Add(fHistDaughterPhi);
300 TH1F* fHistMultK0 = new TH1F("fHistMultK0", "K0 multiplicity", 51, -0.5, 51-0.5);
301 fOutputList->Add(fHistMultK0);
302 TH2F* fHistPtK0 = new TH2F("fHistPtK0", "K0 pt distribution",kCentBins,.5,kCentBins+.5, 100, 0., 10.);
303 fOutputList->Add(fHistPtK0);
304 TH1F* fHistDecayLengthK0 = new TH1F("fHistDecayLengthK0", "K0 decay length", 100, 0., 100.);
305 fOutputList->Add(fHistDecayLengthK0);
306 TH1F* fHistDCAK0 = new TH1F("fHistDCAK0", "DCA of K0 to primary vertex", 40, 0., 0.4);
307 fOutputList->Add(fHistDCAK0);
308 TH2F* fHistKtK0 = new TH2F("fHistKtK0", "Kt distribution of K0 pairs", kCentBins, .5, kCentBins+.5, 300, 0., 3.);
309 fOutputList->Add(fHistKtK0);
311 TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);
312 TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);
313 TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);
314 TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);
315 TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);
316 TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);
317 TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);
318 fOutputList->Add(fHistPx);
319 fOutputList->Add(fHistPy);
320 fOutputList->Add(fHistPz);
321 fOutputList->Add(fHistPxCM);
322 fOutputList->Add(fHistPyCM);
323 fOutputList->Add(fHistPzCM);
324 fOutputList->Add(fHistKsCM);
326 TH1F* fHistPOutLCMS = new TH1F("fHistPOutLCMS","",200,0,2);
327 TH1F* fHistPSideLCMS = new TH1F("fHistPSideLCMS","",200,0,2);
328 TH1F* fHistPLongLCMS = new TH1F("fHistPLongLCMS","",200,0,2);
329 fOutputList->Add(fHistPOutLCMS);
330 fOutputList->Add(fHistPSideLCMS);
331 fOutputList->Add(fHistPLongLCMS);
333 //pair gamma (LCMS to PRF, OSL)
334 TH2F* fHistGamma = new TH2F("fHistGamma","Gamma from LCMS to PRF",500,1,5,100,0,1);
335 fOutputList->Add(fHistGamma);
337 //invariant mass distributions
338 TH3F* fHistMass = new TH3F("fHistMass","",kCentBins,.5,kCentBins+.5,50,0.,5.,400,.3,.7);
339 fOutputList->Add(fHistMass);
340 //TH3F* fHistMassPtCFK0 = new TH3F("fHistMassPtCFK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
341 //fOutputList->Add(fHistMassPtCFK0);
342 //TH3F* fHistMassPtCFBkgK0 = new TH3F("fHistMassPtCFBkgK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
343 //fOutputList->Add(fHistMassPtCFBkgK0);
344 //TH3F* fHistMassQKt = new TH3F("fHistMassQKt","",100,0,1,200,0,2,200,.4,.6);
345 //fOutputList->Add(fHistMassQKt);
346 //TH3F* fHistMassKtK0 = new TH3F("fHistMassKtK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
347 //fOutputList->Add(fHistMassKtK0);
348 //TH3F* fHistMassKtBkgK0 = new TH3F("fHistMassKtBkgK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
349 //fOutputList->Add(fHistMassKtBkgK0);
352 TH1F* fHistSepNumPos = new TH1F("fHistSepNumPos","",200,0,20);
353 fOutputList->Add(fHistSepNumPos);
354 TH1F* fHistSepDenPos = new TH1F("fHistSepDenPos","",200,0,20);
355 fOutputList->Add(fHistSepDenPos);
356 TH1F* fHistSepNumNeg = new TH1F("fHistSepNumNeg","",200,0,20);
357 fOutputList->Add(fHistSepNumNeg);
358 TH1F* fHistSepDenNeg = new TH1F("fHistSepDenNeg","",200,0,20);
359 fOutputList->Add(fHistSepDenNeg);
361 TH2F* fHistSepNumPos2 = new TH2F("fHistSepNumPos2","",100,0,20,100,0,20);
362 TH2F* fHistSepDenPos2 = new TH2F("fHistSepDenPos2","",100,0,20,100,0,20);
363 TH2F* fHistSepNumNeg2 = new TH2F("fHistSepNumNeg2","",100,0,20,100,0,20);
364 TH2F* fHistSepDenNeg2 = new TH2F("fHistSepDenNeg2","",100,0,20,100,0,20);
365 fOutputList->Add(fHistSepNumPos2);
366 fOutputList->Add(fHistSepDenPos2);
367 fOutputList->Add(fHistSepNumNeg2);
368 fOutputList->Add(fHistSepDenNeg2);
370 TH2F* fHistSepDPC = new TH2F("fHistSepDPC","",200,-1,1,50,0,10);
371 TH2F* fHistSepDPCBkg = new TH2F("fHistSepDPCBkg","",200,-1,1,50,0,10);
372 fOutputList->Add(fHistSepDPC);
373 fOutputList->Add(fHistSepDPCBkg);
375 TH1F *fHistPsi = new TH1F("fHistPsi","",90,-PI/2,PI/2);
376 fOutputList->Add(fHistPsi);
377 TH2F *fHistPhi = new TH2F("fHistPhi","",kCentBins,.5,kCentBins+.5,180,0,2*PI);
378 fOutputList->Add(fHistPhi);
379 TH2F *fHistPhiPsi = new TH2F("fHistPhiPsi","",kCentBins,.5,kCentBins+.5,180,0,2*PI);
380 fOutputList->Add(fHistPhiPsi);
383 TH3F *fHistDPhi = new TH3F("fHistDPhi","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
384 TH3F *fHistDPhiBkg = new TH3F("fHistDPhiBkg","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
385 TH3F *fHistDPhiPsi = new TH3F("fHistDPhiPsi","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
386 TH3F *fHistDPhiPsiBkg = new TH3F("fHistDPhiPsiBkg","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
387 fOutputList->Add(fHistDPhi);
388 fOutputList->Add(fHistDPhiBkg);
389 fOutputList->Add(fHistDPhiPsi);
390 fOutputList->Add(fHistDPhiPsiBkg);
392 /////////Pair Distributions///////////////////
395 TH3F* fHistQinvSignal = new TH3F("fHistQinvSignal","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
396 fOutputList->Add(fHistQinvSignal);
397 TH3F* fHistQinvBkg = new TH3F("fHistQinvBkg","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
398 fOutputList->Add(fHistQinvBkg);
401 TH3F* fHistQinvSignalEPIn = new TH3F("fHistQinvSignalEPIn","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
402 fOutputList->Add(fHistQinvSignalEPIn);
403 TH3F* fHistQinvBkgEPIn = new TH3F("fHistQinvBkgEPIn","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
404 fOutputList->Add(fHistQinvBkgEPIn);
405 TH3F* fHistQinvSignalEPOut = new TH3F("fHistQinvSignalEPOut","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
406 fOutputList->Add(fHistQinvSignalEPOut);
407 TH3F* fHistQinvBkgEPOut = new TH3F("fHistQinvBkgEPOut","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
408 fOutputList->Add(fHistQinvBkgEPOut);
410 //mass bins within peak
411 //TH3F* fHistCLCLSignal = new TH3F("fHistCLCLSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
412 //TH3F* fHistCLCLBkg = new TH3F("fHistCLCLBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
413 //TH3F* fHistCLCRSignal = new TH3F("fHistCLCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
414 //TH3F* fHistCLCRBkg = new TH3F("fHistCLCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
415 //TH3F* fHistCRCRSignal = new TH3F("fHistCRCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
416 //TH3F* fHistCRCRBkg = new TH3F("fHistCRCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
417 //fOutputList->Add(fHistCLCLSignal);
418 //fOutputList->Add(fHistCLCLBkg);
419 //fOutputList->Add(fHistCLCRSignal);
420 //fOutputList->Add(fHistCLCRBkg);
421 //fOutputList->Add(fHistCRCRSignal);
422 //fOutputList->Add(fHistCRCRBkg);
425 TH3F *fHist3DOSLSignal[10][4];
426 TH3F *fHist3DOSLBkg[10][4];
429 for(int i3D=0;i3D<10;i3D++){
430 for(int j3D=0;j3D<4;j3D++){
431 TString *histname = new TString("fHist3DOSL");
434 histname->Append("Signal");
435 fHist3DOSLSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
436 fOutputList->Add(fHist3DOSLSignal[i3D][j3D]);
437 histname->Replace(12,6,"Bkg");
438 fHist3DOSLBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
439 fOutputList->Add(fHist3DOSLBkg[i3D][j3D]);
444 //3D Spherical Harmonics
445 //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
446 //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
447 //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
448 //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
449 //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
450 //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
451 //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
452 //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
453 //fOutputList->Add(fHistSHCentLowKt);
454 //fOutputList->Add(fHistSHCentHighKt);
455 //fOutputList->Add(fHistSHSemiCentLowKt);
456 //fOutputList->Add(fHistSHSemiCentHighKt);
457 //fOutputList->Add(fHistSHCentLowKtBkg);
458 //fOutputList->Add(fHistSHCentHighKtBkg);
459 //fOutputList->Add(fHistSHSemiCentLowKtBkg);
460 //fOutputList->Add(fHistSHSemiCentHighKtBkg);
463 //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
464 //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
465 //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
466 //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
467 //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
468 //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
469 //fOutputList->Add(fHistLeftLeftSignal);
470 //fOutputList->Add(fHistLeftRightSignal);
471 //fOutputList->Add(fHistRightRightSignal);
472 //fOutputList->Add(fHistLeftLeftBkg);
473 //fOutputList->Add(fHistLeftRightBkg);
474 //fOutputList->Add(fHistRightRightBkg);
476 //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
477 //fOutputList->Add(fHistSplitK0Sides);
478 //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
479 //fOutputList->Add(fHistSplitK0Centers);
480 //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
481 //fOutputList->Add(fHistQinvSignalNoSplit);
483 PostData(1, fOutputList);
487 //________________________________________________________________________
488 void AliFemtoK0Analysis::Exec(Option_t *)
491 // Called for each event
492 //cout<<"=========== Event # "<<fEventCount+1<<" ==========="<<endl;
494 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
495 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
497 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
498 bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
499 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
501 //cout << "Failed trigger selection." << endl;
505 ///////////////////////////////////////////////////////////
507 unsigned int statusPos=0;
508 unsigned int statusNeg=0;
511 bField = fAOD->GetMagneticField();
512 if(bField == 0) return;
514 if(fFieldPos && bField < 0) return;
515 if(!fFieldPos && bField > 0) return;
519 /////////////////////////////////////////////////
522 //Centrality selection
524 AliCentrality *centrality = fAOD->GetCentrality();
525 float percent = centrality->GetCentralityPercentile("V0M");
527 //Printf("Centrality percent = %f", percent);
529 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
530 float multV0A=aodV0->GetMTotV0A();
531 float multV0C=aodV0->GetMTotV0C();
534 //Printf("No centrality info");
537 if(percent < 0.1 && (multV0A + multV0C < 19500)){
538 //Printf("No centrality info");
541 else if(percent <= 5) centBin=15;
542 else if(percent <= 10) centBin=14;
543 else if(percent <= 15) centBin=13;
544 else if(percent <= 20) centBin=12;
545 else if(percent <= 25) centBin=11;
546 else if(percent <= 30) centBin=10;
547 else if(percent <= 35) centBin=9;
548 else if(percent <= 40) centBin=8;
549 else if(percent <= 45) centBin=7;
550 else if(percent <= 50) centBin=6;
551 else if(percent <= 55) centBin=5;
552 else if(percent <= 60) centBin=4;
553 else if(percent <= 65) centBin=3;
554 else if(percent <= 70) centBin=2;
555 else if(percent <= 75) centBin=1;
556 else if(percent <= 80) centBin=0;
558 //Printf("Skipping Peripheral Event");
561 if(percent > 10 && isCentral) return;
562 ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
564 //flatten centrality dist.
567 if(RejectEventCentFlat(bField,percent)) return;
570 ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
573 AliAODVertex *primaryVertex;
574 double vertex[3]={0};
575 primaryVertex = fAOD->GetPrimaryVertex();
576 vertex[0]=primaryVertex->GetX();
577 vertex[1]=primaryVertex->GetY();
578 vertex[2]=primaryVertex->GetZ();
579 if(vertex[0]<10e-5 && vertex[1]<10e-5 && vertex[2]<10e-5) return;
580 if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
583 double zStep=2*10/double(kZVertexBins), zStart=-10.;
584 for(int i=0; i<kZVertexBins; i++)
586 if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
595 AliEventplane *eventplane = fAOD->GetEventplane();
596 if(fPsiBinning && !eventplane) return;
597 double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
598 ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
600 double psiStep = PI/double(fNPsiBins);
601 double psiStart = -0.5*PI;
602 for(int i=0; i<fNPsiBins; i++)
604 if((psiEP > psiStart+i*psiStep) && (psiEP < psiStart+(i+1)*psiStep))
610 if(!fPsiBinning) psiBin = 0;
612 ////////////////////////////////////////////////////////////////
613 //Cut Values and constants
615 //const bool kMCCase = kFALSE; //switch for MC analysis
616 const int kMaxNumK0 = 300; //maximum number of K0s, array size
617 const float kMinDCAPrimaryPion = 0.4; //minimum dca of pions to primary
618 const float kMaxDCADaughtersK0 = 0.3; //maximum dca of pions to each other - 3D
619 const float kMaxDCAK0 = 0.3; //maximum dca of K0 to primary
620 const float kMaxDLK0 = 30.0; //maximum decay length of K0
621 const float kMinDLK0 = fMinDecayLength; //minimum decay length of K0
622 const float kEtaCut = 0.8; //maximum |pseudorapidity|
623 const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle
625 const float kMinSeparation = fMinSep; //minimum daughter (pair) separation
627 const float kTOFLow = 0.8; //boundary for TOF usage
628 const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas
629 const float kMaxTPCSigmaPion = 3.0; //TPC # of sigmas
631 //const float kMassPion = .13957;
632 const float kMassK0Short = .497614; //true PDG masses
634 ////////////////////////////////////////////////////////////////
636 ////////////////////////////////////////////////////////////////
637 int v0Count = 0; //number of v0s (entries in array)
638 int k0Count = 0; //number of good K0s
640 AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
642 //for daughter sharing studies
643 //int idArray[100] = {0};
647 //TClonesArray *mcArray = 0x0;
649 //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
650 //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
652 for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
654 bool goodK0 = kFALSE;
655 bool goodPiPlus = kFALSE;
656 bool goodPiMinus = kFALSE;
659 AliAODv0* v0 = fAOD->GetV0(i);
662 if(!(v0->GetOnFlyStatus())) continue;
665 if((v0->GetOnFlyStatus())) continue; //for offline
668 //for on-the-fly ordering
669 AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
672 bool orderswitch = kFALSE;
673 if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
674 else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
676 //load daughter tracks
677 AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
678 AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
679 if(!prongTrackPos) continue;
680 if(!prongTrackNeg) continue;
683 if(v0->PtProng(pos0or1) < .15) continue;
684 if(v0->PtProng(neg0or1) < .15) continue;
685 if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
686 if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
688 //load status for PID
689 statusPos=prongTrackPos->GetStatus();
690 if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
691 prongTrackPos->SetAODEvent(fAOD);
692 statusNeg=prongTrackNeg->GetStatus();
693 if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
694 prongTrackNeg->SetAODEvent(fAOD);
697 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
698 if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
700 //Positive daughter identification TOF
702 AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
703 double Ppos = v0->PProng(pos0or1);
704 if(Ppos > kTOFLow) //PiPlus
706 //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
707 if(AliPIDResponse::kDetPidOk == statusPosTOF)
709 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
710 if(probMis < 0.01) //avoid TOF-TPC mismatch
712 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
713 else goodPiPlus = kFALSE;
717 //Negative daughter identification TOF
718 AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
719 double Pneg = v0->PProng(neg0or1);
720 if(Pneg > kTOFLow) //PiMinus
722 //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
723 if(AliPIDResponse::kDetPidOk == statusNegTOF)
725 probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
726 if(probMis < 0.01) //avoid TOF-TPC mismatch
728 if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
729 else goodPiMinus = kFALSE;
735 if(v0->Eta() > kEtaCut) continue;
736 if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;
737 if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;
738 if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;
739 if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue;
740 if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;
741 if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;
742 if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;
743 double v0Dca = v0->DcaV0ToPrimVertex();
744 if(v0Dca > kMaxDCAK0) continue;
745 if(!goodPiMinus || !goodPiPlus) continue;
747 //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
750 //bool MCgood = kFALSE;
752 //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
753 //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
754 //if(mck0dp->GetMother() >= 0){
755 //if(mck0dp->GetMother() == mck0dn->GetMother()){
756 //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
757 //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
758 //if(abs(mck0->GetPdgCode()) == 310){
766 if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
767 //else continue; //removed, Apr 18
769 //Check for shared daughters, using v0 DCA to judge
770 bool v0JudgeNew; //true if new v0 beats old
771 tempK0[v0Count].fSkipShared = kFALSE;
772 double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
774 for(int iID = 0; iID<v0Count; iID++){
775 if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old
776 if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
777 double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
778 v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
779 if(!v0JudgeNew){ //if old beats new...
780 if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing...
781 else{ //but if bad old beats new bad, or good old beats anything, skip new
782 tempK0[v0Count].fSkipShared = kTRUE; //skip new
783 break; //no need to keep checking others
786 else{ //if new beats old...
787 if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...
788 else{ //but if bad new beats bad old, or good new beats anything, skip old
789 tempK0[iID].fSkipShared = kTRUE; //skip old
790 if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
796 if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0
799 //load parameters into temporary class instance
800 if(v0Count < kMaxNumK0)
803 tempK0[v0Count].fK0 = kTRUE;
806 else tempK0[v0Count].fK0 = kFALSE;
808 //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
809 //else tempK0[v0Count].fSideLeft = kFALSE;
810 //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
811 //else tempK0[v0Count].fSideRight = kFALSE;
812 //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
814 tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();
815 tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();
816 tempK0[v0Count].fMomentum[0] = v0->Px();
817 tempK0[v0Count].fMomentum[1] = v0->Py();
818 tempK0[v0Count].fMomentum[2] = v0->Pz();
819 tempK0[v0Count].fPt = v0->Pt();
820 tempK0[v0Count].fMass = v0->MassK0Short();
821 tempK0[v0Count].fV0Dca = v0Dca;
824 tempK0[v0Count].fDDDca = v0->DcaV0Daughters();
825 tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);
826 tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);
827 tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);
828 tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);
829 tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);
831 tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();
832 tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();
835 tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();
836 tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();
840 double v0Phi = v0->Phi(); //between [0,2pi]
841 double v0PhiPsi = v0Phi-psiEP;
842 if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
843 else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
845 tempK0[v0Count].fPhi = v0Phi;
846 tempK0[v0Count].fPhiPsi = v0PhiPsi;
849 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
850 GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
852 prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
853 prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
857 // idArray[idCount*2] = prongTrackPos->GetID();
858 // idArray[idCount*2+1] = prongTrackNeg->GetID();
867 if(k0Count<2) return; //only keep events with more than 1 good K0
869 //Add Event to buffer - this is for event mixing
870 fEC[zBin][centBin][psiBin]->FIFOShift();
871 (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
872 (fEvt)->fFillStatus = 1;
873 int unskippedCount = 0;
874 for(int i=0;i<v0Count;i++)
876 if(!tempK0[i].fSkipShared) //don't include skipped v0s (from shared daughters)
878 ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
879 if(tempK0[i].fK0) //make sure particle is good (mass)
881 (fEvt)->fK0Particle[unskippedCount] = tempK0[i]; //load good, unskipped K0s
882 unskippedCount++; //count good, unskipped K0s
886 (fEvt)->fNumV0s = unskippedCount;
887 //Printf("Number of v0s: %d", v0Count);
888 //Printf("Number of K0s: %d", k0Count);
889 delete [] tempK0; tempK0 = NULL;
891 ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount); // changed 3/25, used to be "k0Count"
892 ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
894 //Printf("Reconstruction Finished. Starting pair studies.");
896 //////////////////////////////////////////////////////////////////////
898 //////////////////////////////////////////////////////////////////////
900 float px1, py1, pz1, px2, py2, pz2; //single kaon values
901 float en1, en2; //single kaon values
902 //float pt1, pt2; //single kaon values
903 float pairPx, pairPy, pairPz, pairP0; //pair momentum values
904 float pairPt, pairMt, pairKt; //pair momentum values
905 float pairMInv, pairPDotQ;
906 float qinv, q0, qx, qy, qz; //pair q values
907 //float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values
908 float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF
910 float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame
912 float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
913 float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
916 for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
918 //single particle histograms (done here to avoid "skipped" v0s
919 ((TH1F*)fOutputList->FindObject("fHistDCADaughters")) ->Fill((fEvt)->fK0Particle[i].fDDDca);
920 ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0")) ->Fill((fEvt)->fK0Particle[i].fDecayLength);
921 ((TH1F*)fOutputList->FindObject("fHistDCAK0")) ->Fill((fEvt)->fK0Particle[i].fV0Dca);
922 ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
923 ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus")) ->Fill((fEvt)->fK0Particle[i].fPosDca);
924 ((TH2F*)fOutputList->FindObject("fHistPtK0")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
925 ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
926 ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
927 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fPosPhi);
928 ((TH1F*)fOutputList->FindObject("fHistDaughterPhi")) ->Fill((fEvt)->fK0Particle[i].fNegPhi);
930 ((TH1F*)fOutputList->FindObject("fHistPx")) ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
931 ((TH1F*)fOutputList->FindObject("fHistPy")) ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
932 ((TH1F*)fOutputList->FindObject("fHistPz")) ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
934 ((TH2F*)fOutputList->FindObject("fHistPhi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhi);
935 ((TH2F*)fOutputList->FindObject("fHistPhiPsi")) ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
937 for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
940 if(evnum==0) startbin=i+1;
942 for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
944 if(evnum==0) // Get rid of shared tracks
946 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
947 if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
948 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
949 if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
952 px1 = (fEvt)->fK0Particle[i].fMomentum[0];
953 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
954 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
955 //pt1 = (fEvt)->fK0Particle[i].fPt;
956 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
957 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
958 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
959 //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
960 if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
962 tempvar = px1; px1 = px2; px2 = tempvar;
963 tempvar = py1; py1 = py2; py2 = tempvar;
964 tempvar = pz1; pz1 = pz2; pz2 = tempvar;
967 en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
968 en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
974 qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
980 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
981 pairKt = pairPt/2.; //used for KT binning
982 pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)
983 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
984 pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF
986 //PRF (this section will probably be removed in favor of later boosting section)
987 p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length
988 am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)
989 epm = en1+en2+am12; //"energy plus mass"
990 p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP
991 if(am12 == 0) continue;
992 h1 = (p112/epm - en1)/am12;
993 ppx = px1+pairPx*h1; //px in PRF
994 ppy = py1+pairPy*h1; //py in PRF
995 ppz = pz1+pairPz*h1; //pz in PRF
996 ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*
997 ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
998 ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
999 ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
1000 ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
1002 //relative momentum in out-side-long for LCMS and PRF
1003 if(pairMt == 0 || pairPt == 0) continue;
1004 qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
1005 qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
1006 //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
1007 qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
1009 //finding gamma for gamma binning/hists (likely will be removed after tests)
1010 p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;
1011 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
1012 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
1013 p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;
1014 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
1015 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
1016 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
1017 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
1018 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
1019 gamma = 1./sqrt(1-betasq);
1020 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
1021 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
1022 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
1023 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
1024 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
1025 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1026 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1027 //getting bin numbers and names for 3D histogram
1028 TString *histname3D = new TString("fHist3DOSL");
1030 if(pairKt < 0.6) ktBin = 0;
1031 else if(pairKt < 0.8) ktBin = 1;
1032 else if(pairKt < 1.0) ktBin = 2;
1034 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1035 *histname3D += ktBin;
1037 //Spherical harmonics
1038 //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1039 //thetaSHCos = qLong/qLength;
1040 //thetaSH = acos(thetaSHCos);
1041 //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1043 //Finding average separation of daughters throughout TPC - two-track cut
1044 float posPositions1[9][3] = {{0}};
1045 float negPositions1[9][3] = {{0}};
1046 float posPositions2[9][3] = {{0}};
1047 float negPositions2[9][3] = {{0}};
1048 for(int iPos = 0; iPos < 9; iPos++){
1049 for(int jPos = 0; jPos < 3; jPos++){
1050 posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1051 negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1052 posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1053 negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1056 float pMean = 0.; //average separation for positive daughters
1057 float nMean = 0.; //average separation for negative daughters
1060 float pMin = 9999.; //minimum separation (updates) - pos
1061 float nMin = 9999.; //minimum separation (updates) - neg
1062 double pCount=0; //counter for number of points used - pos
1063 double nCount=0; //counter for number of points used - neg
1064 for(int ss=0;ss<9;ss++){
1065 if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1067 pDiff = sqrt(pow(posPositions1[ss][0]-posPositions2[ss][0],2)+pow(posPositions1[ss][1]-posPositions2[ss][1],2)+pow(posPositions1[ss][2]-posPositions2[ss][2],2));
1068 pMean = pMean + pDiff;
1069 if(pDiff < pMin) pMin = pDiff;
1071 if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1073 nDiff = sqrt(pow(negPositions1[ss][0]-negPositions2[ss][0],2)+pow(negPositions1[ss][1]-negPositions2[ss][1],2)+pow(negPositions1[ss][2]-negPositions2[ss][2],2));
1074 nMean = nMean + nDiff;
1075 if(nDiff < nMin) nMin = nDiff;
1078 pMean = pMean/pCount;
1079 nMean = nMean/nCount;
1082 ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1083 ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1084 ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1085 ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1088 ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1089 ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1090 ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1091 ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1094 //Decay plane coincidence
1096 float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1097 float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1098 float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1099 float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1100 float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1101 float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1102 float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1103 float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1104 float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1105 float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1106 float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1107 float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1111 cross1[0] = b1*f1-c1*e1;
1112 cross1[1] = c1*d1-a1*f1;
1113 cross1[2] = a1*e1-b1*d1;
1114 cross2[0] = b2*f2-c2*e2;
1115 cross2[1] = c2*d2-a2*f2;
1116 cross2[2] = a2*e2-b2*d2;
1117 float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1118 float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1119 float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1121 if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1122 else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1124 if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1125 //end separation studies
1128 bool center1K0 = kFALSE; //accepted mass K0
1129 bool center2K0 = kFALSE;
1130 if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1131 if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1132 //bool CL1 = kFALSE;
1133 //bool CL2 = kFALSE;
1134 //bool CR1 = kFALSE;
1135 //bool CR2 = kFALSE;
1136 //if(center1K0 && center2K0){
1137 // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1138 // else CR1 = kTRUE;
1139 // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1140 // else CR2 = kTRUE;
1143 //bool SideLeft1 = kFALSE;
1144 //bool SideLeft2 = kFALSE;
1145 //bool SideRight1 = kFALSE;
1146 //bool SideRight2 = kFALSE;
1147 //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1148 //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1149 //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1150 //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1153 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1154 float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1155 bool inPlane1 = kFALSE;
1156 bool inPlane2 = kFALSE;
1157 if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1158 if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1159 if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1160 if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
1163 float dPhi = fabs((fEvt)->fK0Particle[i].fPhi - (fEvt+evnum)->fK0Particle[j].fPhi);
1164 if(dPhi > PI) dPhi = 2*PI-dPhi;
1165 float dPhiPsi = fabs((fEvt)->fK0Particle[i].fPhiPsi - (fEvt+evnum)->fK0Particle[j].fPhiPsi);
1166 if(dPhiPsi > PI) dPhiPsi = 2*PI-dPhiPsi;
1168 if(evnum==0) //Same Event
1170 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1171 //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1172 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1173 //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1174 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1175 //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1177 if(center1K0 && center2K0){
1179 ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1180 //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1181 ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1184 if(inPlane1 && inPlane2)
1185 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1186 else if(!inPlane1 && !inPlane2)
1187 ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
1190 ((TH3F*)fOutputList->FindObject("fHistDPhi"))->Fill(centBin+1,pairKt,dPhi);
1191 ((TH3F*)fOutputList->FindObject("fHistDPhiPsi"))->Fill(centBin+1,pairKt,dPhiPsi);
1194 //for mass bin study
1195 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1196 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1197 //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1201 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1202 histname3D->Append("Signal");
1203 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1208 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1209 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1210 else if(centBin > 9){
1211 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1212 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1213 else if(pairKt < 2.0){
1215 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1216 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1217 else if(centBin > 9){
1218 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1220 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1224 //side-side correlations
1225 //if(!splitK0sides){
1226 // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1227 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1228 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1235 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1236 //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1237 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1238 //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1240 if(center1K0 && center2K0){
1242 ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1245 if(inPlane1 && inPlane2)
1246 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1247 else if(!inPlane1 && !inPlane2)
1248 ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1250 ((TH3F*)fOutputList->FindObject("fHistDPhiBkg"))->Fill(centBin+1,pairKt,dPhi);
1251 ((TH3F*)fOutputList->FindObject("fHistDPhiPsiBkg"))->Fill(centBin+1,pairKt,dPhiPsi);
1253 //for mass bin study
1254 //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1255 //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1256 //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1260 if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1261 histname3D->Replace(12,6,"Bkg");
1262 ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1267 ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1268 ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1269 else if(centBin > 9){
1270 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1271 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1272 else if(pairKt < 2.0){
1274 ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1275 ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1276 else if(centBin > 9){
1277 ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1278 ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1281 //side-side correlations
1282 //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1283 //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1284 //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1292 // Post output data.
1293 PostData(1, fOutputList);
1296 //________________________________________________________________________
1297 void AliFemtoK0Analysis::Terminate(Option_t *)
1299 // Called once at the end of the query
1304 //_________________________________________________________________________
1305 void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1306 // Gets the global position of the track at nine different radii in the TPC
1307 // track is the track you want to propagate
1308 // bfield is the magnetic field of your event
1309 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1311 // Initialize the array to something indicating there was no propagation
1312 for(Int_t i=0;i<9;i++){
1313 for(Int_t j=0;j<3;j++){
1314 globalPositionsAtRadii[i][j]=-9999.;
1318 // Make a copy of the track to not change parameters of the track
1319 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1320 //printf("\nAfter CopyFromVTrack\n");
1323 // The global position of the the track
1324 Double_t xyz[3]={-9999.,-9999.,-9999.};
1326 // Counter for which radius we want
1328 // The radii at which we get the global positions
1329 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1330 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1331 // The global radius we are at
1332 Float_t globalRadius=0;
1334 // Propagation is done in local x of the track
1335 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1336 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1337 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1338 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1339 // adds to the global radius
1341 // Stop if the propagation was not succesful. This can happen for low pt tracks
1342 // that don't reach outer radii
1343 if(!etp.PropagateTo(x,bfield))break;
1344 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1345 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1347 // Roughly reached the radius we want
1348 if(globalRadius > Rwanted[iR]){
1350 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1351 while (globalRadius>Rwanted[iR]){
1353 // printf("propagating to x %5.2f\n",x);
1354 if(!etp.PropagateTo(x,bfield))break;
1355 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1356 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1358 //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
1359 globalPositionsAtRadii[iR][0]=xyz[0];
1360 globalPositionsAtRadii[iR][1]=xyz[1];
1361 globalPositionsAtRadii[iR][2]=xyz[2];
1362 //subtract primary vertex, "zero" track for correct mixed-event comparison
1363 globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1364 globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1365 globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1367 // Indicate we want the next radius
1377 bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1378 //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1379 //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1381 bool newV0Wins = kFALSE;
1382 double pardiff[3] = {newPars[0]-oldPars[0],
1383 newPars[1]-oldPars[1],
1384 newPars[2]-oldPars[2]};
1385 if(cutChoice > 0 && cutChoice < 4){
1386 if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1388 else if(cutChoice == 4){
1389 int newWinCount = 0;
1390 for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1391 if(newWinCount > 1) newV0Wins = kTRUE;
1397 bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1398 { // to flatten centrality distribution
1399 bool RejectEvent = kFALSE;
1401 if(MagField > 0) weightBinSign = 0;
1402 else weightBinSign = 1;
1403 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1404 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1405 int weightBinCent = (int) CentPercent;
1406 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;