]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGLF/STRANGENESS/Correlations/AliAnalysisTaskV0ChCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
From Rongrong: minor changes (binning etc)
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Correlations / AliAnalysisTaskV0ChCorrelations.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-2013, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Marek Bombara *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* The task selects candidates for K0s, Lambdas and AntiLambdas (trigger particles)
17 * and calculates correlations with charged unidentified particles (associated particles) in phi and eta.
18 * The task works with AOD (with or without MC info) events only and containes also mixing for acceptance corrections.
19 * Last update edited by Marek Bombara, August 2013, Marek.Bombara@cern.ch
20 */
21
22#include <TCanvas.h>
23#include <TList.h>
24#include <TH1.h>
25#include <TH2.h>
26#include <TH3.h>
27#include <THnSparse.h>
28#include <TObjArray.h>
29#include <TObject.h>
30
31#include "AliLog.h"
32#include "AliAnalysisManager.h"
33#include "AliAODTrack.h"
34#include "AliAODEvent.h"
35#include "AliAODv0.h"
36#include "AliAODcascade.h"
37#include "AliAODVertex.h"
38
39#include "AliMCEvent.h"
40#include "AliMCVertex.h"
41#include "AliGenEventHeader.h"
42#include "AliAODMCHeader.h"
43#include "AliAODMCParticle.h"
44#include "AliAnalyseLeadingTrackUE.h"
45
46#include "AliAODPid.h"
47#include "AliPIDResponse.h"
48#include "AliEventPoolManager.h"
49#include "AliCentrality.h"
50
51#include "AliAnalysisTaskV0ChCorrelations.h"
52#include "AliPhysicsSelectionTask.h"
53
54#include <AliMultiInputEventHandler.h>
55#include <AliMixInputEventHandler.h>
56
57ClassImp(AliAnalysisTaskV0ChCorrelations)
58ClassImp(AliV0ChBasicParticle)
59
60//________________________________________________________________________
61AliAnalysisTaskV0ChCorrelations::AliAnalysisTaskV0ChCorrelations(const char *name) // All data members should be initialised here
62 : AliAnalysisTaskSE(name),
63 fAnalysisMC(kFALSE),
64 fFillMixed(kTRUE),
65 fMixingTracks(50000),
66 fPoolMgr(0x0),
67
68 fDcaDToPV(0.5),
69 fDcaV0D(0.1),
70 fWithChCh(kFALSE),
71 fOStatus(1),
72
73 fOutput(0),
74 fPIDResponse(0),
75 fHistCentVtx(0),
76 fHistMultiMain(0),
77 fHistMassK0(0),
78 fHistMassLambda(0),
79 fHistMassAntiLambda(0),
80
81 fHistdPhidEtaSib(0),
82 fHistdPhidEtaMix(0),
83 fHistNTrigSib(0),
84
85 fHistMCPtCentTrig(0),
86 fHistRCPtCentTrig(0),
87 fHistMCPtCentAs(0),
88 fHistRCPtCentAs(0),
89 fHistRCPtCentAll(0),
90
91 fHistTemp(0),
92 fHistTemp2(0)// The last in the above list should not have a comma after it
93{
94 // Constructor
95 // Define input and output slots here (never in the dummy constructor)
96 // Input slot #0 works with a TChain - it is connected to the default input container
97 // Output slot #1 writes into a TH1 container
98 DefineOutput(1, TList::Class()); // for output list
99}
100
101//________________________________________________________________________
102AliAnalysisTaskV0ChCorrelations::~AliAnalysisTaskV0ChCorrelations()
103{
104 // Destructor. Clean-up the output list, but not the histograms that are put inside
105 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
106 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
107 delete fOutput;
108 }
109}
110
111//________________________________________________________________________
112void AliAnalysisTaskV0ChCorrelations::UserCreateOutputObjects()
113{
114 // Create histograms
115
116 fOutput = new TList();
117 fOutput->SetOwner(); // IMPORTANT!
118 // defining bins for centrality
119 Int_t nCentralityBins = 9;
120 Double_t centBins[] = {0.,10.,20.,30.,40.,50.,60.,70.,80.,90.};
121 const Double_t* centralityBins = centBins;
122 // defining bins for Z vertex
123 Int_t nZvtxBins = 7;
124 //Double_t vertexBins[] = {-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.};
125 Double_t vertexBins[] = {-7.,-5.,-3.,-1.,1.,3.,5.,7.};
126 const Double_t* zvtxBins = vertexBins;
127 // pt bins of associated particles for the analysis
128 Int_t nPtBins = 7;
129 const Double_t PtBins[8] = {2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0};
130 //Int_t nPtBins = 1;
131 //const Double_t PtBins[2] = {3.0,15.0};
132 // pt bins of trigger particles for the analysis
133 Int_t nPtBinsV0 = 11;
134 //const Double_t PtBinsV0[2] = {6.0,15.0};
135 const Double_t PtBinsV0[12] = {4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0};
136 // V0 candidate: 1 - K0sig, 2 - Lamsig, 3 - Alamsig, 4 - K0bg, 5 - Lambg, 6 - Alambg
137 Int_t nTrigC = 7;
138 const Double_t TrigC[8] = {0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5};
139
140 // Create general histograms
141 fHistCentVtx = new TH2F("fHistCentVtx", "Centrality vs. Z vertex", nCentralityBins, centralityBins, nZvtxBins, zvtxBins);
142 fHistMultiMain = new TH1F("fHistMultiMain", "Multiplicity of main events", 2000, 0, 2000);
143
144 const Int_t mrBins[3] = {nPtBinsV0, nCentralityBins, nTrigC};
145 const Double_t mrMin[3] = {PtBinsV0[0], centralityBins[0], TrigC[0]};
146 const Double_t mrMax[3] = {PtBinsV0[11], centralityBins[9], TrigC[6]};
147
148 // Create histograms for reconstruction track and V0 efficiency
149 fHistMCPtCentTrig = new THnSparseF("fHistMCPtCentTrig", "MC Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);
150 fHistRCPtCentTrig = new THnSparseF("fHistRCPtCentTrig", "Rec Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);
151
152 // pt bins of associated particles for the efficiency
153 Int_t nPtBinsAs = 13;
154 const Double_t PtBinsAs[14] = {2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0};
155
156 fHistMCPtCentAs = new TH2D("fHistMCPtCentAs", "MC Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
157 fHistRCPtCentAs = new TH2D("fHistRCPtCentAs", "Rec Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
158 fHistRCPtCentAll = new TH2D("fHistRCPtCentAll", "Rec Pt vs. Cent. All Prim+Sec", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
159
160 // defining bins for mass distributions
161 Int_t nBins = 200;
162 Double_t mk0Min = 0.40;
163 Double_t mk0Max = 0.60;
164 Double_t mlaMin = 1.07;
165 Double_t mlaMax = 1.17;
166 Double_t malMin = 1.07;
167 Double_t malMax = 1.17;
168
169 const Int_t spBins[3] = {nBins, nPtBinsV0, nCentralityBins};
170 const Double_t spMinK0[3] = {mk0Min, PtBinsV0[0], centralityBins[0]};
171 const Double_t spMaxK0[3] = {mk0Max, PtBinsV0[11], centralityBins[9]};
172 const Double_t spMinLa[3] = {mlaMin, PtBinsV0[0], centralityBins[0]};
173 const Double_t spMaxLa[3] = {mlaMax, PtBinsV0[11], centralityBins[9]};
174 const Double_t spMinAl[3] = {malMin, PtBinsV0[0], centralityBins[0]};
175 const Double_t spMaxAl[3] = {malMax, PtBinsV0[11], centralityBins[9]};
176 // Create mass histograms
177 fHistMassK0 = new THnSparseF("fHistMassK0","V0 mass for K0 hypothesis", 3, spBins, spMinK0, spMaxK0);
178 fHistMassLambda = new THnSparseF("fHistMassLambda","V0 mass for Lambda hypothesis", 3, spBins, spMinLa, spMaxLa);
179 fHistMassAntiLambda = new THnSparseF("fHistMassAntiLambda","V0 mass for AntiLambda hypothesis", 3, spBins, spMinAl, spMaxAl);
180
181 // defining bins for dPhi distributions
182 const Int_t nbPhiBins = 72;
183 //const Double_t kPi = TMath::Pi();
184 //Double_t PhiMin = -kPi/2.;
185 //Double_t PhiMax = -kPi/2. + 2*kPi;
186 Double_t PhiMin = -1.57;
187 Double_t PhiMax = 4.71;
188 Double_t PhiBins[nbPhiBins+1] = {0.};
189 PhiBins[0] = PhiMin;
190 for (Int_t i=0; i<nbPhiBins; i++) { PhiBins[i+1] = PhiBins[i] + (PhiMax - PhiMin)/nbPhiBins; }
191
192 // defining bins for dEta distributions
193 const Int_t nbEtaBins = 40;
194 Double_t EtaMin = -2.0;
195 Double_t EtaMax = 2.0;
196 Double_t EtaBins[nbEtaBins+1] = {0.};
197 EtaBins[0] = EtaMin;
198 for (Int_t i=0; i<nbEtaBins; i++) { EtaBins[i+1] = EtaBins[i] + (EtaMax - EtaMin)/nbEtaBins; }
199
200 const Int_t corBins[7] = {nbPhiBins, nbEtaBins, nPtBinsV0, nPtBins, nCentralityBins, nZvtxBins, nTrigC};
201 const Double_t corMin[7] = {PhiBins[0], EtaBins[0], PtBinsV0[0], PtBins[0], centralityBins[0], zvtxBins[0], TrigC[0]};
202 const Double_t corMax[7] = {PhiBins[72], EtaBins[40], PtBinsV0[11], PtBins[7], centralityBins[9], zvtxBins[7], TrigC[7]};
203 // Create correlation histograms
204 fHistdPhidEtaSib = new THnSparseF("fHistdPhidEtaSib","dPhi vs. dEta siblings", 7, corBins, corMin, corMax);
205 fHistdPhidEtaMix = new THnSparseF("fHistdPhidEtaMix","dPhi vs. dEta mixed", 7, corBins, corMin, corMax);
206
207 // Create histograms for counting the numbers of trigger particles
208 const Int_t corNTrigBins[5] = {nPtBinsV0, nCentralityBins, nZvtxBins, nTrigC};
209 const Double_t corNTrigMin[5] = {PtBinsV0[0], centBins[0], vertexBins[0], TrigC[0]};
210 const Double_t corNTrigMax[5] = {PtBinsV0[11], centBins[9], vertexBins[7], TrigC[7]};
211 fHistNTrigSib = new THnSparseF("fHistNTrigSib","Number trigger sib", 4, corNTrigBins, corNTrigMin, corNTrigMax);
212
213 // Histograms for debugging
214 fHistTemp = new TH1D("fHistTemp", "Temporary", 100, -10., 10.);
215 fHistTemp2 = new TH1D("fHistTemp2", "Temporary2", 100, -10., 10.);
216
217 fOutput->Add(fHistCentVtx);
218
219 fOutput->Add(fHistMultiMain);
220 fOutput->Add(fHistMassK0);
221 fOutput->Add(fHistMassLambda);
222 fOutput->Add(fHistMassAntiLambda);
223
224 fOutput->Add(fHistdPhidEtaSib);
225 fOutput->Add(fHistdPhidEtaMix);
226 fOutput->Add(fHistNTrigSib);
227
228 fOutput->Add(fHistMCPtCentTrig);
229 fOutput->Add(fHistRCPtCentTrig);
230 fOutput->Add(fHistMCPtCentAs);
231 fOutput->Add(fHistRCPtCentAs);
232 fOutput->Add(fHistRCPtCentAll);
233
234 fOutput->Add(fHistTemp);
235 fOutput->Add(fHistTemp2);
236
237 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
238
239 // Settings for event mixing -------------------------------------
240 Int_t trackDepth = fMixingTracks;
241 Int_t poolSize = 200; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
242
243 fPoolMgr = new AliEventPoolManager(poolSize, trackDepth, nCentralityBins, centBins, nZvtxBins, vertexBins);
244 //----------------------------------------------
245}
246
247//________________________________________________________________________
248void AliAnalysisTaskV0ChCorrelations::Terminate(Option_t *)
249{
250 // Draw result to screen, or perform fitting, normalizations
251 // Called once at the end of the query
252
253 fOutput = dynamic_cast<TList*>(GetOutputData(1));
254 if (!fOutput) { AliError("Could not retrieve TList fOutput"); return; }
255
256 // NEW HISTO should be retrieved from the TList container in the above way,
257 // so it is available to draw on a canvas such as below
258}
259
260//_________________________________________________________________________
261void AliAnalysisTaskV0ChCorrelations::UserExec(Option_t *)
262{
263 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
264 AliInputEventHandler *inEvMain = (AliInputEventHandler*)mgr->GetInputEventHandler();
265
266 // physics selection
267 UInt_t maskIsSelected = inEvMain->IsEventSelected();
268 // 2010 data trigger selection
269 //Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
270 // 2011 data trigger selection
271 Bool_t isSelected = ((maskIsSelected & AliVEvent::kMB) || (maskIsSelected & AliVEvent::kCentral) || (maskIsSelected & AliVEvent::kSemiCentral));
272 if (!isSelected) return;
273 // calculating the event types
274 if (maskIsSelected & AliVEvent::kMB) fHistTemp->Fill(2);
275 if (maskIsSelected & AliVEvent::kCentral) fHistTemp->Fill(4);
276 if (maskIsSelected & AliVEvent::kSemiCentral) fHistTemp->Fill(6);
277
278 AliAODEvent* aod = (AliAODEvent*)inEvMain->GetEvent();
279 fPIDResponse = inEvMain->GetPIDResponse();
280
281 // pt intervals for trigger particles
282 const Double_t kPi = TMath::Pi();
283 Double_t PtTrigMin = 4.0;
284 Double_t PtTrigMax = 15.0;
285 // pt interval for associated particles
286 Double_t PtAssocMin = 2.0;
287 fHistMultiMain->Fill(aod->GetNumberOfTracks());
288
289 // Vertex cut
290 Double_t cutPrimVertex = 7.0;
291 AliAODVertex *myPrimVertex = aod->GetPrimaryVertex();
292 if (!myPrimVertex) return;
293 if ( ( TMath::Abs(myPrimVertex->GetZ()) ) >= cutPrimVertex) return ;
294 Double_t lPVx = myPrimVertex->GetX();
295 Double_t lPVy = myPrimVertex->GetY();
296 Double_t lPVz = myPrimVertex->GetZ();
297
298 if (TMath::Abs(lPVx)<10e-5 && TMath::Abs(lPVy)<10e-5 && TMath::Abs(lPVz)<10e-5) return;
299
300 // Centrality definition
301 Double_t lCent = 0.0;
302 AliCentrality *centralityObj = 0;
303 centralityObj = aod->GetHeader()->GetCentralityP();
304 lCent = centralityObj->GetCentralityPercentile("V0M");
305 if ((lCent < 0.)||(lCent > 90.)) return;
306 fHistCentVtx->Fill(lCent,lPVz);
307
308//=========== MC loop ===============================
309 if (fAnalysisMC)
310 {
311 AliAODMCHeader *aodMCheader = (AliAODMCHeader*)aod->FindListObject(AliAODMCHeader::StdBranchName());
312 Float_t vzMC = aodMCheader->GetVtxZ();
313 if (TMath::Abs(vzMC) >= 7.) return;
314 //retrieve MC particles from event
315 TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
316 if(!mcArray){
317 Printf("No MC particle branch found");
318 return;
319 }
320
321 Int_t nMCAllTracks = mcArray->GetEntriesFast();
322 // new tracks array - without injected signal
323 TObjArray * mcTracks = new TObjArray;
324 //selectedMCTracks->SetOwner(kTRUE);
325
326 for (Int_t i = 0; i < nMCAllTracks; i++)
327 {
328 AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcArray->At(i);
329 if (!mcTrack) {
330 Error("ReadEventAODMC", "Could not receive particle %d", i);
331 continue;
332 }
333 mcTracks->Add(mcTrack);
334 }
335
336 // get labels for injected particles -------
337 TObject* mc = mcArray;
338 //TObjArray * mcTracks = mcArray;
339 Int_t skipParticlesAbove = 0;
340 //AliGenEventHeader* eventHeader = 0;
341 //Int_t headers = 0;
342 //headers = aodMCheader->GetNCocktailHeaders();
343 //eventHeader = aodMCheader->GetCocktailHeader(0);
344 AliGenEventHeader* eventHeader = aodMCheader->GetCocktailHeader(0);
345 Int_t headers = aodMCheader->GetNCocktailHeaders();
346 if (!eventHeader)
347 {
348 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
349 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
350 AliError("First event header not found. Skipping this event.");
351 return;
352 }
353 skipParticlesAbove = eventHeader->NProduced();
354 //cout << "ski label = " << skipParticlesAbove << endl;
355 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s.Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
356 //cout << "before MC compressing = " << mcTracks->GetEntriesFast() << endl;
357 RemovingInjectedSignal(mcTracks,mc,skipParticlesAbove);
358 //cout << "after MC compressing = " << mcTracks->GetEntriesFast() << endl;
359 // -------------
360
361 Int_t nMCTracks = mcTracks->GetEntriesFast();
362 //cout << "number of MC tracks = " << nMCTracks << endl;
363 for (Int_t iMC = 0; iMC<nMCTracks; iMC++)
364 {
365 AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcTracks->At(iMC);
366 if (!mcTrack) {
367 Error("ReadEventAODMC", "Could not receive particle %d", iMC);
368 continue;
369 }
370 // track part
371 Double_t mcTrackEta = mcTrack->Eta();
372 Double_t mcTrackPt = mcTrack->Pt();
373 Bool_t TrIsPrim = mcTrack->IsPhysicalPrimary();
374 Bool_t TrEtaMax = TMath::Abs(mcTrackEta)<0.9;
375 Bool_t TrPtMin = mcTrackPt>PtAssocMin;
376 Bool_t TrCharge = (mcTrack->Charge())!=0;
377 if (TrIsPrim && TrEtaMax && TrPtMin && TrCharge) fHistMCPtCentAs->Fill(mcTrackPt,lCent);
378 // V0 part
379 Int_t mcMotherPdg = 0;
380 Int_t mcPartPdg = mcTrack->GetPdgCode();
381
382 // Keep only K0s, Lambda and AntiLambda
383 if ((mcPartPdg != 310) && (mcPartPdg != 3122) && (mcPartPdg != (-3122))) continue;
384
385 //cout << " mc pdg is " << mcPartPdg << endl;
386
387 Bool_t IsK0 = mcPartPdg==310;
388 Bool_t IsLambda = mcPartPdg==3122;
389 Bool_t IsAntiLambda = mcPartPdg==-3122;
390 Bool_t IsSigma = kFALSE;
391 Int_t mcMotherLabel = mcTrack->GetMother();
392 AliAODMCParticle *mcMother = (AliAODMCParticle*)mcArray->At(mcMotherLabel);
393 if (mcMotherLabel < 0) {mcMotherPdg = 0;} else {mcMotherPdg = mcMother->GetPdgCode();}
394 //if ((mcMotherLabel >= 0) && mcMother)
395 //if ((mcMotherLabel >= 0) && mcMother)
396 //{
397 // Bool_t IsFromCascade = ((mcMotherPdg==3312)||(mcMotherPdg==-3312)||(mcMotherPdg==3334)||(mcMotherPdg==-3334));
398 Bool_t IsFromSigma = ((mcMotherPdg==3212)||(mcMotherPdg==-3212));
399 IsFromSigma = IsFromSigma || ((mcMotherPdg==3224)||(mcMotherPdg==-3224));
400 IsFromSigma = IsFromSigma || ((mcMotherPdg==3214)||(mcMotherPdg==-3214));
401 IsFromSigma = IsFromSigma || ((mcMotherPdg==3114)||(mcMotherPdg==-3114));
402 if ((IsFromSigma) && (mcMother->IsPhysicalPrimary()) && (IsLambda || IsAntiLambda)) IsSigma = kTRUE;
403 Double_t mcRapidity = mcTrack->Y();
404 Bool_t V0RapMax = TMath::Abs(mcRapidity)<0.75;
405 Bool_t PtInterval = ((mcTrackPt>PtTrigMin)&&(mcTrackPt<PtTrigMax));
406 //Bool_t PtInterval = kTRUE;
407 IsK0 = IsK0 && (mcTrack->IsPhysicalPrimary());
408 IsLambda = IsLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
409 IsAntiLambda = IsAntiLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
410 Double_t mcK0[3] = {mcTrackPt, lCent, 1};
411 Double_t mcLa[3] = {mcTrackPt, lCent, 2};
412 Double_t mcAl[3] = {mcTrackPt, lCent, 3};
413 if (IsK0 && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcK0);
414 if (IsLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcLa);
415 if (IsAntiLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcAl);
416 //}
417 }
418 // ------- access the real data -----------
419 Int_t nTracks = aod->GetNumberOfTracks();
420 // new tracks array - without injected signal
421 TObjArray * selectedMCTracks = new TObjArray;
422 //selectedMCTracks->SetOwner(kTRUE);
423
424 for (Int_t i = 0; i < nTracks; i++)
425 {
426 AliAODTrack* tr = aod->GetTrack(i);
427 selectedMCTracks->Add(tr);
428 }
429 // cout << "before compressing = " << selectedMCTracks->GetEntriesFast() << endl;
430 RemovingInjectedSignal(selectedMCTracks,mc,skipParticlesAbove);
431 //AliAnalyseLeadingTrackUE::RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
432 //RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
433 // cout << "after compressing = " << selectedMCTracks->GetEntriesFast() << endl;
434 //-----------------------
435 //cout << "before getting label -4" << endl;
436 Int_t nRecTracks = selectedMCTracks->GetEntriesFast();
437 //cout << "before getting label -3" << endl;
438 for (Int_t i = 0; i < nRecTracks; i++)
439 {
440 //AliAODTrack* tras = aod->GetTrack(i);
441 AliAODTrack* tras = (AliAODTrack*)selectedMCTracks->At(i);
442 //cout << "before getting label -2" << endl;
443 //cout << " and i = " << i << endl;
444 //cout << "pt = " << tras->Pt() << " and i = " << i << endl;
445 if ((tras->Pt())<PtAssocMin) continue;
446 //cout << "before getting label -1" << endl;
447 if (!(IsMyGoodPrimaryTrack(tras))) continue;
448 //cout << "before getting label 0" << endl;
449 Int_t AssocLabel = tras->GetLabel();
450 if (AssocLabel<=0) continue;
451 //cout << "before getting label 1" << endl;
452 Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
453 //Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
454 //if(!(static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary())) continue;
455 //cout << "before getting label 2" << endl;
456 Double_t mcPt = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Pt();
457 //cout << "before getting label 3" << endl;
458 //Double_t mcEta = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Eta();
459 //if (mcPt<PtAssocMin) continue;
460 //if (TMath::Abs(mcEta)>0.8) continue;
461 //Double_t trPt = tras->Pt();
462 fHistRCPtCentAll->Fill(mcPt,lCent);
463 if (isPhyPrim) fHistRCPtCentAs->Fill(mcPt,lCent);
464 }
465 // ------- end of real data access, for V0s see the main V0 loop --------
466 }
467//============= End of MC loop ======================
468
469 // Track selection loop
470 //--------------------------------
471 Int_t nTracks = aod->GetNumberOfTracks();
472 // new tracks array
473 TObjArray * selectedTracks = new TObjArray;
474 selectedTracks->SetOwner(kTRUE);
475
476 Bool_t ChChWith = GetWithChCh();
477 TObjArray * selectedChargedTriggers = new TObjArray;
478 selectedChargedTriggers->SetOwner(kTRUE);
479 for (Int_t i = 0; i < nTracks; i++)
480 {
481 AliAODTrack* tr = aod->GetTrack(i);
482 if ((tr->Pt())<PtAssocMin) continue;
483 if (!(IsMyGoodPrimaryTrack(tr))) continue;
484 selectedTracks->Add(tr);
485 // saving the Charged trigger particles
486 if ((tr->Pt()>4.)&&(tr->Pt()<15.))
487 {
488 selectedChargedTriggers->Add(new AliV0ChBasicParticle(tr->Eta(), tr->Phi(), tr->Pt(), 7));
489 }
490 }
491
492 //---------------------------------
493
494 // V0 loop for reconstructed event
495 TObjArray * selectedV0s = new TObjArray;
496 selectedV0s->SetOwner(kTRUE);
497
498 Int_t nV0s = aod->GetNumberOfV0s();
499
500 for (Int_t i = 0; i < nV0s; i++)
501 { // start of V0 slection loop
502 AliAODv0* aodV0 = dynamic_cast<AliAODv0 *>(aod->GetV0(i));
503 if (!aodV0) {
504 AliError(Form("ERROR: Could not retrieve aodv0 %d", i));
505 continue;
506 }
507 if (aodV0->GetOnFlyStatus()) fHistTemp->Fill(6.);
508 if (!aodV0->GetOnFlyStatus()) fHistTemp->Fill(8.);
509 //cout << "pt of v0: " << aodV0->Pt() << endl;
510 if (((aodV0->Pt())<PtTrigMin)||((aodV0->Pt())>PtTrigMax)) continue;
511 // get daughters
512 const AliAODTrack *myTrackPos;
513 const AliAODTrack *myTrackNeg;
514 AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
515 AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));
516
517 if (!myTrackPosTest || !myTrackNegTest) {
518 Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
519 continue;
520 }
521
522 if( myTrackPosTest->Charge() ==1){
523 myTrackPos = myTrackPosTest;
524 myTrackNeg = myTrackNegTest;
525 }
526
527 if( myTrackPosTest->Charge() ==-1){
528 myTrackPos = myTrackNegTest;
529 myTrackNeg = myTrackPosTest;
530 }
531
532 // effective mass calculations for each hypothesis
533 Double_t lInvMassK0 = aodV0->MassK0Short();
534 Double_t lInvMassAntiLambda = aodV0->MassAntiLambda();
535 Double_t lInvMassLambda = aodV0->MassLambda();
536
537 // calculations for c*tau cut--------------------------------------
538 // Double_t lDVx = aodV0->GetSecVtxX();
539 // Double_t lDVy = aodV0->GetSecVtxY();
540 // Double_t lDVz = aodV0->GetSecVtxZ();
541 // const Double_t kLambdaMass = 1.115683;
542 // const Double_t kK0Mass = 0.497648;
543 // Double_t cutcTauLam = 3*7.89;
544 // Double_t cutcTauK0 = 3*2.68;
545 // Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lDVx - lPVx,2) + TMath::Power(lDVy- lPVy,2) + TMath::Power(lDVz - lPVz,2 ));
546 // Double_t lPV0 = TMath::Sqrt((aodV0->Pt())*(aodV0->Pt())+(aodV0->Pz())*(aodV0->Pz()));
547 // Double_t lcTauLam = (lV0DecayLength*kLambdaMass)/lPV0;
548 // Double_t lcTauK0 = (lV0DecayLength*kK0Mass)/lPV0;
549 // sc - standard cuts
550 //Bool_t cutK0sc = (lcTauK0<cutcTauK0);
551 //Bool_t cutLambdasc = (lcTauLam<cutcTauLam);
552 //Bool_t cutAntiLambdasc = (lcTauLam<cutcTauLam);
553 Bool_t cutK0sc = kTRUE;
554 Bool_t cutLambdasc = kTRUE;
555 Bool_t cutAntiLambdasc = kTRUE;
556
557 //------------------------------------------------
558
559 // preparations for daughter's PID cut------------------------------
560 Float_t nSigmaPosPion = 0.;
561 Float_t nSigmaNegPion = 0.;
562 Float_t nSigmaPosProton = 0.;
563 Float_t nSigmaNegProton = 0.;
564
565 const AliAODPid *pPid = myTrackPos->GetDetPid();
566 const AliAODPid *nPid = myTrackNeg->GetDetPid();
567
568 if (pPid)
569 {
570 Double_t pdMom = pPid->GetTPCmomentum();
571 if (pdMom<1.)
572 {
573 nSigmaPosPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
574 nSigmaPosProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
575 }
576 }
577
578 if (nPid)
579 {
580 Double_t ndMom = nPid->GetTPCmomentum();
581 if (ndMom<1.)
582 {
583 nSigmaNegPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
584 nSigmaNegProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
585 }
586 }
587 Bool_t bpPion = TMath::Abs(nSigmaPosPion) <= 3.;
588 Bool_t bpProton = TMath::Abs(nSigmaPosProton) <= 3.;
589 Bool_t bnPion = TMath::Abs(nSigmaNegPion) <= 3.;
590 Bool_t bnProton = TMath::Abs(nSigmaNegProton) <= 3.;
591
592 Bool_t cutK0Pid = (bpPion && bnPion);
593 Bool_t cutLambdaPid = (bpProton && bnPion);
594 Bool_t cutAntiLambdaPid = (bpPion && bnProton);
595 //--------------------------------------------------
596 // mass cuts
597 Bool_t cutMassLambda = ((lInvMassLambda>1.10) && (lInvMassLambda<1.13));
598 Bool_t cutMassAntiLambda = ((lInvMassAntiLambda>1.10) && (lInvMassAntiLambda<1.13));
599 Bool_t cutMassK0 = (lInvMassK0>0.47) && (lInvMassK0<0.53);
600
601 cutK0sc = cutK0Pid && (!cutMassLambda) && (!cutMassAntiLambda);
602 cutLambdasc = cutLambdaPid && (!cutMassK0);
603 cutAntiLambdasc = cutAntiLambdaPid && (!cutMassK0);
604 // special cut related to AP diagram for K0s
605 Bool_t k0APcut = (aodV0->PtArmV0()>(TMath::Abs(0.2*aodV0->AlphaV0())));
606 cutK0sc = cutK0sc && k0APcut;
607 // fill the mass histograms
608
609 Int_t oStatus = GetOStatus();
610
611 if (!IsMyGoodV0(aod,aodV0,myTrackPos,myTrackNeg,oStatus)) continue;
612 Double_t spK0[3] = {lInvMassK0, aodV0->Pt(), lCent};
613 Double_t spLa[3] = {lInvMassLambda, aodV0->Pt(), lCent};
614 Double_t spAl[3] = {lInvMassAntiLambda, aodV0->Pt(), lCent};
615 if (cutK0sc) fHistMassK0->Fill(spK0);
616 if (cutLambdasc) fHistMassLambda->Fill(spLa);
617 if (cutAntiLambdasc) fHistMassAntiLambda->Fill(spAl);
618 // select final V0s for correlation, selected background to study its contribution to correlation function
619 // the values for signal might change in the future ...
620 Bool_t K0Signal = (lInvMassK0>0.48)&&(lInvMassK0<0.52);
621 Bool_t K0Bckg = ((lInvMassK0>0.40)&&(lInvMassK0<0.44)) || ((lInvMassK0>0.56)&&(lInvMassK0<0.60));
622
623 Bool_t LamSignal = (lInvMassLambda>1.108)&&(lInvMassLambda<1.125);
624 Bool_t LamBckg = ((lInvMassLambda>1.090)&&(lInvMassLambda<1.100)) || ((lInvMassLambda>1.135)&&(lInvMassLambda<1.145));
625
626 Bool_t ALamSignal = (lInvMassAntiLambda>1.108)&&(lInvMassAntiLambda<1.125);
627 Bool_t ALamBckg = ((lInvMassAntiLambda>1.090)&&(lInvMassAntiLambda<1.100)) || ((lInvMassAntiLambda>1.135)&&(lInvMassAntiLambda<1.145));
628
629 // Fill selected V0 particle array
630 if ((cutK0sc)&&(K0Signal))
631 {
632 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 1.));
633 Double_t nTrigK0Sig[4] = {aodV0->Pt(), lCent, lPVz, 1.};
634 fHistNTrigSib->Fill(nTrigK0Sig);
635 }
636 if ((cutK0sc)&&(K0Bckg))
637 {
638 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 4.));
639 Double_t nTrigK0Bkg[4] = {aodV0->Pt(), lCent, lPVz, 4.};
640 fHistNTrigSib->Fill(nTrigK0Bkg);
641 }
642 if ((cutLambdasc)&&(LamSignal))
643 {
644 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 2.));
645 Double_t nTrigLaSig[4] = {aodV0->Pt(), lCent, lPVz, 2.};
646 fHistNTrigSib->Fill(nTrigLaSig);
647 }
648 if ((cutLambdasc)&&(LamBckg))
649 {
650 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 5.));
651 Double_t nTrigLaBkg[4] = {aodV0->Pt(), lCent, lPVz, 5.};
652 fHistNTrigSib->Fill(nTrigLaBkg);
653 }
654 if ((cutAntiLambdasc)&&(ALamSignal))
655 {
656 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 3.));
657 Double_t nTrigAlSig[4] = {aodV0->Pt(), lCent, lPVz, 3.};
658 fHistNTrigSib->Fill(nTrigAlSig);
659 }
660 if ((cutAntiLambdasc)&&(ALamBckg))
661 {
662 selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 6.));
663 Double_t nTrigAlBkg[4] = {aodV0->Pt(), lCent, lPVz, 6.};
664 fHistNTrigSib->Fill(nTrigAlBkg);
665 }
666
667
668 //===== MC part for V0 reconstruction efficiency ==============
669 if (fAnalysisMC)
670 {
671 TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
672 if(!mcArray){
673 Printf("No MC particle branch found");
674 return;
675 }
676
677 Int_t MotherOfMotherPdg = 0;
678
679 Int_t myTrackPosLabel = TMath::Abs(myTrackPos->GetLabel());
680 Int_t myTrackNegLabel = TMath::Abs(myTrackNeg->GetLabel());
681 AliAODMCParticle *mcPosTrack = (AliAODMCParticle*)mcArray->At(myTrackPosLabel);
682 if (!mcPosTrack) continue;
683 AliAODMCParticle *mcNegTrack = (AliAODMCParticle*)mcArray->At(myTrackNegLabel);
684 if (!mcNegTrack) continue;
685
686 Int_t PosTrackPdg = mcPosTrack->GetPdgCode();
687 Int_t NegTrackPdg = mcNegTrack->GetPdgCode();
688 //if (!mcPosTrack->IsPrimary()) continue;
689 //if (!mcNegTrack->IsPrimary()) continue;
690
691 Int_t myTrackPosMotherLabel = mcPosTrack->GetMother();
692 Int_t myTrackNegMotherLabel = mcNegTrack->GetMother();
693
694 if ((myTrackPosMotherLabel==-1)||(myTrackNegMotherLabel==-1)) continue;
695
696 AliAODMCParticle *mcPosMother = (AliAODMCParticle*)mcArray->At(myTrackPosMotherLabel);
697 if (!mcPosMother) continue;
698 Int_t MotherPdg = mcPosMother->GetPdgCode();
699 Int_t MotherOfMother = mcPosMother->GetMother();
700 //if (MotherOfMother == -1) MotherOfMotherPdg = 0;
701
702 if (myTrackPosMotherLabel!=myTrackNegMotherLabel) continue;
703 //if (!mcPosMother->IsPrimary()) continue;
704 Bool_t IsK0FromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==310)&&(PosTrackPdg==211)&&(NegTrackPdg==-211);
705 Bool_t IsLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
706 Bool_t IsAntiLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);
707
708 Bool_t ComeFromSigma = kFALSE;
709 Bool_t ComeFromSigmaLa = kFALSE;
710 Bool_t ComeFromSigmaAl = kFALSE;
711
712 if (MotherOfMother != -1)
713 {
714 AliAODMCParticle *mcPosMotherOfMother = (AliAODMCParticle*)mcArray->At(MotherOfMother);
715 MotherOfMotherPdg = mcPosMotherOfMother->GetPdgCode();
716 Int_t MoMPdg = TMath::Abs(MotherOfMotherPdg);
717 ComeFromSigma = (mcPosMotherOfMother->IsPhysicalPrimary())&&((MoMPdg==3212)||(MoMPdg==3224)||(MoMPdg==3214)||(MoMPdg==3114));
718 ComeFromSigmaLa = ComeFromSigma && (MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
719 ComeFromSigmaAl = ComeFromSigma && (MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);
720 }
721
722 IsLambdaFromMC = IsLambdaFromMC || ComeFromSigmaLa;
723 IsAntiLambdaFromMC = IsAntiLambdaFromMC || ComeFromSigmaAl;
724
725 Double_t RecMotherPt = aodV0->Pt();
726 //cout << "Pt of rec v0 = " << RecMotherPt << " nMC = " << mcArray->GetEntries() << endl;
727 //cout << "Pos Label = " << myTrackPosLabel << " Neg Label = " << myTrackNegLabel << endl;
728 Double_t rcK0[3] = {RecMotherPt, lCent, 1};
729 Double_t rcLa[3] = {RecMotherPt, lCent, 2};
730 Double_t rcAl[3] = {RecMotherPt, lCent, 3};
731 if ((cutK0sc)&&(K0Signal)&&(IsK0FromMC)) fHistRCPtCentTrig->Fill(rcK0);
732 if ((cutLambdasc)&&(LamSignal)&&(IsLambdaFromMC)) fHistRCPtCentTrig->Fill(rcLa);
733 if ((cutAntiLambdasc)&&(ALamSignal)&&(IsAntiLambdaFromMC)) fHistRCPtCentTrig->Fill(rcAl);
734
735 }
736
737 //===== End of the MC part for V0 reconstruction efficiency ===
738
739 Int_t nSelectedTracks = selectedTracks->GetEntries();
740 // Correlation part
741 //===================================
742 for (Int_t j = 0; j < nSelectedTracks; j++)
743 {
744 AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
745 if ((atr->Pt())>=(aodV0->Pt())) continue;
746 Double_t dEta = atr->Eta() - aodV0->Eta();
747 Double_t dPhi = atr->Phi() - aodV0->Phi();
748 if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
749 if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;
750 // removing autocorrelations
751 //----------------------------------
752 Int_t negID = myTrackNeg->GetID();
753 Int_t posID = myTrackPos->GetID();
754 Int_t atrID = atr->GetID();
755
756 if ((TMath::Abs(negID)+1)==(TMath::Abs(atrID))) continue;
757 if ((TMath::Abs(posID)+1)==(TMath::Abs(atrID))) continue;
758 //----------------------------------
759
760 fHistNTrigSib->Sumw2();
761 fHistdPhidEtaSib->Sumw2();
762 // Filling correlation histograms and histograms for triggers counting
763 //----------------- K0 ---------------------
764 if ((cutK0sc)&&(K0Signal))
765 {
766 Double_t spK0Sig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 1.};
767 fHistdPhidEtaSib->Fill(spK0Sig);
768 }
769
770 if ((cutK0sc)&&(K0Bckg))
771 {
772 Double_t spK0Bkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 4.};
773 fHistdPhidEtaSib->Fill(spK0Bkg);
774 }
775 //---------------- Lambda -------------------
776 if ((cutLambdasc)&&(LamSignal))
777 {
778 Double_t spLaSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 2.};
779 fHistdPhidEtaSib->Fill(spLaSig);
780 }
781
782 if ((cutLambdasc)&&(LamBckg))
783 {
784 Double_t spLaBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 5.};
785 fHistdPhidEtaSib->Fill(spLaBkg);
786 }
787 //------------- AntiLambda -------------------
788 if ((cutAntiLambdasc)&&(ALamSignal))
789 {
790 Double_t spAlSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 3.};
791 fHistdPhidEtaSib->Fill(spAlSig);
792 }
793
794 if ((cutAntiLambdasc)&&(ALamBckg))
795 {
796 Double_t spAlBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 6.};
797 fHistdPhidEtaSib->Fill(spAlBkg);
798 }
799
800 } // end of correlation loop
801 //===================================
802
803
804 } // end of V0 selection loop
805
806 // ===========================================
807 // Ch-Ch correlation part
808
809 if (ChChWith)
810 {
811 Int_t nSelectedChargedTriggers = selectedChargedTriggers->GetEntries();
812 for (Int_t i = 0; i < nSelectedChargedTriggers; i++)
813 {
814 AliV0ChBasicParticle* chTrig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
815 Double_t chTrigPhi = chTrig->Phi();
816 Double_t chTrigEta = chTrig->Eta();
817 Double_t chTrigPt = chTrig->Pt();
818
819 Double_t nTrigChSig[4] = {chTrigPt, lCent, lPVz, 7.};
820 fHistNTrigSib->Fill(nTrigChSig);
821
822 Int_t nSelectedTracks = selectedTracks->GetEntries();
823 for (Int_t j = 0; j < nSelectedTracks; j++)
824 {
825 AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
826 if ((atr->Pt())>=chTrigPt) continue;
827 Double_t dEta = atr->Eta() - chTrigEta;
828 Double_t dPhi = atr->Phi() - chTrigPhi;
829 if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
830 if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;
831
832 Double_t spCh[7] = {dPhi, dEta, chTrigPt, atr->Pt(), lCent, lPVz, 7.};
833 fHistdPhidEtaSib->Fill(spCh);
834 }
835 }
836 }
837 // end of Ch-Ch correlation part
838
839 // Mixing ==============================================
840
841 fHistdPhidEtaMix->Sumw2();
842 AliEventPool* pool = fPoolMgr->GetEventPool(lCent, lPVz);
843 if (!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", lCent, lPVz));
844 //pool->SetDebug(1);
845 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
846 {
847 Int_t nMix = pool->GetCurrentNEvents();
848 for (Int_t jMix=0; jMix<nMix; jMix++)
849 {// loop through mixing events
850 TObjArray* bgTracks = pool->GetEvent(jMix);
851 for (Int_t i=0; i<selectedV0s->GetEntriesFast(); i++)
852 {// loop through selected V0 particles
853 AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedV0s->At(i);
854 Double_t trigPhi = trig->Phi();
855 Double_t trigEta = trig->Eta();
856 Double_t trigPt = trig->Pt();
857 Short_t trigC = trig->WhichCandidate();
858 for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
859 { // mixing tracks loop
860 AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
861 // be careful tracks may have bigger pt than v0s.
862 if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;
863 Double_t dEtaMix = assoc->Eta() - trigEta;
864 Double_t dPhiMix = assoc->Phi() - trigPhi;
865 if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
866 if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;
867
868 Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
869 fHistdPhidEtaMix->Fill(spMix);
870
871 } // end of mixing track loop
872 }// end of loop through selected V0 particles
873 if (ChChWith)
874 {
875 for (Int_t i=0; i<selectedChargedTriggers->GetEntriesFast(); i++)
876 {// loop through selected charged trigger particles
877 AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
878 Double_t trigPhi = trig->Phi();
879 Double_t trigEta = trig->Eta();
880 Double_t trigPt = trig->Pt();
881 Short_t trigC = trig->WhichCandidate();
882 for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
883 { // mixing tracks loop
884 AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
885 // be careful tracks may have bigger pt than v0s.
886 if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;
887
888 Double_t dEtaMix = assoc->Eta() - trigEta;
889 Double_t dPhiMix = assoc->Phi() - trigPhi;
890 if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
891 if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;
892
893 Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
894 fHistdPhidEtaMix->Fill(spMix);
895
896 } // end of mixing track loop
897 }// end of loop through selected Ch particles
898 } // end of ChCh mixing
899
900 }// end of loop of mixing events
901 }
902
903 TObjArray* tracksClone = (TObjArray*) selectedTracks->Clone();
904 tracksClone->SetOwner(kTRUE);
905 pool->UpdatePool(tracksClone);
906 //delete selectedtracks;
907
908 PostData(1, fOutput);
909
910}
911//___________________________________________
912Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodPrimaryTrack(const AliAODTrack *t)
913{
914 // Pseudorapidity cut
915 if (TMath::Abs(t->Eta())>0.9) return kFALSE;
916 // Should correspond to set of cuts suitable for correlation analysis, 768 - hybrid tracks for 2011 data
917 if (!t->TestFilterBit(768)) return kFALSE;
918
919 return kTRUE;
920}
921//_____________________________________________
922Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodDaughterTrack(const AliAODTrack *t)
923{
924 // TPC refit
925 if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
926 // Minimum number of clusters
927 Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
928 if (nCrossedRowsTPC < 70) return kFALSE;
929 Int_t findable=t->GetTPCNclsF();
930 if (findable <= 0) return kFALSE;
931 if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
932
933 return kTRUE;
934}
935//______________________________________________
936Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodV0(const AliAODEvent* aod, const AliAODv0* aodV0, const AliAODTrack* myTrackPos, const AliAODTrack* myTrackNeg, Int_t oSta)
937{
938 if (!aodV0) {
939 AliError(Form("ERROR: Could not retrieve aodV0"));
940 return kFALSE;
941 }
942
943 // Offline reconstructed V0 only
944 if (oSta==1) {if (aodV0->GetOnFlyStatus()) return kFALSE;}
945 if (oSta==3) {if (!aodV0->GetOnFlyStatus()) return kFALSE;}
946
947 if (oSta==2)
948 {
949 if (aodV0->GetOnFlyStatus())
950 {
951 return kTRUE;
952 } else {
953 return kFALSE;
954 }
955 }
956 if (oSta==4)
957 {
958 if (!aodV0->GetOnFlyStatus())
959 {
960 return kTRUE;
961 } else {
962 return kFALSE;
963 }
964 }
965
966 // Get daughters and check them
967 AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
968 AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));
969
970 if (!myTrackPosTest || !myTrackNegTest) {
971 Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
972 return kFALSE;
973 }
974
975 if( myTrackPosTest->Charge() ==1){
976 myTrackPos = myTrackPosTest;
977 myTrackNeg = myTrackNegTest;
978 }
979
980 if( myTrackPosTest->Charge() ==-1){
981 myTrackPos = myTrackNegTest;
982 myTrackNeg = myTrackPosTest;
983 }
984
985 // Track cuts for daughter tracks
986 if ( !(IsMyGoodDaughterTrack(myTrackPos)) || !(IsMyGoodDaughterTrack(myTrackNeg)) ) return kFALSE;
987
988 // Unlike signs of daughters
989 if (myTrackNegTest->Charge() == myTrackPosTest->Charge()) return kFALSE;
990
991 // Rapidity cut
992 //Double_t lCutRap = 0.75;
993 //Double_t lRapK0s = aodV0->Y(310);
994 //Double_t lRapLambda = aodV0->Y(3122);
995 //Double_t lRapAntiLambda = aodV0->Y(-3122);
996
997 // Pseudorapidity cut - there are high pt V0
998 Double_t lCutEta = 0.75;
999 Double_t lEtaV0 = aodV0->Eta();
1000 if (TMath::Abs(lEtaV0)>=lCutEta) return kFALSE;
1001 //if (TMath::Abs(lRapK0s)>=lCutRap) return kFALSE;
1002 //if (TMath::Abs(lRapLambda)>=lCutRap) return kFALSE;
1003 //if (TMath::Abs(lRapAntiLambda)>=lCutRap) return kFALSE;
1004
1005 // getting global variables
1006 Float_t dcaDaughtersToPrimVtx = GetDcaDToPV();
1007 Float_t dcaBetweenDaughters = GetDcaV0D();
1008
1009 // DCA of daughter track to Primary Vertex
1010 Float_t xyn=aodV0->DcaNegToPrimVertex();
1011 if (TMath::Abs(xyn)<dcaDaughtersToPrimVtx) return kFALSE;
1012 Float_t xyp=aodV0->DcaPosToPrimVertex();
1013 if (TMath::Abs(xyp)<dcaDaughtersToPrimVtx) return kFALSE;
1014
1015 // DCA of daughter tracks
1016 Double_t dca=aodV0->DcaV0Daughters();
1017 if (dca>dcaBetweenDaughters) return kFALSE;
1018
1019 // Cosinus of pointing angle
1020 Double_t cpa=aodV0->CosPointingAngle(aod->GetPrimaryVertex());
1021 if (cpa<0.998) return kFALSE;
1022
1023 // Fiducial volume cut
1024 Double_t xyz[3]; aodV0->GetSecondaryVtx(xyz);
1025 Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
1026 if (r2<0.9*0.9) return kFALSE;
1027 if (r2>100*100) return kFALSE;
1028
1029 // c*tau cut - in main V0 loop - depends on particle hypothesis
1030
1031 // Minimum pt of daughters
1032 Double_t lMomPos[3] = {999,999,999};
1033 Double_t lMomNeg[3] = {999,999,999};
1034
1035 lMomPos[0] = aodV0->MomPosX();
1036 lMomPos[1] = aodV0->MomPosY();
1037 lMomPos[2] = aodV0->MomPosZ();
1038
1039 lMomNeg[0] = aodV0->MomNegX();
1040 lMomNeg[1] = aodV0->MomNegY();
1041 lMomNeg[2] = aodV0->MomNegZ();
1042
1043 Double_t lPtPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1]);
1044 Double_t lPtNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1]);
1045
1046 Double_t cutMinPtDaughter = 0.160;
1047 if (lPtPos<cutMinPtDaughter || lPtNeg<cutMinPtDaughter) return kFALSE;
1048
1049 // Daughter PID cut - in main V0 loop - depends on particle hypothesis
1050
1051 return kTRUE;
1052}
1053
1054void AliAnalysisTaskV0ChCorrelations::RemovingInjectedSignal(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
1055{
1056 // remove injected signals (primaries above <maxLabel>)
1057 // <tracks> can be the following cases:
1058 // a. tracks: in this case the label is taken and then case b.
1059 // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
1060 // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
1061
1062 TClonesArray* arrayMC = 0;
1063 AliMCEvent* mcEvent = 0;
1064 if (mcObj->InheritsFrom("AliMCEvent"))
1065 mcEvent = static_cast<AliMCEvent*>(mcObj);
1066 else if (mcObj->InheritsFrom("TClonesArray"))
1067 arrayMC = static_cast<TClonesArray*>(mcObj);
1068 else
1069 {
1070 mcObj->Dump();
1071 AliFatal("Invalid object passed");
1072 }
1073
1074 Int_t before = tracks->GetEntriesFast();
1075
1076 for (Int_t i=0; i<before; ++i)
1077 {
1078 AliVParticle* part = (AliVParticle*) tracks->At(i);
1079
1080 if (part->InheritsFrom("AliESDtrack")) part = mcEvent->GetTrack(TMath::Abs(part->GetLabel()));
1081 if (part->InheritsFrom("AliAODTrack"))
1082 {
1083 part =(AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel()));
1084 //cout << "toto musi byt len pri reco trackoch" << endl;
1085 }
1086
1087 AliVParticle* mother = part;
1088 if (mcEvent)
1089 {
1090 while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
1091 {
1092 if (((AliMCParticle*)mother)->GetMother() < 0)
1093 {
1094 mother = 0;
1095 break;
1096 }
1097
1098 mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
1099 if (!mother)
1100 break;
1101 }
1102 }
1103 else
1104 {
1105 // find the primary mother
1106 while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
1107 {
1108 if (((AliAODMCParticle*)mother)->GetMother() < 0)
1109 {
1110 mother = 0;
1111 break;
1112 }
1113
1114 mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1115 if (!mother)
1116 break;
1117 }
1118 }
1119
1120 if (!mother)
1121 {
1122 //Printf("WARNING: No mother found for particle %d:", part->GetLabel());
1123 continue;
1124 }
1125
1126 if (mother->GetLabel() > maxLabel)
1127 {
1128// Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
1129 TObject* object = tracks->RemoveAt(i);
1130 if (tracks->IsOwner())
1131 delete object;
1132 }
1133 }
1134
1135 tracks->Compress();
1136
1137 AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
1138
1139}
1140