]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliAnalysisTaskFlowStrange.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add setters to change azimuthal acceltance when semi-good tpc runs have been detecte...
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskFlowStrange.cxx
CommitLineData
24373b38 1/*************************************************************************
2* Copyright(c) 1998-2008,ALICE Experiment at CERN,All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use,copy,modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee,provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
16/////////////////////////////////////////////////////
17// AliAnalysisTaskFlowStrange:
18// Analysis task to select K0/Lambda candidates for flow analysis.
19// Authors: Cristian Ivan (civan@cern.ch)
20// Carlos Perez (cperez@cern.ch)
21// Pawel Debski (pdebski@cern.ch)
22//////////////////////////////////////////////////////
23
24#include "TChain.h"
25#include "TList.h"
26#include "TH1D.h"
27#include "TH2D.h"
28#include "TH3D.h"
29#include "TF1.h"
30#include "TProfile.h"
31#include "TProfile2D.h"
32#include "TVector3.h"
33#include "TStopwatch.h"
34#include "TFile.h"
35
36#include "TRandom3.h"
37
38#include "AliAnalysisManager.h"
39#include "AliInputEventHandler.h"
40
41#include "AliVVertex.h"
42#include "AliVVZERO.h"
43#include "AliStack.h"
44#include "AliMCEvent.h"
45
46#include "AliESDEvent.h"
47#include "AliESDtrack.h"
48#include "AliESDVertex.h"
49#include "AliESDv0.h"
50#include "AliESDtrackCuts.h"
51
52#include "AliAODEvent.h"
53#include "AliAODTrack.h"
54#include "AliAODVertex.h"
55#include "AliAODv0.h"
56#include "AliAODTracklets.h"
57#include "AliAODHeader.h"
58
59#include "AliAODMCHeader.h"
60#include "AliAODMCParticle.h"
61#include "TClonesArray.h"
62#include "TDatabasePDG.h"
63#include "TParticlePDG.h"
64
65#include "TMath.h"
66#include "TObjArray.h"
67#include "AliFlowCandidateTrack.h"
68
69#include "AliFlowTrackCuts.h"
70#include "AliFlowEventCuts.h"
71#include "AliFlowEvent.h"
72#include "AliFlowBayesianPID.h"
73#include "AliFlowCommonConstants.h"
74#include "AliFlowVector.h"
75
76#include "AliAnalysisTaskFlowStrange.h"
77
78ClassImp(AliAnalysisTaskFlowStrange)
79
80//=======================================================================
81AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange() :
82 AliAnalysisTaskSE(),
83 fPIDResponse(NULL),
84 fFB1(NULL),
85 fFB1024(NULL),
86 fTPCevent(NULL),
87 fVZEevent(NULL),
88 fCandidates(NULL),
89 fList(NULL),
90 fRunNumber(-1),
91 fDebug(0),
92 fQAlevel(0),
93 fReadESD(kFALSE),
94 fReadMC(kFALSE),
6c0d2d67 95 fAddPiToMCReactionPlane(kTRUE),
512ced40 96 fPostMatched(0),
24373b38 97 fAvoidExec(kFALSE),
98 fSkipSelection(kFALSE),
24373b38 99 fUseFP(kFALSE),
100 fRunOnpA(kFALSE),
101 fRunOnpp(kFALSE),
102 fExtraEventRejection(kFALSE),
6c0d2d67 103 fSkipCentralitySelection(kFALSE),
24373b38 104 fCentMethod("V0MTRK"),
105 fCentPerMin(0),
106 fCentPerMax(100),
107 fThisCent(-1.0),
6c0d2d67 108 fV0M(0.0),
109 fTRK(0.0),
110 fPriVtxZ(0.0),
111 fSPDVtxZ(0.0),
112 fSPDtracklets(0),
113 fVZETotM(0.0),
114 fRefMultTPC(0),
115 fRefMultHyb(0),
512ced40 116 fVertexZcut(10.0),
24373b38 117 fExcludeTPCEdges(kFALSE),
118 fSpecie(0),
119 fOnline(kFALSE),
120 fHomemade(kFALSE),
121 fWhichPsi(1),
122 fVZEsave(kFALSE),
123 fVZEload(NULL),
124 fVZEResponse(NULL),
125 fVZEmb(kFALSE),
126 fVZEByDisk(kTRUE),
127 fVZECa(0),
128 fVZECb(3),
129 fVZEAa(0),
130 fVZEAb(3),
131 fVZEQA(NULL),
6c0d2d67 132 fHarmonic(2),
24373b38 133 fPsi2(0.0),
134 fMCEP(0.0),
6c0d2d67 135 fQVZEACos(0.0),
136 fQVZEASin(0.0),
137 fQVZECCos(0.0),
138 fQVZECSin(0.0),
139 fQVZEA(0.0),
140 fQVZEC(0.0),
141 fQTPCACos(0.0),
142 fQTPCASin(0.0),
143 fQTPCCCos(0.0),
144 fQTPCCSin(0.0),
145 fQTPC2hCos(0.0),
146 fQTPC2hSin(0.0),
147 fQTPCA(0.0),
148 fQTPCC(0.0),
149 fQTPCA_nTracks(0),
150 fQTPCC_nTracks(0),
151 fSkipTerminate(kTRUE),
24373b38 152 fMassBins(0),
153 fMinMass(0.0),
154 fMaxMass(0.0),
512ced40
RAB		 155 fMinMassX(-1.0),
156 fMaxMassX(-1.0),
24373b38 157 fRFPFilterBit(1),
158 fRFPminPt(0.2),
159 fRFPmaxPt(5.0),
6c0d2d67 160 fRFPAminEta(0.0),
161 fRFPAmaxEta(+0.8),
162 fRFPCminEta(-0.8),
163 fRFPCmaxEta(0.0),
24373b38 164 fRFPTPCsignal(10.0),
165 fRFPmaxIPxy(2.4),
166 fRFPmaxIPz(3.2),
167 fRFPTPCncls(70),
168 fDecayMass(0.0),
169 fDecayPhi(0.0),
170 fDecayEta(0.0),
171 fDecayPt(0.0),
172 fDecayDCAdaughters(0.0),
173 fDecayCosinePointingAngleXY(0.0),
174 fDecayRadXY(0.0),
175 fDecayDecayLength(0.0),
176 fDecayQt(0.0),
177 fDecayAlpha(0.0),
178 fDecayRapidity(0.0),
179 fDecayProductIPXY(0.0),
512ced40
RAB		 180 fDecayIPneg(0.0),
181 fDecayIPpos(0.0),
182 fDecayXneg(0.0),
183 fDecayXpos(0.0),
24373b38 184 fDecayIDneg(-1),
185 fDecayIDpos(-1),
186 fDecayID(-1),
512ced40
RAB		 187 fDecayMatchOrigin(0.0),
188 fDecayMatchPhi(0.0),
189 fDecayMatchEta(0.0),
190 fDecayMatchPt(0.0),
191 fDecayMatchRadXY(0.0),
24373b38 192 fDecayMinEta(0.0),
193 fDecayMaxEta(0.0),
194 fDecayMinPt(0.0),
195 fDecayMaxDCAdaughters(0.0),
196 fDecayMinCosinePointingAngleXY(0.0),
197 fDecayMinQt(0.0),
198 fDecayAPCutPie(kTRUE),
199 fDecayMinRadXY(0.0),
200 fDecayMaxDecayLength(0.0),
201 fDecayMaxProductIPXY(0.0),
202 fDecayMaxRapidity(0.0),
203 fDaughterPhi(0.0),
204 fDaughterEta(0.0),
205 fDaughterPt(0.0),
206 fDaughterNClsTPC(0),
207 fDaughterCharge(0),
208 fDaughterNFClsTPC(0),
209 fDaughterNSClsTPC(0),
210 fDaughterChi2PerNClsTPC(0.0),
211 fDaughterXRows(0.0),
212 fDaughterImpactParameterXY(0.0),
213 fDaughterImpactParameterZ(0.0),
214 fDaughterStatus(0),
512ced40 215 fDaughterITScm(0),
24373b38 216 fDaughterNSigmaPID(0.0),
217 fDaughterKinkIndex(0),
6c0d2d67 218 fDaughterAtSecPhi(0.0),
219 fDaughterAtSecEta(0.0),
220 fDaughterAtSecPt(0.0),
512ced40
RAB		 221 fDaughterMatchPhi(0.0),
222 fDaughterMatchEta(0.0),
223 fDaughterMatchPt(0.0),
224 fDaughterMatchImpactParameterXY(0.0),
225 fDaughterMatchImpactParameterZ(0.0),
24373b38 226 fDaughterUnTag(kTRUE),
227 fDaughterMinEta(0.0),
228 fDaughterMaxEta(0.0),
229 fDaughterMinPt(0.0),
230 fDaughterMinNClsTPC(0),
231 fDaughterMinXRows(0),
232 fDaughterMaxChi2PerNClsTPC(0.0),
233 fDaughterMinXRowsOverNClsFTPC(0.0),
234 fDaughterMinImpactParameterXY(0.0),
235 fDaughterMaxNSigmaPID(0.0) {
236 //ctor
512ced40 237 for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
6c0d2d67 238 for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
239 for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
24373b38 240}
241//=======================================================================
242AliAnalysisTaskFlowStrange::AliAnalysisTaskFlowStrange(const char *name) :
243 AliAnalysisTaskSE(name),
244 fPIDResponse(NULL),
245 fFB1(NULL),
246 fFB1024(NULL),
247 fTPCevent(NULL),
248 fVZEevent(NULL),
249 fCandidates(NULL),
250 fList(NULL),
251 fRunNumber(-1),
252 fDebug(0),
253 fQAlevel(0),
254 fReadESD(kFALSE),
255 fReadMC(kFALSE),
6c0d2d67 256 fAddPiToMCReactionPlane(kTRUE),
512ced40 257 fPostMatched(0),
24373b38 258 fAvoidExec(kFALSE),
259 fSkipSelection(kFALSE),
24373b38 260 fUseFP(kFALSE),
261 fRunOnpA(kFALSE),
262 fRunOnpp(kFALSE),
263 fExtraEventRejection(kFALSE),
6c0d2d67 264 fSkipCentralitySelection(kFALSE),
24373b38 265 fCentMethod("V0MTRK"),
266 fCentPerMin(0),
267 fCentPerMax(100),
268 fThisCent(-1.0),
6c0d2d67 269 fV0M(0.0),
270 fTRK(0.0),
271 fPriVtxZ(0.0),
272 fSPDVtxZ(0.0),
273 fSPDtracklets(0),
274 fVZETotM(0.0),
275 fRefMultTPC(0),
276 fRefMultHyb(0),
512ced40 277 fVertexZcut(10.0),
24373b38 278 fExcludeTPCEdges(kFALSE),
279 fSpecie(0),
280 fOnline(kFALSE),
281 fHomemade(kFALSE),
282 fWhichPsi(1),
283 fVZEsave(kFALSE),
284 fVZEload(NULL),
285 fVZEResponse(NULL),
286 fVZEmb(kFALSE),
287 fVZEByDisk(kTRUE),
288 fVZECa(0),
289 fVZECb(3),
290 fVZEAa(0),
291 fVZEAb(3),
292 fVZEQA(NULL),
6c0d2d67 293 fHarmonic(2),
24373b38 294 fPsi2(0.0),
295 fMCEP(0.0),
6c0d2d67 296 fQVZEACos(0.0),
297 fQVZEASin(0.0),
298 fQVZECCos(0.0),
299 fQVZECSin(0.0),
300 fQVZEA(0.0),
301 fQVZEC(0.0),
302 fQTPCACos(0.0),
303 fQTPCASin(0.0),
304 fQTPCCCos(0.0),
305 fQTPCCSin(0.0),
306 fQTPC2hCos(0.0),
307 fQTPC2hSin(0.0),
308 fQTPCA(0.0),
309 fQTPCC(0.0),
310 fQTPCA_nTracks(0),
311 fQTPCC_nTracks(0),
312 fSkipTerminate(kTRUE),
24373b38 313 fMassBins(0),
314 fMinMass(0.0),
315 fMaxMass(0.0),
512ced40
RAB		 316 fMinMassX(-1.0),
317 fMaxMassX(-1.0),
24373b38 318 fRFPFilterBit(1),
319 fRFPminPt(0.2),
320 fRFPmaxPt(5.0),
6c0d2d67 321 fRFPAminEta(0.0),
322 fRFPAmaxEta(+0.8),
323 fRFPCminEta(-0.8),
324 fRFPCmaxEta(0.0),
24373b38 325 fRFPTPCsignal(10.0),
326 fRFPmaxIPxy(2.4),
327 fRFPmaxIPz(3.2),
328 fRFPTPCncls(70),
329 fDecayMass(0.0),
330 fDecayPhi(0.0),
331 fDecayEta(0.0),
332 fDecayPt(0.0),
333 fDecayDCAdaughters(0.0),
334 fDecayCosinePointingAngleXY(0.0),
335 fDecayRadXY(0.0),
336 fDecayDecayLength(0.0),
337 fDecayQt(0.0),
338 fDecayAlpha(0.0),
339 fDecayRapidity(0.0),
340 fDecayProductIPXY(0.0),
512ced40
RAB		 341 fDecayIPneg(0.0),
342 fDecayIPpos(0.0),
343 fDecayXneg(0.0),
344 fDecayXpos(0.0),
24373b38 345 fDecayIDneg(-1),
346 fDecayIDpos(-1),
347 fDecayID(-1),
512ced40
RAB		 348 fDecayMatchOrigin(0.0),
349 fDecayMatchPhi(0.0),
350 fDecayMatchEta(0.0),
351 fDecayMatchPt(0.0),
352 fDecayMatchRadXY(0.0),
24373b38 353 fDecayMinEta(0.0),
354 fDecayMaxEta(0.0),
355 fDecayMinPt(0.0),
356 fDecayMaxDCAdaughters(0.0),
357 fDecayMinCosinePointingAngleXY(0.0),
358 fDecayMinQt(0.0),
359 fDecayAPCutPie(kTRUE),
360 fDecayMinRadXY(0.0),
361 fDecayMaxDecayLength(0.0),
362 fDecayMaxProductIPXY(0.0),
363 fDecayMaxRapidity(0.0),
364 fDaughterPhi(0.0),
365 fDaughterEta(0.0),
366 fDaughterPt(0.0),
367 fDaughterNClsTPC(0),
368 fDaughterCharge(0),
369 fDaughterNFClsTPC(0),
370 fDaughterNSClsTPC(0),
371 fDaughterChi2PerNClsTPC(0.0),
372 fDaughterXRows(0.0),
373 fDaughterImpactParameterXY(0.0),
374 fDaughterImpactParameterZ(0.0),
375 fDaughterStatus(0),
512ced40 376 fDaughterITScm(0),
24373b38 377 fDaughterNSigmaPID(0.0),
378 fDaughterKinkIndex(0),
6c0d2d67 379 fDaughterAtSecPhi(0.0),
380 fDaughterAtSecEta(0.0),
381 fDaughterAtSecPt(0.0),
512ced40
RAB		 382 fDaughterMatchPhi(0.0),
383 fDaughterMatchEta(0.0),
384 fDaughterMatchPt(0.0),
385 fDaughterMatchImpactParameterXY(0.0),
386 fDaughterMatchImpactParameterZ(0.0),
24373b38 387 fDaughterUnTag(kTRUE),
388 fDaughterMinEta(0.0),
389 fDaughterMaxEta(0.0),
390 fDaughterMinPt(0.0),
391 fDaughterMinNClsTPC(0),
392 fDaughterMinXRows(0),
393 fDaughterMaxChi2PerNClsTPC(0.0),
394 fDaughterMinXRowsOverNClsFTPC(0.0),
395 fDaughterMinImpactParameterXY(0.0),
396 fDaughterMaxNSigmaPID(0.0) {
397 //ctor
512ced40 398 for(Int_t i=0; i!=6; ++i) fDaughterITSConfig[i]=-1;
6c0d2d67 399 for(Int_t i=0; i!=2000; ++i) fQTPCA_fID[i]=-1;
400 for(Int_t i=0; i!=2000; ++i) fQTPCC_fID[i]=-1;
24373b38 401 DefineInput( 0,TChain::Class());
402 DefineOutput(1,TList::Class());
403 DefineOutput(2,AliFlowEventSimple::Class()); // TPC object
404 DefineOutput(3,AliFlowEventSimple::Class()); // VZE object
405}
406//=======================================================================
407AliAnalysisTaskFlowStrange::~AliAnalysisTaskFlowStrange() {
408 //dtor
409 if (fCandidates) delete fCandidates;
410 if (fTPCevent) delete fTPCevent;
411 if (fVZEevent) delete fVZEevent;
412 if (fList) delete fList;
413}
414//=======================================================================
6c0d2d67 415TList* AliAnalysisTaskFlowStrange::RunTerminateAgain(TList *lst) {
416 fList = lst;
417 fSpecie = Int_t( ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSpecie) );
418 fSkipSelection = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kSkipSelection);
419 fReadMC = ((TProfile*)((TList*)fList->FindObject("Event"))->FindObject("Configuration"))->GetBinContent(kReadMC);
420 Terminate(NULL);
421 return fList;
422}
423//=======================================================================
24373b38 424void AliAnalysisTaskFlowStrange::UserCreateOutputObjects() {
425 //UserCreateOutputObjects
426 fList=new TList();
427 fList->SetOwner();
428 AddQAEvents();
429 AddQACandidates();
430
431 if(fReadESD) MakeFilterBits();
432
433 AliFlowCommonConstants *cc = AliFlowCommonConstants::GetMaster();
434 cc->SetNbinsMult(3000); cc->SetMultMin(0); cc->SetMultMax(30000);
6c0d2d67 435 cc->SetNbinsPt(100); cc->SetPtMin(0.0); cc->SetPtMax(20.0);
24373b38 436 cc->SetNbinsPhi(100); cc->SetPhiMin(0.0); cc->SetPhiMax(TMath::TwoPi());
437 cc->SetNbinsEta(100); cc->SetEtaMin(-5.0); cc->SetEtaMax(+5.0);
438 cc->SetNbinsQ(100); cc->SetQMin(0.0); cc->SetQMax(3.0);
439 cc->SetNbinsMass(fMassBins);
440 cc->SetMassMin(fMinMass);
441 cc->SetMassMax(fMaxMass);
442
512ced40
RAB		 443 if(fMinMassX<0) {
444 if(fSpecie==0) {
445 fMinMassX = 0.494;
446 fMaxMassX = 0.502;
447 // double lowEdge[14]={0.398, 0.420, 0.444, 0.468, 0.486,
448 // 0.490, 0.494, 0.498, 0.502, 0.506,
449 // 0.524, 0.548, 0.572, 0.598};
450 } else {
451 fMinMassX = 1.114;
452 fMaxMassX = 1.118;
453 // double lowEdge[13]={1.084, 1.094, 1.104, 1.110, 1.114,
454 // 1.116, 1.118, 1.122, 1.128, 1.138,
455 // 1.148, 1.158, 1.168};} else {
456 }
457 }
24373b38 458 //loading pid response
459 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
460 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
461 fPIDResponse = inputHandler->GetPIDResponse();
462
6c0d2d67 463 if(fUseFP) {
464 fTPCevent = new AliFlowEvent(100);
465 fVZEevent = new AliFlowEvent(100);
466 //array of candidates
467 fCandidates = new TObjArray(100);
468 fCandidates->SetOwner();
469 }
24373b38 470 PostData(1,fList);
471 if(fUseFP) { // for connection to the flow package
472 PostData(2,fTPCevent);
473 PostData(3,fVZEevent);
474 }
475
476 gRandom->SetSeed();
477}
478//=======================================================================
479void AliAnalysisTaskFlowStrange::AddQAEvents() {
480 // function to add event qa
481 TH1D *tH1D;
482 TProfile *tProfile;
483 TList *tQAEvents=new TList();
484 tQAEvents->SetName("Event");
485 tQAEvents->SetOwner();
486 tH1D = new TH1D("Events","Number of Events",6,0,6); tQAEvents->Add(tH1D);
487 tH1D->GetXaxis()->SetBinLabel(1,"exec");
488 tH1D->GetXaxis()->SetBinLabel(2,"userexec");
489 tH1D->GetXaxis()->SetBinLabel(3,"reached");
490 tH1D->GetXaxis()->SetBinLabel(4,"selected");
491 tH1D->GetXaxis()->SetBinLabel(5,"rejectedByLowQw");
492 tH1D->GetXaxis()->SetBinLabel(6,"rejectedByErrorLoadVZEcal");
6c0d2d67 493 tProfile = new TProfile("Configuration","Configuration",10,0.5,10.5); tQAEvents->Add(tProfile);
494 tProfile->Fill(kSpecie,fSpecie,1);
495 tProfile->GetXaxis()->SetBinLabel(kSpecie,"fSpecie");
496 tProfile->Fill(kHarmonic,fHarmonic,1);
497 tProfile->GetXaxis()->SetBinLabel(kHarmonic,"fHarmonic");
498 tProfile->Fill(kReadMC,fReadMC,1);
499 tProfile->GetXaxis()->SetBinLabel(kReadMC,"fReadMC");
500 tProfile->Fill(kSkipSelection,fSkipSelection,1);
501 tProfile->GetXaxis()->SetBinLabel(kSkipSelection,"fSkipSelection");
24373b38 502 tH1D = new TH1D("POI","POIs;multiplicity",800,0,800); tQAEvents->Add(tH1D);
503 tH1D = new TH1D("UNTAG","UNTAG;Untagged Daughters",800,0,800);tQAEvents->Add(tH1D);
504 tH1D = new TH1D("RealTime","RealTime;LogT sec",2000,-10,+10); tQAEvents->Add(tH1D);
505 fList->Add(tQAEvents);
6c0d2d67 506 AddEventSpy("EventsRaw");
c41a93af 507 AddEventSpy("EventsReached");
508 AddEventSpy("EventsSelected");
24373b38 509 AddMakeQSpy();
510}
511//=======================================================================
c41a93af 512void AliAnalysisTaskFlowStrange::AddEventSpy(TString name) {
24373b38 513 TH1D *tH1D;
514 TH2D *tH2D;
515 TList *tList=new TList();
c41a93af 516 tList->SetName(name.Data());
24373b38 517 tList->SetOwner();
6c0d2d67 518 tH2D = new TH2D("VTXZ","VTXZ;PriVtxZ;SPDVtxZ",60,-25,+25,60,-25,+25); tList->Add( tH2D );
c41a93af 519 tH2D = new TH2D("CCCC","CCCC;V0M;TRK",60,-10,110,60,-10,110); tList->Add( tH2D );
6c0d2d67 520 tH2D = new TH2D("REFM","REFM;TPC ONLY;HYBRID",100,0,3000,100,0,3000); tList->Add( tH2D );
521 tH2D = new TH2D("SPDVZE","SPDVZE;SPD Tracklets;Total Multiplicity in VZERO",100,0,3500,100,0,25000); tList->Add( tH2D );
c41a93af 522 if(fReadMC) {
c41a93af 523 tH1D = new TH1D("MCEP","MCEP;MCEP",100,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
24373b38 524 }
525 fList->Add(tList);
526}
527//=======================================================================
6c0d2d67 528void AliAnalysisTaskFlowStrange::FillEventSpy(TString name) {
529 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("VTXZ"))->Fill( fPriVtxZ, fSPDVtxZ );
530 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("CCCC"))->Fill( fV0M, fTRK );
531 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("REFM"))->Fill( fRefMultTPC, fRefMultHyb );
532 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject("SPDVZE"))->Fill( fSPDtracklets, fVZETotM );
533 if(fReadMC) {
534 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject("MCEP"))->Fill( fMCEP );
535 }
536}
537//=======================================================================
24373b38 538void AliAnalysisTaskFlowStrange::AddMakeQSpy() {
24373b38 539 TH1D *tH1D;
540 TH2D *tH2D;
541 TProfile *tPF1;
542 TList *tList=new TList();
543 tList->SetName("MakeQSpy");
544 tList->SetOwner();
6c0d2d67 545 fList->Add(tList);
24373b38 546 tH1D = new TH1D("RFPTPC","TPC Refrence Multiplicity;multiplicity",3000,0,3000); tList->Add( tH1D );
547 tH1D = new TH1D("RFPVZE","VZERO Reference Multiplicity;multiplicity",3000,0,30000); tList->Add( tH1D );
512ced40
RAB		 548 tH1D = new TH1D("QmTPC","TPC Normalized Q vector;|Q|/M",3000,0,1); tList->Add( tH1D );
549 tH1D = new TH1D("QmVZE","VZERO Normalized Q vector;|Q|/M",3000,0,1); tList->Add( tH1D );
24373b38 550 tH2D = new TH2D("TPCAllPhiEta","TPCall;Phi;Eta",180,0,TMath::TwoPi(),80,-0.9,+0.9); tList->Add( tH2D );
551 tH2D = new TH2D("VZEAllPhiEta","VZEall;Phi;Eta",20,0,TMath::TwoPi(),40,-4.0,+6.0); tList->Add( tH2D );
552 tH1D = new TH1D("TPCPSI","TPCPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
553 tH1D = new TH1D("TPCPSIA","TPCPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
6c0d2d67 554 tH1D = new TH1D("TPCPSIC","TPCPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
24373b38 555 tH1D = new TH1D("VZEPSI","VZEPSI;PSI",72,0,TMath::Pi()); tList->Add( tH1D );
556 tH1D = new TH1D("VZEPSIA","VZEPSIA;PSIA",72,0,TMath::Pi()); tList->Add( tH1D );
6c0d2d67 557 tH1D = new TH1D("VZEPSIC","VZEPSIC;PSIC",72,0,TMath::Pi()); tList->Add( tH1D );
558 tPF1 = new TProfile("TPCQm","TPCQm",6,0.5,6.5); tList->Add( tPF1 );
559 tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
560 tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
24373b38 561 tPF1->GetXaxis()->SetBinLabel(5,"Qy"); tPF1->GetXaxis()->SetBinLabel(6,"Qx");
6c0d2d67 562 tPF1 = new TProfile("VZEQm","VZEQm",6,0.5,6.5); tList->Add( tPF1 );
563 tPF1->GetXaxis()->SetBinLabel(1,"Qcy"); tPF1->GetXaxis()->SetBinLabel(2,"Qcx");
564 tPF1->GetXaxis()->SetBinLabel(3,"Qay"); tPF1->GetXaxis()->SetBinLabel(4,"Qax");
24373b38 565 tPF1->GetXaxis()->SetBinLabel(5,"Qy"); tPF1->GetXaxis()->SetBinLabel(6,"Qx");
6c0d2d67 566 tPF1 = new TProfile("QmVZEAQmVZEC","QmVZEAQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
567 tPF1 = new TProfile("QmVZEASQUARED","QmVZEASQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
568 tPF1 = new TProfile("QmVZECSQUARED","QmVZECSQUARED",1,0.5,1.5,"s"); tList->Add( tPF1 );
569 tPF1 = new TProfile("QmTPCQmVZEA","QmTPCQmVZEA",1,0.5,1.5,"s"); tList->Add( tPF1 );
570 tPF1 = new TProfile("QmTPCQmVZEC","QmTPCQmVZEC",1,0.5,1.5,"s"); tList->Add( tPF1 );
571 tH1D = new TH1D("ChiSquaredVZEA","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
572 tH1D = new TH1D("ChiSquaredVZEC","ChiSquaredVZEC",1,0.5,1.5); tList->Add( tH1D );
573 if(fReadMC) {
574 tH1D = new TH1D("PSIMCDIFFTPC","PSIMCDIFFTPC;MC-TPC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
575 tH1D = new TH1D("PSIMCDIFFTPCA","PSIMCDIFFTPCA;MC-TPCA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
576 tH1D = new TH1D("PSIMCDIFFTPCC","PSIMCDIFFTPCC;MC-TPCC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
577 tH1D = new TH1D("PSIMCDIFFVZE","PSIMCDIFFVZE;MC-VZE",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
578 tH1D = new TH1D("PSIMCDIFFVZEA","PSIMCDIFFVZEA;MC-VZEA",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
579 tH1D = new TH1D("PSIMCDIFFVZEC","PSIMCDIFFVZEC;MC-VZEC",72,-TMath::TwoPi(),TMath::TwoPi()); tList->Add( tH1D );
24373b38 580 }
6c0d2d67 581 tList=new TList(); tList->SetName("TPCRFPall"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
582 tList=new TList(); tList->SetName("TPCRFPsel"); tList->SetOwner(); AddTPCRFPSpy(tList); fList->Add(tList);
24373b38 583}
584//=======================================================================
585void AliAnalysisTaskFlowStrange::AddQACandidates() {
586 // function to add histogramming for candidates
587 if(fSkipSelection) return;
24373b38 588 TList *tList;
589 TH1D *tH1D;
590 TH2D *tH2D;
24373b38 591
592 //reconstruction
593 if(fReadESD) {
6c0d2d67 594 tList=new TList(); tList->SetName("ESD_TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
595 tList=new TList(); tList->SetName("ESD_TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
24373b38 596 tH2D = new TH2D("NPAIR", "NPAIR;NPOS;NNEG",1000,0,5000,1000,0,5000); tList->Add(tH2D);
597 tH2D = new TH2D("PtIPXY","PtIPXY;Pt;IPxy", 100,0,10,200,-10,+10); tList->Add(tH2D);
598 }
6c0d2d67 599 //aod prefilter
600 if(fSpecie<10) {
601 //candidates
602 tList=new TList(); tList->SetName("V0SAll"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
603 tH2D = new TH2D("V0SADC","V0S AFTER DAUGHTER CUTS;V0ALL;V0IMW",100,0,1000,100,0,1000); tList->Add(tH2D);
604 tList=new TList(); tList->SetName("V0SSel"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
605 //daughters
606 tList=new TList(); tList->SetName("V0SDau"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
607 //flow
608 tList=new TList(); tList->SetName("V0SAllVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
609 tList=new TList(); tList->SetName("V0SSelVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
24373b38 610 // IN-OUT
6c0d2d67 611 if(fQAlevel>1) {
612 tList=new TList(); tList->SetName("V0SAllIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
613 tList=new TList(); tList->SetName("V0SAllOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
614 tList=new TList(); tList->SetName("V0SSelIP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
615 tList=new TList(); tList->SetName("V0SSelOP"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
616 }
617 //match
618 if(fReadMC) {
619 tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
620 tH1D = new TH1D("Events", "Events",5,0.5,5.5); tList->Add(tH1D);
621 tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
622 tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
623 tH1D->GetXaxis()->SetBinLabel(3,"Daughters in stack");
624 tH1D->GetXaxis()->SetBinLabel(4,"Correspond to decay");
625 tH1D->GetXaxis()->SetBinLabel(5,"Decay has mother");
626 tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
627 tList=new TList(); tList->SetName("MthDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
628 tList=new TList(); tList->SetName("MthPosPos"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
629 tList=new TList(); tList->SetName("MthNegNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
630 tList=new TList(); tList->SetName("MthPosNeg"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
631 tList=new TList(); tList->SetName("MthNegDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
632 tList=new TList(); tList->SetName("MthPosDau"); tList->SetOwner(); AddTrackSpy(tList,true); fList->Add(tList);
633 tList=new TList(); tList->SetName("MthFeedDown"); tList->SetOwner(); AddCandidatesSpy(tList,true); fList->Add(tList);
634 tList=new TList(); tList->SetName("UnMth"); tList->SetOwner(); AddCandidatesSpy(tList,false); fList->Add(tList);
635 tList=new TList(); tList->SetName("UnMthDau"); tList->SetOwner(); AddTrackSpy(tList,false); fList->Add(tList);
636 tList=new TList(); tList->SetName("V0SMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
637 tList=new TList(); tList->SetName("V0SMthPosPosVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
638 tList=new TList(); tList->SetName("V0SMthNegNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
639 tList=new TList(); tList->SetName("V0SMthPosNegVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
640 tList=new TList(); tList->SetName("V0SUnMthVn"); tList->SetOwner(); AddDecayVn(tList); fList->Add(tList);
641 }
642 } else {
643 //candidates
644 tList=new TList(); tList->SetName("TrkAll"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
645 tList=new TList(); tList->SetName("TrkSel"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
646 tList=new TList(); tList->SetName("TrkAllVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
647 tList=new TList(); tList->SetName("TrkSelVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
648 //match
649 if(fReadMC) {
650 tList=new TList(); tList->SetName("STATMC"); tList->SetOwner(); fList->Add(tList);
651 tH1D = new TH1D("Events", "Events",3,0.5,3.5); tList->Add(tH1D);
652 tH1D->GetXaxis()->SetBinLabel(1,"Selected events");
653 tH1D->GetXaxis()->SetBinLabel(2,"Stack found");
654 tH1D->GetXaxis()->SetBinLabel(3,"Track in stack");
655 tList=new TList(); tList->SetName("Mth"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
656 tList=new TList(); tList->SetName("MthPos"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
657 tList=new TList(); tList->SetName("MthNeg"); tList->SetOwner(); AddTrackSpy(tList); fList->Add(tList);
658 tList=new TList(); tList->SetName("MthVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
659 tList=new TList(); tList->SetName("MthPosVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
660 tList=new TList(); tList->SetName("MthNegVn"); tList->SetOwner(); AddTrackVn(tList); fList->Add(tList);
661 }
24373b38 662 }
663 //stack
664 if(fReadMC) {
6c0d2d67 665 //tList=new TList(); tList->SetName("GenTru"); tList->SetOwner(); AddCandidatesSpy(tList); fList->Add(tList);
666 tList=new TList(); tList->SetName("MCTPionGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
667 tList=new TList(); tList->SetName("MCTKaonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
24373b38 668 tList=new TList(); tList->SetName("MCTK0sGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
6c0d2d67 669 tList=new TList(); tList->SetName("MCTProtonGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
24373b38 670 tList=new TList(); tList->SetName("MCTLdaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
671 tList=new TList(); tList->SetName("MCTPhiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
672 tList=new TList(); tList->SetName("MCTXiGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
6c0d2d67 673 tList=new TList(); tList->SetName("MCTOmegaGenAcc"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
674 tList=new TList(); tList->SetName("MCTPion"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
675 tList=new TList(); tList->SetName("MCTKaon"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
24373b38 676 tList=new TList(); tList->SetName("MCTK0s"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
677 tList=new TList(); tList->SetName("MCTLda"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
6c0d2d67 678 tList=new TList(); tList->SetName("MCTProton"); tList->SetOwner(); AddMCParticleSpy(tList); fList->Add(tList);
24373b38 679 }
680}
681//=======================================================================
682void AliAnalysisTaskFlowStrange::Exec(Option_t* option) {
683 // bypassing ::exec (needed because of AMPT)
684 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(0);
685 if(fAvoidExec) {
686 AliAnalysisTaskFlowStrange::UserExec(option);
687 } else {
688 AliAnalysisTaskSE::Exec(option);
689 }
690}
691//=======================================================================
692void AliAnalysisTaskFlowStrange::UserExec(Option_t *option) {
693 // bridge
694 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(1);
695 AliAnalysisTaskFlowStrange::MyUserExec(option);
696}
697//=======================================================================
698void AliAnalysisTaskFlowStrange::MyNotifyRun() {
699 if(fQAlevel>5 && !fReadESD) AddVZEQA();
700 if(fVZEsave) AddVZEROResponse();
701}
702//=======================================================================
703Bool_t AliAnalysisTaskFlowStrange::CalibrateEvent() {
704 if(fVZEsave) SaveVZEROResponse();
705 if(fQAlevel>5 && !fReadESD) SaveVZEROQA(); // 2BIMPROVED
706 Bool_t okay=kTRUE;
707 if(fVZEload) {
708 LoadVZEROResponse();
709 if(!fVZEResponse) okay = kFALSE;
710 }
711 return okay;
712}
713//=======================================================================
6c0d2d67 714Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliESDEvent*) {
715 // ESD reading discontinued: TO BE UPDATED
716 /*
24373b38 717 Double_t acceptEvent=kTRUE;
718 Double_t tTPCVtxZ = tESD->GetPrimaryVertexTPC()->GetZ();
719 if(tESD->GetPrimaryVertexTPC()->GetNContributors()<=0) return kFALSE;
720 Double_t tSPDVtxZ = tESD->GetPrimaryVertexSPD()->GetZ();
721 if(tESD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
722 // EventCuts
723 AliCentrality *cent = tESD->GetCentrality();
724 Double_t cc1, cc2;
725 cc1 = cent->GetCentralityPercentile("V0M");
726 cc2 = cent->GetCentralityPercentile("TRK");
727 TString mycent = fCentMethod;
728 if(fCentMethod.Contains("V0MTRK")) {
729 acceptEvent = TMath::Abs(cc1-cc2)>5.0?kFALSE:acceptEvent; // a la Alex
730 mycent = "V0M";
731 }
732 fThisCent = cent->GetCentralityPercentile( mycent );
733 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
734 acceptEvent = TMath::Abs(tTPCVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
512ced40 735 acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
c41a93af 736 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
737 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
24373b38 738 // EndOfCuts
c41a93af 739 if(acceptEvent) {
740 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
741 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
742 }
24373b38 743 return acceptEvent;
6c0d2d67 744 */
745 return kFALSE;
24373b38 746}
747//=======================================================================
6c0d2d67 748Bool_t AliAnalysisTaskFlowStrange::MinimumRequirementsAA(AliAODEvent *tAOD) {
749 fRunNumber = tAOD->GetRunNumber();
24373b38 750 AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
6c0d2d67 751 fV0M = cent->GetCentralityPercentile("V0M");
752 fTRK = cent->GetCentralityPercentile("TRK");
24373b38 753 TString mycent = fCentMethod;
754 if(fCentMethod.Contains("V0MTRK")) {
24373b38 755 mycent = "V0M";
756 }
757 fThisCent = cent->GetCentralityPercentile( mycent );
6c0d2d67 758 fPriVtxZ = tAOD->GetPrimaryVertex()->GetZ();
759 fSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
760 fSPDtracklets = tAOD->GetTracklets()->GetNumberOfTracklets();
761 fVZETotM = tAOD->GetVZEROData()->GetMTotV0A() + tAOD->GetVZEROData()->GetMTotV0C();
762 int hyb_fb = 272; // for 2010h::AOD086
763 if(fRunNumber>=166529&&fRunNumber<=170593) {
764 hyb_fb = 768; // for 2011h::AOD145
765 }
766 fRefMultTPC = RefMult(tAOD,128);
767 fRefMultHyb = RefMult(tAOD,hyb_fb);
24373b38 768 if(fReadMC) {
6c0d2d67 769 fMCEP = -999;
24373b38 770 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(tAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
6c0d2d67 771 if(mcHeader) {
772 fMCEP = mcHeader->GetReactionPlaneAngle();
773 if(fAddPiToMCReactionPlane) fMCEP += (gRandom->Rndm()>0.5)*TMath::Pi();
24373b38 774 }
24373b38 775 }
6c0d2d67 776 // centrality selection health
777 // cut in Vtx 10 & NContributors
778 if(!fSkipCentralitySelection) if(fThisCent<0||fThisCent>100) return kFALSE;
779 // vtx z position compatibility within 5 mm
780 if(TMath::Abs(fPriVtxZ-fSPDVtxZ)>0.5) return kFALSE;
781 // specific cuts for 2010h
782 if(fRunNumber>=136851&&fRunNumber<=139517) {
783 //if( tpc>400.0+26.0/15.0*hyb ) return kFALSE;
784 //if( mvze>3000.0+22.0/3.0*ntrklets ) return kFALSE;
785 }
786 // specific cuts for 2011h
787 if(fRunNumber>=166529&&fRunNumber<=170593) {
788 //if( tpc>400.0+26.0/15.0*hyb ) return kFALSE;
789 //if( mvze>3000.0+22.0/3.0*ntrklets ) return kFALSE;
c41a93af 790 }
6c0d2d67 791 return kTRUE;
792}
793//=======================================================================
794Bool_t AliAnalysisTaskFlowStrange::AcceptAAEvent(AliAODEvent *tAOD) {
795 Bool_t minimum = MinimumRequirementsAA(tAOD);
796 FillEventSpy("EventsRaw");
797 if(!minimum) return kFALSE;
798
799 Double_t acceptEvent=kTRUE;
800 TString mycent = fCentMethod;
801 if(fCentMethod.Contains("V0MTRK")) {
802 acceptEvent = TMath::Abs(fV0M-fTRK)>5.0?kFALSE:acceptEvent;
803 mycent = "V0M";
24373b38 804 }
6c0d2d67 805 if(!fSkipCentralitySelection) acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
806 acceptEvent = TMath::Abs(fPriVtxZ)>fVertexZcut?kFALSE:acceptEvent;
807 // HISTOGRAMMING
808 FillEventSpy("EventsReached");
809 if(acceptEvent) FillEventSpy("EventsSelected");
24373b38 810 return acceptEvent;
811}
812//=======================================================================
6c0d2d67 813Bool_t AliAnalysisTaskFlowStrange::AcceptPPEvent(AliAODEvent*) {
814 // PP reading discontinued: TO BE UPDATED
815 /*
24373b38 816 Double_t acceptEvent=kTRUE;
817 Double_t tVtxZ = tAOD->GetPrimaryVertex()->GetZ();
818 if(tAOD->GetPrimaryVertex()->GetNContributors()<=0) return kFALSE;
819 Double_t tSPDVtxZ = tAOD->GetPrimaryVertexSPD()->GetZ();
820 if(tAOD->GetPrimaryVertexSPD()->GetNContributors()<=0) return kFALSE;
821 // EventCuts
822 AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
823 Double_t cc1, cc2;
824 cc1 = cent->GetCentralityPercentile("V0M");
825 cc2 = cent->GetCentralityPercentile("TRK");
826 fThisCent = GetReferenceMultiplicity();
827 //for pp i use fCentPerXXX to select on multiplicity
828 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
829 acceptEvent = TMath::Abs(tVtxZ-tSPDVtxZ)>0.5?kFALSE:acceptEvent;
512ced40 830 acceptEvent = TMath::Abs(tVtxZ)>fVertexZcut?kFALSE:acceptEvent;
c41a93af 831 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
832 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
24373b38 833 // EndOfCuts
c41a93af 834 if(acceptEvent) {
835 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tVtxZ, tSPDVtxZ );
836 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
837 }
24373b38 838 return acceptEvent;
6c0d2d67 839 */
840 return kFALSE;
24373b38 841}
842//=======================================================================
843Int_t AliAnalysisTaskFlowStrange::GetReferenceMultiplicity() { //toberefined
844 AliAODEvent *tAOD = (AliAODEvent *) InputEvent();
845 if(!tAOD) return -1;
846 AliAODTrack *track;
847 Int_t rawN = tAOD->GetNumberOfTracks();
848 Int_t ref=0;
849 for(Int_t id=0; id!=rawN; ++id) {
850 track = tAOD->GetTrack(id);
851 if(!track->TestFilterBit(fRFPFilterBit)) continue;
852 ++ref;
853 }
854 return ref;
855}
856//=======================================================================
6c0d2d67 857Bool_t AliAnalysisTaskFlowStrange::AcceptPAEvent(AliAODEvent*) {
858 // PA reading discontinued: TO BE UPDATED
859 /*
24373b38 860 //if(aod->GetHeader()->GetEventNumberESDFile() == 0) return; //rejecting first chunk NOT NEEDED ANYMORE
861 Int_t bc2 = tAOD->GetHeader()->GetIRInt2ClosestInteractionMap();
862 if(bc2!=0) return kFALSE;
863 Int_t bc1 = tAOD->GetHeader()->GetIRInt1ClosestInteractionMap();
864 if(bc1!=0) return kFALSE;
865 Short_t isPileup = tAOD->IsPileupFromSPD(5);
866 if(isPileup!=0) return kFALSE;
867 if(tAOD->GetHeader()->GetRefMultiplicityComb08()<0) return kFALSE;
868
869 const AliAODVertex* spdVtx = tAOD->GetPrimaryVertexSPD();
870 if(!spdVtx) return kFALSE;
871 if(spdVtx->GetNContributors()<=0) return kFALSE;
872
873 const AliAODVertex* tpcVtx=NULL;
874 Int_t nVertices = tAOD->GetNumberOfVertices();
875 for(Int_t iVertices = 0; iVertices < nVertices; iVertices++){
876 const AliAODVertex* vertex = tAOD->GetVertex(iVertices);
877 if (vertex->GetType() != AliAODVertex::kMainTPC) continue;
878 tpcVtx = vertex;
879 }
880 if(!tpcVtx) return kFALSE;
881 if(tpcVtx->GetNContributors()<=0) return kFALSE;
882 Double_t tTPCVtxZ = tpcVtx->GetZ();
883 Double_t tSPDVtxZ = spdVtx->GetZ();
884 if (TMath::Abs(tSPDVtxZ - tTPCVtxZ)>2.0) return kFALSE;
885 if(plpMV(tAOD)) return kFALSE;
886
887 Double_t acceptEvent=kTRUE;
888 // EventCuts
889 AliCentrality *cent = tAOD->GetHeader()->GetCentralityP();
890 Double_t cc1, cc2;
891 cc1 = cent->GetCentralityPercentile("V0M");
892 cc2 = cent->GetCentralityPercentile("TRK");
893 if(fCentMethod.Contains("V0MTRK")) fCentMethod = "V0M";
894 fThisCent = cent->GetCentralityPercentile( fCentMethod );
895 acceptEvent = (fThisCent<fCentPerMin||fThisCent>fCentPerMax)?kFALSE:acceptEvent;
512ced40 896 acceptEvent = TMath::Abs(tTPCVtxZ)>fVertexZcut?kFALSE:acceptEvent;
24373b38 897 // EndOfCuts
c41a93af 898 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
899 ((TH2D*)((TList*)fList->FindObject("EventsReached"))->FindObject("CCCC"))->Fill( cc1, cc2 );
900 if(acceptEvent) {
901 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("VTXZ"))->Fill( tTPCVtxZ, tSPDVtxZ );
902 ((TH2D*)((TList*)fList->FindObject("EventsSelected"))->FindObject("CCCC"))->Fill( cc1, cc2 );
903 }
24373b38 904 return acceptEvent;
6c0d2d67 905 */
906 return kFALSE;
24373b38 907}
908//=======================================================================
909void AliAnalysisTaskFlowStrange::MyUserExec(Option_t *) {
910 // MAIN ROUTINE
911 TStopwatch tTime;
912 tTime.Start();
6c0d2d67 913 if(fDebug) {
914 printf("****************\n");
915 printf("****************\n");
916 printf("**::MyUserExec()\n");
917 }
918 if(fUseFP) fCandidates->SetLast(-1);
24373b38 919 AliESDEvent *tESD=dynamic_cast<AliESDEvent*>(InputEvent());
920 AliAODEvent *tAOD=dynamic_cast<AliAODEvent*>(InputEvent());
6c0d2d67 921 Int_t prevRun = fRunNumber;
24373b38 922 //=>check event
923 Bool_t acceptEvent=kFALSE;
924 if(fReadESD) {
925 if(!tESD) {ResetContainers(); Publish(); return;}
926 acceptEvent = fRunOnpp?kFALSE:fRunOnpA?kFALSE:AcceptAAEvent(tESD);
24373b38 927 } else {
928 if(!tAOD) {ResetContainers(); Publish(); return;}
929 acceptEvent = fRunOnpp?AcceptPPEvent(tAOD):fRunOnpA?AcceptPAEvent(tAOD):AcceptAAEvent(tAOD);
24373b38 930 }
6c0d2d67 931 if(prevRun!=fRunNumber) {
24373b38 932 MyNotifyRun();
933 }
934 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(2);
935 //=>does the event clear?
936 if(!acceptEvent) {ResetContainers(); Publish(); return;}
937 // healthy event incomming
938 if( !CalibrateEvent() ) { // saves/retrieves/qas VZEROCAL
939 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(5);
940 ResetContainers(); Publish(); return; // issue retrieving callibration
941 }
6c0d2d67 942 // loads Q vectors
943 MakeQVectors();
944 if(fPsi2<-0.1) {
945 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(4);
946 ResetContainers(); Publish(); return;
24373b38 947 }
6c0d2d67 948 //}
24373b38 949 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("Events"))->Fill(3);
6c0d2d67 950 //=>great, lets do our stuff!
24373b38 951 //=>load candidates
952 if(!fSkipSelection) {
953 if(fReadESD) {
954 ReadFromESD(tESD);
955 } else {
956 if(fSpecie<10) ReadFromAODv0(tAOD);
957 else ChargeParticles(tAOD);
958 }
6c0d2d67 959 if(fUseFP) AddCandidates();
24373b38 960 //=>flow
961 //=>done
962 }
963 tTime.Stop();
964 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("RealTime"))->Fill( TMath::Log( tTime.RealTime() ) );
965 Publish();
966}
967//=======================================================================
968void AliAnalysisTaskFlowStrange::Publish() {
969 PostData(1,fList);
970 if(fUseFP) {
971 PostData(2,fTPCevent);
972 PostData(3,fVZEevent);
973 }
974}
975//=======================================================================
976void AliAnalysisTaskFlowStrange::ReadFromESD(AliESDEvent *tESD) {
977 AliStack *stack=NULL;
978 if(fReadMC) {
979 AliMCEvent *mcevent=NULL;
980 mcevent = MCEvent();
981 if(mcevent) stack = mcevent->Stack();
982 }
983
984 Int_t num = tESD->GetNumberOfTracks();
985 AliESDtrack *myTrack;
986 Int_t plist[3000], nlist[3000], np=0, nn=0;
987 Double_t pd0[3000], nd0[3000];
988 for (Int_t i=0; i!=num; ++i) {
989 myTrack = (AliESDtrack*) tESD->GetTrack(i);
990 if(!myTrack) continue;
991 LoadTrack(myTrack);
6c0d2d67 992 FillTrackSpy("ESD_TrkAll");
24373b38 993 if(!AcceptDaughter()) continue;
6c0d2d67 994 FillTrackSpy("ESD_TrkSel");
995 ((TH2D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("PtIPXY" ))->Fill( myTrack->Pt(), fDaughterImpactParameterXY );
24373b38 996 if( myTrack->Charge()>0 ) {
997 pd0[np] = fDaughterImpactParameterXY;
998 plist[np++] = i;
999 } else {
1000 nd0[nn] = fDaughterImpactParameterXY;
1001 nlist[nn++] = i;
1002 }
1003 }
6c0d2d67 1004 ((TH1D*)((TList*)fList->FindObject("ESD_TrkSel"))->FindObject("NPAIR" ))->Fill( np,nn );
24373b38 1005 const AliESDVertex *vtx = tESD->GetPrimaryVertex();
1006 AliESDtrack *pT, *nT;
1007 for(int p=0; p!=np; ++p) {
1008 pT = (AliESDtrack*) tESD->GetTrack( plist[p] );
1009 for(int n=0; n!=nn; ++n) {
1010 nT = (AliESDtrack*) tESD->GetTrack( nlist[n] );
1011 fDecayProductIPXY = pd0[p]*nd0[n];
1012 AliExternalTrackParam pETP(*pT), nETP(*nT);
1013 Double_t xa, xb;
1014 pETP.GetDCA(&nETP,tESD->GetMagneticField(),xa,xb);
1015 fDecayDCAdaughters = pETP.PropagateToDCA(&nETP,tESD->GetMagneticField());
1016 AliESDv0 vertex( nETP,nlist[n], pETP,plist[p] );
1017 fDecayCosinePointingAngleXY = CosThetaPointXY( &vertex, vtx );
1018 fDecayRadXY = DecayLengthXY( &vertex, vtx );
1019 fDecayPt = vertex.Pt();
1020 fDecayPhi = vertex.Phi();
1021 fDecayEta = vertex.Eta();
1022 Double_t pmx, pmy, pmz, nmx, nmy, nmz;
1023 vertex.GetNPxPyPz(nmx,nmy,nmz);
1024 vertex.GetPPxPyPz(pmx,pmy,pmz);
1025 TVector3 mom1(pmx,pmy,pmz), mom2(nmx,nmy,nmz), mom(vertex.Px(),vertex.Py(),vertex.Pz());
1026 Double_t qlpos = mom1.Dot(mom)/mom.Mag();
1027 Double_t qlneg = mom2.Dot(mom)/mom.Mag();
1028 fDecayQt = mom1.Perp(mom);
1029 fDecayAlpha = (qlpos-qlneg)/(qlpos+qlneg);
1030 Double_t mpi = 0.13957018;
1031 if(fSpecie==0) {
1032 Double_t eppi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1033 Double_t enpi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1034 fDecayMass = TMath::Sqrt( mpi*mpi + mpi*mpi + 2*(eppi*enpi - pmx*nmx - pmy*nmy - pmz*nmz ) );
1035 fDecayRapidity = vertex.RapK0Short();
1036 } else {
1037 Double_t mpr = 0.938272013;
1038 Double_t epi, epr;
1039 if(fDecayAlpha>0) {
1040 epr = TMath::Sqrt( mpr*mpr + pmx*pmx + pmy*pmy + pmz*pmz );
1041 epi = TMath::Sqrt( mpi*mpi + nmx*nmx + nmy*nmy + nmz*nmz );
1042 } else {
1043 epi = TMath::Sqrt( mpi*mpi + pmx*pmx + pmy*pmy + pmz*pmz );
1044 epr = TMath::Sqrt( mpr*mpr + nmx*nmx + nmy*nmy + nmz*nmz );
1045 }
1046 fDecayMass = TMath::Sqrt( mpi*mpi + mpr*mpr + 2*(epi*epr - pmx*nmx - pmy*nmy - pmz*nmz ) );
1047 fDecayRapidity = vertex.RapLambda();
1048 }
1049 Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + vertex.Px()*vertex.Px() + vertex.Py()*vertex.Py() + vertex.Pz()*vertex.Pz() );
1050 Double_t gamma = energy/fDecayMass;
1051 fDecayDecayLength = DecayLength( &vertex, vtx )/gamma;
1052 Double_t dPHI = fDecayPhi;
1053 Double_t dDPHI = dPHI - fPsi2;
1054 if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1055 if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1056 if(fQAlevel>1) {
6c0d2d67 1057 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1058 else FillCandidateSpy("V0SAllIP");
24373b38 1059 }
6c0d2d67 1060 FillCandidateSpy("V0SAll");
1061 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1062 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
24373b38 1063 if(!AcceptCandidate()) continue;
1064 if(fDecayMass<fMinMass) continue;
1065 if(fDecayMass>fMaxMass) continue;
1066 // PID missing
1067 if(fQAlevel>1) {
6c0d2d67 1068 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1069 else FillCandidateSpy("V0SSelIP");
24373b38 1070 }
6c0d2d67 1071 FillCandidateSpy("V0SSel");
1072 ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1073 ((TH2D*)((TList*)fList->FindObject("V0SSel"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
24373b38 1074
1075 fDecayIDneg = nT->GetID();
1076 fDecayIDpos = pT->GetID();
6c0d2d67 1077 if(fUseFP) MakeTrack();
1078 LoadTrack(pT); FillTrackSpy("V0SDau");
1079 LoadTrack(nT); FillTrackSpy("V0SDau");
24373b38 1080
1081 //===== BEGIN OF MCMATCH
1082 if(stack) {
1083 bool matched = false;
1084 Int_t labelpos = pT->GetLabel();
1085 Int_t labelneg = nT->GetLabel();
1086 Double_t rOri=-1;
1087 if( labelpos>0 && labelneg>0 ) {
1088 TParticle *mcpos = stack->Particle( labelpos );
1089 TParticle *mcneg = stack->Particle( labelneg );
1090 Int_t pdgRecPos = mcpos->GetPdgCode();
1091 Int_t pdgRecNeg = mcneg->GetPdgCode();
1092 if( pdgRecPos==211&&pdgRecNeg==-211 ) if(mcpos->GetMother(0)>0) {
1093 if( mcpos->GetMother(0)==mcneg->GetMother(0) ) {
1094 TParticle *mcmot = stack->Particle( mcpos->GetMother(0) );
1095 rOri = TMath::Sqrt( mcmot->Vx()*mcmot->Vx() + mcmot->Vy()*mcmot->Vy() );
1096 if( TMath::Abs(mcmot->GetPdgCode())==310) {
1097 if(mcmot->GetNDaughters()==2) matched=true;
1098 }
1099 }
1100 }
1101 }
1102 if(matched) {
6c0d2d67 1103 FillCandidateSpy("Mth");
1104 ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("D0PD0N"))->Fill( pd0[p],nd0[n] );
1105 ((TH2D*)((TList*)fList->FindObject("Mth"))->FindObject("XPOSXNEG"))->Fill( xa, xb );
1106 ((TH1D*)((TList*)fList->FindObject("Mth"))->FindObject("MCOrigin"))->Fill( rOri );
1107 LoadTrack(pT); FillTrackSpy("MthDau");
1108 LoadTrack(nT); FillTrackSpy("MthDau");
24373b38 1109 }
1110 }
1111 //===== END OF MCMATCH
1112 }
1113 }
1114}
1115//=======================================================================
1116void AliAnalysisTaskFlowStrange::ReadStack(TClonesArray* mcArray) {
1117 if(!mcArray) return;
6c0d2d67 1118 AliAODMCParticle *myMCTrack;//, *iMCDau, *jMCDau;
24373b38 1119 for(int i=0; i!=mcArray->GetEntriesFast(); ++i) {
1120 myMCTrack = dynamic_cast<AliAODMCParticle*>(mcArray->At( i ));
1121 if(!myMCTrack) continue;
6c0d2d67 1122 /*
24373b38 1123 int tPDG=310;
1124 if(fSpecie>0) tPDG = 3122;
1125 if( TMath::Abs(myMCTrack->PdgCode())==tPDG )
1126 if( myMCTrack->GetNDaughters() == 2 ) {
1127 Int_t iDau = myMCTrack->GetDaughter(0);
1128 Int_t jDau = myMCTrack->GetDaughter(1);
1129 AliAODMCParticle *posDau=NULL;
1130 AliAODMCParticle *negDau=NULL;
1131 if(iDau>0&&jDau>0) {
1132 iMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( iDau ));
1133 jMCDau = dynamic_cast<AliAODMCParticle*>(mcArray->At( jDau ));
1134 if(iMCDau) {
1135 if(iMCDau->Charge()>0) posDau=iMCDau;
1136 else negDau=iMCDau;
1137 }
1138 if(jMCDau) {
1139 if(jMCDau->Charge()>0) posDau=jMCDau;
1140 else negDau=jMCDau;
1141 }
1142 } //got two daughters
1143 if(posDau&&negDau) {
1144 Double_t dx = myMCTrack->Xv() - posDau->Xv();
1145 Double_t dy = myMCTrack->Yv() - posDau->Yv();
1146 Double_t dz = myMCTrack->Zv() - posDau->Zv();
1147 fDecayRadXY = TMath::Sqrt( dx*dx + dy*dy );
1148 TVector3 momPos(posDau->Px(),posDau->Py(),posDau->Pz());
1149 TVector3 momNeg(negDau->Px(),negDau->Py(),negDau->Pz());
1150 TVector3 momTot(myMCTrack->Px(),myMCTrack->Py(),myMCTrack->Pz());
1151 Double_t qlpos = momPos.Dot(momTot)/momTot.Mag();
1152 Double_t qlneg = momNeg.Dot(momTot)/momTot.Mag();
1153 fDecayQt = momPos.Perp(momTot);
1154 fDecayAlpha = 1.-2./(1.+qlpos/qlneg);
1155 fDecayMass = myMCTrack->GetCalcMass();
1156 Double_t energy = myMCTrack->E();
1157 Double_t gamma = energy/fDecayMass;
1158 fDecayDecayLength = TMath::Sqrt(dx*dx+dy*dy+dz*dz)/gamma;
1159 fDecayPt = myMCTrack->Pt();
1160 fDecayPhi = myMCTrack->Phi();
1161 fDecayEta = myMCTrack->Eta();
1162 fDecayRapidity = myMCTrack->Y();
1163 fDecayDCAdaughters = 0;
1164 fDecayCosinePointingAngleXY = 1;
1165 fDecayProductIPXY = -1;
1166 if(AcceptCandidate()) FillCandidateSpy("GenTru");
1167 }
1168 } // k0/lda with two daughters
6c0d2d67 1169 */
24373b38 1170 //==== BEGIN TRACK CUTS
1171 if(myMCTrack->Eta()<-0.8) continue;
1172 if(myMCTrack->Eta()>+0.8) continue;
1173 //==== END TRACK CUTS
1174 switch( TMath::Abs(myMCTrack->PdgCode()) ) {
6c0d2d67 1175 case (211): //pi
1176 FillMCParticleSpy( "MCTPion", myMCTrack );
1177 if( myMCTrack->IsPrimary() )
1178 FillMCParticleSpy( "MCTPionGenAcc", myMCTrack );
1179 break;
1180 case (321): //kaon
1181 FillMCParticleSpy( "MCTKaon", myMCTrack );
1182 if( myMCTrack->IsPrimary() )
1183 FillMCParticleSpy( "MCTKaonGenAcc", myMCTrack );
1184 break;
24373b38 1185 case (310): //k0s
1186 FillMCParticleSpy( "MCTK0s", myMCTrack );
1187 if( myMCTrack->IsPrimary() )
1188 FillMCParticleSpy( "MCTK0sGenAcc", myMCTrack );
1189 break;
6c0d2d67 1190 case (2212): //proton
1191 FillMCParticleSpy( "MCTProton", myMCTrack );
1192 if( myMCTrack->IsPrimary() )
1193 FillMCParticleSpy( "MCTProtonGenAcc", myMCTrack );
1194 break;
24373b38 1195 case (3122): //lda
1196 FillMCParticleSpy( "MCTLda", myMCTrack );
1197 if( myMCTrack->IsPrimary() )
1198 FillMCParticleSpy( "MCTLdaGenAcc", myMCTrack );
1199 break;
1200 case (333): //phi
1201 if( myMCTrack->IsPrimary() )
1202 FillMCParticleSpy( "MCTPhiGenAcc", myMCTrack );
1203 break;
1204 case (3312): //xi
1205 if( myMCTrack->IsPrimary() )
1206 FillMCParticleSpy( "MCTXiGenAcc", myMCTrack );
1207 break;
6c0d2d67 1208 case (3334): //omega
1209 if( myMCTrack->IsPrimary() )
1210 FillMCParticleSpy( "MCTOmegaGenAcc", myMCTrack );
1211 break;
24373b38 1212 }
24373b38 1213 }
1214}
1215//=======================================================================
1216Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliESDv0 *me, const AliVVertex *vtx) {
1217 TVector3 mom( me->Px(), me->Py(), 0 );
1218 TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1219 Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1220 return ctp;
1221}
1222//=======================================================================
1223Double_t AliAnalysisTaskFlowStrange::CosThetaPointXY(AliAODv0 *me, const AliVVertex *vtx) {
1224 TVector3 mom( me->Px(), me->Py(), 0 );
1225 TVector3 fli( me->Xv()-vtx->GetX(), me->Yv()-vtx->GetY(), 0 );
1226 Double_t ctp = mom.Dot(fli) / mom.Mag() / fli.Mag();
1227 return ctp;
1228}
1229//=======================================================================
1230Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliESDv0 *me, const AliVVertex *vtx) {
1231 Double_t dx = me->Xv()-vtx->GetX();
1232 Double_t dy = me->Yv()-vtx->GetY();
1233 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1234 return dxy;
1235}
1236//=======================================================================
1237Double_t AliAnalysisTaskFlowStrange::DecayLengthXY(AliAODv0 *me, const AliVVertex *vtx) {
1238 Double_t dx = me->Xv()-vtx->GetX();
1239 Double_t dy = me->Yv()-vtx->GetY();
1240 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy );
1241 return dxy;
1242}
1243//=======================================================================
1244Double_t AliAnalysisTaskFlowStrange::DecayLength(AliESDv0 *me, const AliVVertex *vtx) {
1245 Double_t dx = me->Xv()-vtx->GetX();
1246 Double_t dy = me->Yv()-vtx->GetY();
1247 Double_t dz = me->Zv()-vtx->GetZ();
1248 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1249 return dxy;
1250}
1251//=======================================================================
1252Double_t AliAnalysisTaskFlowStrange::DecayLength(AliAODv0 *me, const AliVVertex *vtx) {
1253 Double_t dx = me->Xv()-vtx->GetX();
1254 Double_t dy = me->Yv()-vtx->GetY();
1255 Double_t dz = me->Zv()-vtx->GetZ();
1256 Double_t dxy = TMath::Sqrt( dx*dx + dy*dy + dz*dz );
1257 return dxy;
1258}
1259//=======================================================================
1260void AliAnalysisTaskFlowStrange::ReadFromAODv0(AliAODEvent *tAOD) {
1261 TClonesArray* mcArray=NULL;
1262 if(fReadMC) {
1263 mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1264 ReadStack(mcArray);
1265 }
1266
1267 Int_t nV0s = tAOD->GetNumberOfV0s();
1268 AliAODv0 *myV0;
1269 Int_t v0all=0, v0imw=0;
1270 for (Int_t i=0; i!=nV0s; ++i) {
1271 myV0 = (AliAODv0*) tAOD->GetV0(i);
1272 if(!myV0) continue;
1273 if(!fOnline) if(myV0->GetOnFlyStatus() ) continue;
1274 if(fOnline) if(!myV0->GetOnFlyStatus() ) continue;
1275 AliAODTrack *iT, *jT;
1276 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1277 Double_t pos[3],cov[6];
1278 vtx->GetXYZ(pos);
1279 vtx->GetCovarianceMatrix(cov);
1280 const AliESDVertex vESD(pos,cov,100.,100);
1281 // TESTING CHARGE
1282 int iPos, iNeg;
1283 iT=(AliAODTrack*) myV0->GetDaughter(0);
1284 if(iT->Charge()>0) {
1285 iPos = 0; iNeg = 1;
1286 } else {
1287 iPos = 1; iNeg = 0;
1288 }
1289 // END OF TEST
1290
1291 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1292 AliESDtrack ieT( iT );
1293 ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1294 ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1295 ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1296 ieT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1297 LoadTrack(&ieT,iT->Chi2perNDF());
1298 Float_t ip[2];
1299 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1300 fDaughterImpactParameterXY = ip[0];
1301 fDaughterImpactParameterZ = ip[1];
512ced40 1302 fDecayIPpos = fDaughterImpactParameterXY; //ieT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
24373b38 1303 if(!AcceptDaughter()) continue;
1304
1305 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1306 AliESDtrack jeT( jT );
1307 jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1308 jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1309 jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1310 jeT.PropagateToDCA(&vESD, tAOD->GetMagneticField(), 100);
1311 LoadTrack(&jeT,jT->Chi2perNDF());
1312 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1313 fDaughterImpactParameterXY = ip[0];
1314 fDaughterImpactParameterZ = ip[1];
512ced40 1315 fDecayIPneg = fDaughterImpactParameterXY; //jeT.GetD(pos[0], pos[1], tAOD->GetMagneticField());
24373b38 1316 if(!AcceptDaughter()) continue;
512ced40 1317
24373b38 1318 if( fExcludeTPCEdges ) {
1319 if( IsAtTPCEdge(iT->Phi(),iT->Pt(),+1,tAOD->GetMagneticField()) ) continue;
1320 if( IsAtTPCEdge(jT->Phi(),jT->Pt(),-1,tAOD->GetMagneticField()) ) continue;
1321 }
512ced40 1322 ieT.GetDCA(&jeT,tAOD->GetMagneticField(),fDecayXpos,fDecayXneg);
24373b38 1323 /*
1324 // cutting out population close to TPC edges :: strange excess saw in 2010
1325 if( fExcludeTPCEdges ) {
1326 Double_t phimod = myV0->Phi();
1327 int sectors[6] = {5,6,9,10,11,12};
1328 for(int ii=0; ii!=6; ++ii)
1329 if( (phimod<(sectors[ii]+1)*TMath::Pi()/9.0) && (phimod>sectors[ii]*TMath::Pi()/9.0) )
1330 return 0;
1331 }
1332 */
1333 if(fSpecie==0)
1334 fDecayRapidity = myV0->RapK0Short();
1335 else
1336 fDecayRapidity = myV0->RapLambda();
1337 fDecayDCAdaughters = myV0->DcaV0Daughters();
1338 fDecayCosinePointingAngleXY = CosThetaPointXY( myV0, vtx );
1339 fDecayRadXY = DecayLengthXY( myV0, vtx );
1340 fDecayPt = myV0->Pt();
1341 fDecayPhi = myV0->Phi();
1342 fDecayEta = myV0->Eta();
512ced40 1343 fDecayProductIPXY = fDecayIPpos*fDecayIPneg;
24373b38 1344 fDecayQt = myV0->PtArmV0();
1345 fDecayAlpha = myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1346 if(myV0->ChargeProng(iPos)<0) fDecayAlpha = -fDecayAlpha; // protects for a change in convention
1347 fDecayPt = myV0->Pt();
1348 fDecayEta = myV0->Eta();
1349 if( fSpecie==0 ) {
1350 fDecayMass = myV0->MassK0Short();
1351 } else {
1352 if(fDecayAlpha>0) fDecayMass = myV0->MassLambda();
1353 else fDecayMass = myV0->MassAntiLambda();
1354 }
1355 v0all++;
1356 if(fDecayMass<fMinMass) continue;
1357 if(fDecayMass>fMaxMass) continue;
1358 v0imw++;
1359 Double_t energy = TMath::Sqrt( fDecayMass*fDecayMass + myV0->Px()*myV0->Px() + myV0->Py()*myV0->Py() + myV0->Pz()*myV0->Pz() );
1360 Double_t gamma = energy/fDecayMass;
1361 fDecayDecayLength = DecayLength( myV0, vtx )/gamma;
1362 Double_t dPHI = fDecayPhi;
1363 Double_t dDPHI = dPHI - fPsi2;
1364 if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
1365 if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
1366 if(fQAlevel>1) {
6c0d2d67 1367 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SAllOP");
1368 else FillCandidateSpy("V0SAllIP");
24373b38 1369 }
6c0d2d67 1370 FillCandidateSpy("V0SAll");
1371 FillDecayVn("V0SAllVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
24373b38 1372 if(!AcceptCandidate()) continue;
1373 //PID for lambda::proton
1374 if( fSpecie>0 )
1375 if(fDecayPt<1.2) {
1376 if(fDecayAlpha>0) {
1377 if( !PassesPIDCuts(&ieT) ) continue; //positive track
1378 } else {
1379 if( !PassesPIDCuts(&jeT) ) continue; //negative track
1380 }
1381 }
1382 if(fQAlevel>1) {
6c0d2d67 1383 if( (dDPHI>TMath::PiOver4()) && (dDPHI<3*TMath::PiOver4()) ) FillCandidateSpy("V0SSelOP");
1384 else FillCandidateSpy("V0SSelIP");
24373b38 1385 }
6c0d2d67 1386 FillCandidateSpy("V0SSel");
1387 FillDecayVn("V0SSelVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
512ced40
RAB		 1388 // ============================
1389 // Posting for FlowAnalysis
1390 if(!fPostMatched) {
1391 fDecayIDneg = iT->GetID();
1392 fDecayIDpos = jT->GetID();
6c0d2d67 1393 if(fUseFP) MakeTrack();
512ced40
RAB		 1394 }
1395 // ============================
1396
24373b38 1397 LoadTrack(&ieT,iT->Chi2perNDF());
1398 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1399 fDaughterImpactParameterXY = ip[0];
1400 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1401 FillTrackSpy("V0SDau");
24373b38 1402 LoadTrack(&jeT,jT->Chi2perNDF());
1403 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1404 fDaughterImpactParameterXY = ip[0];
1405 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1406 FillTrackSpy("V0SDau");
512ced40
RAB		 1407 //===== BEGIN OF MCMATCH
1408 if(fReadMC) ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event
24373b38 1409 if(mcArray) {
512ced40 1410 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
24373b38 1411 bool matched = false;
1412 bool feeddown = false;
1413 Int_t labelpos = iT->GetLabel();
1414 Int_t labelneg = jT->GetLabel();
512ced40
RAB		 1415 AliAODMCParticle *mcpos = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelpos) );
1416 AliAODMCParticle *mcneg = (AliAODMCParticle*) mcArray->At( TMath::Abs(labelneg) );
1417 if( mcpos && mcneg ) {
1418 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Daughters in stack
24373b38 1419 Int_t pdgRecPos = mcpos->GetPdgCode();
1420 Int_t pdgRecNeg = mcneg->GetPdgCode();
1421 int pospdg=211, negpdg=211;
1422 int mompdg=310, fdwpdg=333;
1423 if(fSpecie>0) {
1424 mompdg=3122;
1425 fdwpdg=3312;
1426 if(fDecayAlpha>0) {
1427 pospdg=2212; negpdg=211;
1428 } else {
1429 negpdg=2212; pospdg=211;
1430 }
1431 }
1432 if( TMath::Abs(pdgRecPos)==pospdg&&TMath::Abs(pdgRecNeg)==negpdg )
512ced40 1433 if(mcpos->GetMother()>-1)
24373b38 1434 if( mcpos->GetMother()==mcneg->GetMother() ) {
1435 AliAODMCParticle *mcmot = (AliAODMCParticle*) mcArray->At( mcpos->GetMother() );
512ced40
RAB		 1436 fDecayMatchOrigin = TMath::Sqrt( mcmot->Xv()*mcmot->Xv() + mcmot->Yv()*mcmot->Yv() );
1437 fDecayMatchPt = mcmot->Pt();
1438 fDecayMatchEta = mcmot->Eta();
1439 fDecayMatchPhi = mcmot->Phi();
24373b38 1440 if( TMath::Abs(mcmot->GetPdgCode())==mompdg) {
1441 if(mcmot->GetNDaughters()==2) {
512ced40 1442 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 4 ); // Correspond to decay
24373b38 1443 matched=true;
1444 Double_t dx = mcmot->Xv() - mcpos->Xv();
1445 Double_t dy = mcmot->Yv() - mcpos->Yv();
512ced40 1446 fDecayMatchRadXY = TMath::Sqrt( dx*dx + dy*dy );
24373b38 1447 }
512ced40
RAB		 1448 if(mcmot->GetMother()>-1) {
1449 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 5 ); // Decay has mother
24373b38 1450 AliAODMCParticle *mcfdw = (AliAODMCParticle*) mcArray->At( mcmot->GetMother() );
1451 if( TMath::Abs(mcfdw->GetPdgCode())==fdwpdg)
1452 feeddown=true;
1453 } // k0/lda have mother
1454 } // mother matches k0/lda
1455 } // both have same mother
512ced40 1456 }
24373b38 1457 if(matched) {
6c0d2d67 1458 FillCandidateSpy("Mth",true);
1459 FillDecayVn("V0SMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
512ced40
RAB		 1460 if(fPostMatched>0) {
1461 fDecayIDneg = iT->GetID();
1462 fDecayIDpos = jT->GetID();
6c0d2d67 1463 if(fUseFP) MakeTrack();
512ced40
RAB		 1464 }
1465 if(labelpos<0&&labelneg<0) {
6c0d2d67 1466 FillCandidateSpy("MthNegNeg",true);
1467 FillDecayVn("V0SMthNegNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
1468 } else if(labelpos>0&&labelneg>0) {
1469 FillDecayVn("V0SMthPosPosVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
512ced40 1470 } else if(labelpos*labelneg<0) {
6c0d2d67 1471 FillCandidateSpy("MthPosNeg",true);
1472 FillDecayVn("V0SMthPosNegVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
512ced40 1473 }
6c0d2d67 1474 AliAODVertex *secvtx = myV0->GetSecondaryVtx();
1475 Double_t possec[3],covsec[6];
1476 secvtx->GetXYZ(possec);
1477 secvtx->GetCovarianceMatrix(covsec);
1478 const AliESDVertex vSecVtx(possec,covsec,100.,100);
1479 AliESDtrack trackAtSecI( iT );
1480 trackAtSecI.SetTPCClusterMap( iT->GetTPCClusterMap() );
1481 trackAtSecI.SetTPCSharedMap( iT->GetTPCSharedMap() );
1482 trackAtSecI.SetTPCPointsF( iT->GetTPCNclsF() );
1483 trackAtSecI.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1484 fDaughterAtSecPhi = trackAtSecI.Phi();
1485 fDaughterAtSecEta = trackAtSecI.Eta();
1486 fDaughterAtSecPt = trackAtSecI.Pt();
24373b38 1487 LoadTrack(&ieT,iT->Chi2perNDF());
512ced40
RAB		 1488 fDaughterMatchPhi=mcpos->Phi();
1489 fDaughterMatchEta=mcpos->Eta();
1490 fDaughterMatchPt=mcpos->Pt();
24373b38 1491 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1492 fDaughterImpactParameterXY = ip[0];
1493 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1494 FillTrackSpy("MthDau",true);
1495 if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1496 else FillTrackSpy("MthPosDau",true);
1497 AliESDtrack trackAtSecJ( jT );
1498 trackAtSecJ.SetTPCClusterMap( jT->GetTPCClusterMap() );
1499 trackAtSecJ.SetTPCSharedMap( jT->GetTPCSharedMap() );
1500 trackAtSecJ.SetTPCPointsF( jT->GetTPCNclsF() );
1501 trackAtSecJ.PropagateToDCA(&vSecVtx, tAOD->GetMagneticField(), 100);
1502 fDaughterAtSecPhi = trackAtSecJ.Phi();
1503 fDaughterAtSecEta = trackAtSecJ.Eta();
1504 fDaughterAtSecPt = trackAtSecJ.Pt();
24373b38 1505 LoadTrack(&jeT,jT->Chi2perNDF());
512ced40
RAB		 1506 fDaughterMatchPhi=mcneg->Phi();
1507 fDaughterMatchEta=mcneg->Eta();
1508 fDaughterMatchPt=mcneg->Pt();
24373b38 1509 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1510 fDaughterImpactParameterXY = ip[0];
1511 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1512 FillTrackSpy("MthDau",true);
1513 if(labelpos<0||labelneg<0) FillTrackSpy("MthNegDau",true);
1514 else FillTrackSpy("MthPosDau",true);
512ced40 1515 } else {
6c0d2d67 1516 FillCandidateSpy("UnMth",false);
1517 FillDecayVn("V0SUnMthVn",fDecayMass,fDecayPt,fDecayPhi,fDecayEta,fDecayIDpos,fDecayIDneg);
512ced40
RAB		 1518 if(fPostMatched<0) {
1519 fDecayIDneg = iT->GetID();
1520 fDecayIDpos = jT->GetID();
6c0d2d67 1521 if(fUseFP) MakeTrack();
512ced40
RAB		 1522 }
1523 LoadTrack(&ieT,iT->Chi2perNDF());
512ced40
RAB		 1524 ieT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1525 fDaughterImpactParameterXY = ip[0];
1526 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1527 FillTrackSpy("UnMthDau",false);
512ced40 1528 LoadTrack(&jeT,jT->Chi2perNDF());
512ced40
RAB		 1529 jeT.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1530 fDaughterImpactParameterXY = ip[0];
1531 fDaughterImpactParameterZ = ip[1];
6c0d2d67 1532 FillTrackSpy("UnMthDau",false);
24373b38 1533 }
1534 if(feeddown) {
6c0d2d67 1535 FillCandidateSpy("MthFeedDown",true);
24373b38 1536 }
1537 }
1538 //===== END OF MCMATCH
1539 }
6c0d2d67 1540 ((TH2D*)((TList*)fList->FindObject("V0SAll"))->FindObject("V0SADC"))->Fill( v0all,v0imw );
1541 QCStoreDecayVn("V0SAllVn");
1542 QCStoreDecayVn("V0SSelVn");
1543 if(fReadMC) {
1544 QCStoreDecayVn("V0SMthVn");
1545 QCStoreDecayVn("V0SMthNegNegVn");
1546 QCStoreDecayVn("V0SMthPosPosVn");
1547 QCStoreDecayVn("V0SMthPosNegVn");
1548 }
24373b38 1549 return;
1550}
1551//=======================================================================
1552Bool_t AliAnalysisTaskFlowStrange::PassesPIDCuts(AliESDtrack *myTrack, AliPID::EParticleType pid) {
1553 Bool_t pass=kTRUE;
1554 if(fPIDResponse) {
1555 fDaughterNSigmaPID = fPIDResponse->NumberOfSigmasTPC(myTrack,pid);
1556 if( TMath::Abs(fDaughterNSigmaPID) > fDaughterMaxNSigmaPID )
1557 pass = kFALSE;
1558 }
1559 return pass;
1560}
1561//=======================================================================
1562void AliAnalysisTaskFlowStrange::ChargeParticles(AliAODEvent *tAOD) {
6c0d2d67 1563 //benchmark purposes
24373b38 1564 if(!tAOD) return;
6c0d2d67 1565 TClonesArray* mcArray=NULL;
1566 if(fReadMC) {
1567 mcArray = dynamic_cast<TClonesArray*>(tAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1568 ReadStack(mcArray);
1569 }
24373b38 1570 for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
1571 AliAODTrack *t = tAOD->GetTrack( i );
1572 if(!t) continue;
6c0d2d67 1573 if( !t->TestFilterBit(1) ) continue;
24373b38 1574 fDecayMass=0.0; // using mass as nsigmas control plot
1575 if(fPIDResponse) { // PID
1576 switch(fSpecie) { // TPC PID only
1577 case(kPION):
1578 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kPion);
1579 break;
1580 case(kKAON):
1581 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kKaon);
1582 break;
1583 case(kPROTON):
1584 fDecayMass = fPIDResponse->NumberOfSigmasTPC(t,AliPID::kProton);
1585 break;
1586 }
1587 }
1588 Bool_t pass = kTRUE;
1589 if( TMath::Abs(fDecayMass) > 3.0 ) pass=kFALSE;
6c0d2d67 1590 if( t->Eta()<-0.5 || t->Eta()>+0.5 ) pass=kFALSE;
24373b38 1591 if( t->Pt()<0.2 || t->Pt()>20.0 ) pass=kFALSE;
24373b38 1592 AliESDtrack et( t );
1593 et.SetTPCClusterMap( t->GetTPCClusterMap() );
1594 et.SetTPCSharedMap( t->GetTPCSharedMap() );
1595 et.SetTPCPointsF( t->GetTPCNclsF() );
1596 Float_t ip[2];
1597 LoadTrack(&et,t->Chi2perNDF());
1598 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
1599 Double_t pos[3];
1600 vtx->GetXYZ(pos);
1601 et.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
1602 fDaughterImpactParameterXY = ip[0];
1603 fDaughterImpactParameterZ = ip[1];
24373b38 1604
6c0d2d67 1605 FillTrackSpy("TrkAll");
1606 FillTrackVn("TrkAllVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1607 if(!pass) continue;
1608 FillTrackSpy("TrkSel");
1609 FillTrackVn("TrkSelVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1610 if(fReadMC) {
1611 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 1 ); // Selected event
1612 if(mcArray) {
1613 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 2 ); // Stack found
1614 bool matched = false;
1615 Int_t label = t->GetLabel();
1616 AliAODMCParticle *mcpar = (AliAODMCParticle*) mcArray->At( TMath::Abs(label) );
1617 if( mcpar ) {
1618 ((TH1D*)((TList*)fList->FindObject("STATMC"))->FindObject("Events"))->Fill( 3 ); // Particle in stack
1619 Int_t pdgmcpar = TMath::Abs(mcpar->GetPdgCode());
1620 switch(fSpecie) {
1621 case(kPION):
1622 if(pdgmcpar==211) matched = true;
1623 break;
1624 case(kKAON):
1625 if(pdgmcpar==211) matched = true;
1626 break;
1627 case(kPROTON):
1628 if(pdgmcpar==2212) matched = true;
1629 break;
1630 }
1631 if(!mcpar->IsPrimary()) matched = false;
1632 }
1633 if(matched) {
1634 FillTrackSpy("Mth");
1635 FillTrackVn("MthVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1636 if(label<0) {
1637 FillTrackSpy("MthNeg");
1638 FillTrackVn("MthNegVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1639 } else {
1640 FillTrackSpy("MthPos");
1641 FillTrackVn("MthPosVn",t->Pt(),t->Phi(),t->Eta(),t->GetID());
1642 }
1643 }
1644 }
1645 }
1646 if(fUseFP) {
1647 fDecayPt=t->Pt();
1648 fDecayPhi=t->Phi();
1649 fDecayEta=t->Eta();
1650 fDecayID=t->GetID();
1651 MakeTrack();
1652 }
1653 }
1654 QCStoreTrackVn("TrkAllVn");
1655 QCStoreTrackVn("TrkSelVn");
1656 if(fReadMC) {
1657 QCStoreTrackVn("MthVn");
1658 QCStoreTrackVn("MthNegVn");
1659 QCStoreTrackVn("MthPosVn");
24373b38 1660 }
1661 return;
1662}
1663//=======================================================================
1664void AliAnalysisTaskFlowStrange::Terminate(Option_t *) {
1665 //terminate
6c0d2d67 1666 if(fSkipTerminate) return;
1667 Double_t MeanQaQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->GetBinContent( 1 );
1668 Double_t MeanQaQa = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->GetBinContent( 1 );
1669 Double_t MeanQcQc = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->GetBinContent( 1 );
1670 Double_t MeanQaQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->GetBinContent( 1 );
1671 Double_t MeanQcQt = ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->GetBinContent( 1 );
1672 if(!TMath::AreEqualAbs(MeanQaQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQcQt,0,1e-10)&&!TMath::AreEqualAbs(MeanQaQc,0,1e-10)) {
1673 Double_t OneOverChiSquaredVZEA = MeanQaQa*MeanQcQt/MeanQaQc/MeanQaQt-1;
1674 Double_t OneOverChiSquaredVZEC = MeanQcQc*MeanQaQt/MeanQaQc/MeanQcQt-1;
1675 if(!TMath::AreEqualAbs(OneOverChiSquaredVZEA,0,1e-10)&&!TMath::AreEqualAbs(OneOverChiSquaredVZEC,0,1e-10)) {
1676 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->SetBinContent( 1, 1/OneOverChiSquaredVZEA );
1677 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->SetBinContent( 1, 1/OneOverChiSquaredVZEC );
1678 }
1679 }
1680 if(fSkipSelection) return;
1681 if(fSpecie<10) {
1682 ComputeDecayVn("V0SAllVn");
1683 ComputeDecayVn("V0SSelVn");
1684 if(fReadMC) {
1685 ComputeDecayVn("V0SMthVn");
1686 ComputeDecayVn("V0SMthPosPosVn");
1687 ComputeDecayVn("V0SMthNegNegVn");
1688 ComputeDecayVn("V0SMthPosNegVn");
1689 ComputeDecayVn("V0SUnMthVn");
1690 }
1691 } else {
1692 ComputeTrackVn("TrkAllVn");
1693 ComputeTrackVn("TrkSelVn");
1694 if(fReadMC) {
1695 ComputeTrackVn("MthVn");
1696 ComputeTrackVn("MthPosVn");
1697 ComputeTrackVn("MthNegVn");
1698 }
1699 }
24373b38 1700}
1701//=======================================================================
1702void AliAnalysisTaskFlowStrange::MakeTrack() {
1703 // create track for flow tasks
1704 if(fCandidates->GetLast()+5>fCandidates->GetSize()) {
1705 fCandidates->Expand( fCandidates->GetSize()+20 );
1706 }
1707 Bool_t overwrite = kTRUE;
1708 AliFlowCandidateTrack *oTrack = (static_cast<AliFlowCandidateTrack*> (fCandidates->At( fCandidates->GetLast()+1 )));
1709 if( !oTrack ) { // creates new
1710 oTrack = new AliFlowCandidateTrack();
1711 overwrite = kFALSE;
1712 } else { // overwrites
1713 oTrack->ClearMe();
1714 }
1715 oTrack->SetMass(fDecayMass);
1716 oTrack->SetPt(fDecayPt);
1717 oTrack->SetPhi(fDecayPhi);
1718 oTrack->SetEta(fDecayEta);
1719 if(fSpecie<10) {
1720 oTrack->AddDaughter(fDecayIDpos);
1721 oTrack->AddDaughter(fDecayIDneg);
1722 } else {
1723 oTrack->SetID( fDecayID );
1724 }
1725 oTrack->SetForPOISelection(kTRUE);
1726 oTrack->SetForRPSelection(kFALSE);
1727 if(overwrite) {
1728 fCandidates->SetLast( fCandidates->GetLast()+1 );
1729 } else {
1730 fCandidates->AddLast(oTrack);
1731 }
1732 return;
1733}
1734//=======================================================================
1735void AliAnalysisTaskFlowStrange::AddCandidates() {
1736 // adds candidates to flow events (untaging if necessary)
1737 if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1738 if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1739 if(fDebug) printf("I received %d candidates\n",fCandidates->GetEntriesFast());
1740 Int_t untagged=0;
1741 Int_t poi=0;
1742 for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
1743 AliFlowCandidateTrack *cand = static_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
1744 if(!cand) continue;
1745 cand->SetForPOISelection(kTRUE);
1746 cand->SetForRPSelection(kFALSE);
1747 poi++;
6c0d2d67 1748 //if(fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
24373b38 1749 if(fSpecie<10) { // DECAYS
1750 // untagging ===>
1751 if(fDaughterUnTag) {
1752 for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
1753 if(fDebug) printf(" >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
1754 for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1755 AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1756 if(!iRP) continue;
1757 if(!iRP->InRPSelection()) continue;
1758 if(cand->GetIDDaughter(iDau) == iRP->GetID()) {
1759 if(fDebug) printf(" was in RP set");
1760 ++untagged;
1761 iRP->SetForRPSelection(kFALSE);
1762 fTPCevent->SetNumberOfRPs( fTPCevent->GetNumberOfRPs() -1 );
1763 }
1764 }
1765 if(fDebug) printf("\n");
1766 }
1767 }
1768 // <=== untagging
1769 fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1770 } else { // CHARGED
1771 // adding only new tracks and tagging accordingly ===>
1772 Bool_t found=kFALSE;
1773 for(int iRPs=0; iRPs!=fTPCevent->NumberOfTracks(); ++iRPs ) {
1774 AliFlowTrack *iRP = static_cast<AliFlowTrack*>(fTPCevent->GetTrack( iRPs ));
1775 if(!iRP) continue;
1776 if(!iRP->InRPSelection()) continue;
1777 if(cand->GetID() == iRP->GetID()) {
1778 if(fDebug) printf(" >charged track (%d) was also found in RP set (adding poi tag)\n",cand->GetID());
6c0d2d67 1779 iRP->SetMass( cand->Mass() );
24373b38 1780 iRP->SetForPOISelection(kTRUE);
1781 found = kTRUE;
1782 }
1783 }
1784 if(!found) // not found adding track
1785 fTPCevent->InsertTrack( ((AliFlowTrack*) cand) );
1786 }
1787 fVZEevent->InsertTrack( ((AliFlowTrack*) cand) );
1788 } //END OF LOOP
1789 fTPCevent->SetNumberOfPOIs( poi );
1790 fVZEevent->SetNumberOfPOIs( poi );
1791 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("POI"))->Fill( poi );
1792 ((TH1D*)((TList*)fList->FindObject("Event"))->FindObject("UNTAG"))->Fill( untagged );
1793 if(fDebug) printf("FlowEventTPC %d tracks | %d RFP | %d POI\n",fTPCevent->NumberOfTracks(),fTPCevent->GetNumberOfRPs(),fTPCevent->GetNumberOfPOIs());
1794 if(fDebug) printf("FlowEventVZE %d tracks | %d RFP | %d POI\n",fVZEevent->NumberOfTracks(),fVZEevent->GetNumberOfRPs(),fVZEevent->GetNumberOfPOIs());
1795}
1796//=======================================================================
1797void AliAnalysisTaskFlowStrange::PushBackFlowTrack(AliFlowEvent *flowevent, Double_t pt, Double_t phi, Double_t eta, Double_t w, Int_t id) {
1798 AliFlowTrack rfp;
1799 rfp.SetPt(pt);
1800 rfp.SetPhi(phi);
1801 rfp.SetEta(eta);
1802 rfp.SetWeight(w);
1803 rfp.SetForRPSelection(kTRUE);
1804 rfp.SetForPOISelection(kFALSE);
6c0d2d67 1805 rfp.SetMass(-999);
24373b38 1806 rfp.SetID(id);
1807 flowevent->InsertTrack( &rfp );
1808}
1809//=======================================================================
1810Bool_t AliAnalysisTaskFlowStrange::IsAtTPCEdge(Double_t phi,Double_t pt,Int_t charge,Double_t b) {
1811 // Origin: Alex Dobrin
1812 // Implemented by Carlos Perez
1813 TF1 cutLo("cutLo", "-0.01/x+pi/18.0-0.015", 0, 100);
1814 TF1 cutHi("cutHi", "0.55/x/x+pi/18.0+0.03", 0, 100);
1815 Double_t phimod = phi;
1816 if(b<0) phimod = TMath::TwoPi()-phimod; //for negatve polarity field
1817 if(charge<0) phimod = TMath::TwoPi()-phimod; //for negatve charge
1818 phimod += TMath::Pi()/18.0;
1819 phimod = fmod(phimod, TMath::Pi()/9.0);
1820 if( phimod<cutHi.Eval(pt) && phimod>cutLo.Eval(pt) )
1821 return kTRUE;
1822
1823 return kFALSE;
1824}
1825//=======================================================================
1826void AliAnalysisTaskFlowStrange::MakeQVectors() {
1827 //computes event plane and updates fPsi2
1828 //if there is a problem fPsi->-1
1829 fPsi2=-1;
24373b38 1830 //=>loading event
6c0d2d67 1831 MakeQVZE(InputEvent());
1832 MakeQTPC(InputEvent());
1833 if(fUseFP&&fReadMC) {
c41a93af 1834 fVZEevent->SetMCReactionPlaneAngle( fMCEP );
1835 fTPCevent->SetMCReactionPlaneAngle( fMCEP );
1836 }
6c0d2d67 1837 if(fDebug) {
1838 printf("**::MakeQVectors()");
1839 printf(" fQVZEACos %.16f | fQVZEASin %.16f || fQVZEA %.3f | fQVZEC %.3f \n",fQVZEACos, fQVZEASin, fQVZEA, fQVZEC);
1840 printf(" nQTPA_nTracks %d | fQTPC_nTracks %d || fQTPCA %.3f | fQTPCC %.3f \n",fQTPCA_nTracks, fQTPCC_nTracks, fQTPCA, fQTPCC);
1841 printf(" fQTPCACos %.16f | fQTPCASin %.16f || fQTPC2hCos %.16f | fQTPC2hSin %.16f \n",fQTPCACos, fQTPCASin, fQTPC2hCos, fQTPC2hSin);
1842 }
24373b38 1843 //=>computing psi
1844 //VZERO
6c0d2d67 1845 Double_t qvzecos,qvzesin,psivzea,psivzec,psivze,qvze;
1846 psivzea = ( TMath::Pi()+TMath::ATan2(-fQVZEASin,-fQVZEACos) )/fHarmonic;
1847 psivzec = ( TMath::Pi()+TMath::ATan2(-fQVZECSin,-fQVZECCos) )/fHarmonic;
1848 qvzecos = fQVZEACos + fQVZECCos;
1849 qvzesin = fQVZEASin + fQVZECSin;
1850 qvze = fQVZEA + fQVZEC;
1851 psivze = ( TMath::Pi()+TMath::ATan2(-qvzesin,-qvzecos) )/fHarmonic;
24373b38 1852 //TPC
6c0d2d67 1853 Double_t qtpccos,qtpcsin,psitpca,psitpcc,psitpc,qtpc;
1854 psitpca = ( TMath::Pi()+TMath::ATan2(-fQTPCASin,-fQTPCACos) )/fHarmonic;
1855 psitpcc = ( TMath::Pi()+TMath::ATan2(-fQTPCCSin,-fQTPCCCos) )/fHarmonic;
1856 qtpccos = fQTPCACos + fQTPCCCos;
1857 qtpcsin = fQTPCASin + fQTPCCSin;
1858 qtpc = fQTPCA + fQTPCC;
1859 psitpc = ( TMath::Pi()+TMath::ATan2(-qtpcsin,-qtpccos) )/fHarmonic;
24373b38 1860 //=>does the event clear?
1861 switch(fWhichPsi) {
1862 case(1): //VZERO
6c0d2d67 1863 if(fQVZEA<2||fQVZEC<2) return;
24373b38 1864 fPsi2 = psivze;
1865 break;
1866 case(2): //TPC
6c0d2d67 1867 if(fQTPCA<2||fQTPCC<2) return;
24373b38 1868 fPsi2 = psitpc;
1869 break;
1870 }
6c0d2d67 1871 //computing physical Qm vectors
1872 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
1873 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
1874 Double_t vzec_qmnor = TMath::Sqrt( vzec_qmcos*vzec_qmcos + vzec_qmsin*vzec_qmsin );
1875 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
1876 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
1877 Double_t vzea_qmnor = TMath::Sqrt( vzea_qmcos*vzea_qmcos + vzea_qmsin*vzea_qmsin );
1878 Double_t vze_qmcos = qvzecos/qvze;
1879 Double_t vze_qmsin = qvzesin/qvze;
1880 Double_t vze_qmnor = TMath::Sqrt( vze_qmcos*vze_qmcos + vze_qmsin*vze_qmsin );
1881 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
1882 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
1883 Double_t tpcc_qmnor = TMath::Sqrt( tpcc_qmcos*tpcc_qmcos + tpcc_qmsin*tpcc_qmsin );
1884 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
1885 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
1886 Double_t tpca_qmnor = TMath::Sqrt( tpca_qmcos*tpca_qmcos + tpca_qmsin*tpca_qmsin );
1887 Double_t tpc_qmcos = qtpccos/qtpc;
1888 Double_t tpc_qmsin = qtpcsin/qtpc;
1889 Double_t tpc_qmnor = TMath::Sqrt( tpc_qmcos*tpc_qmcos + tpc_qmsin*tpc_qmsin );
24373b38 1890 //=>great! recording
1891 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSI"))->Fill( psivze );
1892 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIA"))->Fill( psivzea );
6c0d2d67 1893 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEPSIC"))->Fill( psivzec );
1894 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPVZE"))->Fill( qvze );
1895 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZE"))->Fill( vze_qmnor );
24373b38 1896 //------
1897 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSI"))->Fill( psitpc );
1898 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIA"))->Fill( psitpca );
6c0d2d67 1899 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCPSIC"))->Fill( psitpcc );
1900 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("RFPTPC"))->Fill( qtpc );
1901 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPC"))->Fill( tpc_qmnor );
24373b38 1902 //------
6c0d2d67 1903 if(fReadMC) {
1904 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPC"))->Fill( psitpc-fMCEP );
1905 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCA"))->Fill( psitpca-fMCEP );
1906 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFTPCC"))->Fill( psitpcc-fMCEP );
1907 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZE"))->Fill( psivze-fMCEP );
1908 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEA"))->Fill( psivzea-fMCEP );
1909 ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("PSIMCDIFFVZEC"))->Fill( psivzec-fMCEP );
1910 }
24373b38 1911 //------
6c0d2d67 1912 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 1., tpcc_qmsin, tpcc_qmnor );
1913 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 2., tpcc_qmcos, tpcc_qmnor );
1914 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 3., tpca_qmsin, tpca_qmnor );
1915 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 4., tpca_qmcos, tpca_qmnor );
1916 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 5., tpc_qmsin, tpc_qmnor );
1917 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCQm"))->Fill( 6., tpc_qmcos, tpc_qmnor );
1918 //------
1919 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 1., vzec_qmsin, vzec_qmnor );
1920 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 2., vzec_qmcos, vzec_qmnor );
1921 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 3., vzea_qmsin, vzea_qmnor );
1922 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 4., vzea_qmcos, vzea_qmnor );
1923 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 5., vze_qmsin, vze_qmnor );
1924 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEQm"))->Fill( 6., vze_qmcos, vze_qmnor );
1925 //------
1926 Double_t vzeqaqc = vzec_qmcos*vzea_qmcos + vzec_qmsin*vzea_qmsin;
1927 Double_t vzeqatpcq = vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin;
1928 Double_t vzeqctpcq = vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin;
1929 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEAQmVZEC"))->Fill( 1., vzeqaqc, vze_qmnor );
1930 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZEASQUARED"))->Fill( 1., vzea_qmnor*vzea_qmnor, vze_qmnor );
1931 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmVZECSQUARED"))->Fill( 1., vzec_qmnor*vzec_qmnor, vze_qmnor );
1932 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEA"))->Fill( 1., vzeqatpcq, vze_qmnor );
1933 ((TProfile*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("QmTPCQmVZEC"))->Fill( 1., vzeqctpcq, vze_qmnor );
24373b38 1934 return;
1935}
1936//=======================================================================
6c0d2d67 1937void AliAnalysisTaskFlowStrange::MakeQVZE(AliVEvent *tevent) {
24373b38 1938 //=>cleaning
6c0d2d67 1939 if(fUseFP) fVZEevent->ClearFast(); // flowpackage
24373b38 1940 //=>external weights
1941 Double_t extW[64];
1942 for(int i=0;i!=64;++i) extW[i]=1;
1943 if((!fVZEsave)&&(fVZEResponse)) {
1944 Double_t minC = fCentPerMin, maxC = fCentPerMax;
1945 if(fVZEmb) {
1946 minC = 0;
1947 maxC = 80;
1948 }
1949 Int_t ybinmin = fVZEResponse->GetYaxis()->FindBin(minC+1e-6);
1950 Int_t ybinmax = fVZEResponse->GetYaxis()->FindBin(maxC-1e-6);
6c0d2d67 1951 if(fSkipCentralitySelection) {
1952 ybinmin=-1;
1953 ybinmax=-1;
1954 }
24373b38 1955 for(int i=0;i!=64;++i) extW[i] = fVZEResponse->Integral(i+1,i+1,ybinmin,ybinmax)/(maxC-minC);
1956 //ring-wise normalization
1957 Double_t ring[8];
1958 for(int j=0; j!=8; ++j) {
1959 ring[j]=0;
1960 for(int i=0;i!=8;++i) ring[j] += extW[j*8+i]/8;
1961 }
1962 //disk-wise normalization
1963 Double_t disk[2];
1964 int xbinmin, xbinmax;
1965 xbinmin = 1+8*fVZECa;
1966 xbinmax = 8+8*fVZECb;
1967 disk[0] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
1968 xbinmin = 33+8*fVZEAa;
1969 xbinmax = 40+8*fVZEAb;
1970 disk[1] = fVZEResponse->Integral(xbinmin,xbinmax,ybinmin,ybinmax)/(maxC-minC)/(xbinmax-xbinmin+1);
1971 //for(int i=0;i!=64;++i) printf("CELL %d -> W = %f ||",i,extW[i]);
1972
1973 if(fVZEByDisk) {
1974 for(int i=0;i!=64;++i) extW[i] = disk[i/32]/extW[i];
1975 } else {
1976 for(int i=0;i!=64;++i) extW[i] = ring[i/8]/extW[i];
1977 }
1978 //for(int i=0;i!=64;++i) printf(" W = %f \n",extW[i]);
1979 }
1980 //=>computing
6c0d2d67 1981 fQVZEACos=fQVZEASin=fQVZEA=fQVZECCos=fQVZECSin=fQVZEC=0;
24373b38 1982 Int_t rfp=0;
1983 Double_t eta, phi, w;
1984 //v0c -> qa
1985 for(int id=fVZECa*8;id!=8+fVZECb*8;++id) {
1986 eta = -3.45+0.5*(id/8);
1987 phi = TMath::PiOver4()*(0.5+id%8);
1988 w = tevent->GetVZEROEqMultiplicity(id)*extW[id];
6c0d2d67 1989 fQVZECCos += w*TMath::Cos(fHarmonic*phi);
1990 fQVZECSin += w*TMath::Sin(fHarmonic*phi);
1991 fQVZEC += w;
24373b38 1992 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
1993 rfp++;
6c0d2d67 1994 if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
24373b38 1995 }
1996 //v0a -> qb
1997 for(int id=32+fVZEAa*8;id!=40+fVZEAb*8;++id) {
1998 eta = +4.8-0.6*((id/8)-4);
1999 phi = TMath::PiOver4()*(0.5+id%8);
2000 w = tevent->GetVZEROEqMultiplicity(id)*extW[id];
6c0d2d67 2001 fQVZEACos += w*TMath::Cos(fHarmonic*phi);
2002 fQVZEASin += w*TMath::Sin(fHarmonic*phi);
2003 fQVZEA += w;
24373b38 2004 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("VZEAllPhiEta"))->Fill( phi, eta, w );
2005 rfp++;
6c0d2d67 2006 if(fUseFP) PushBackFlowTrack(fVZEevent,0,phi,eta,w,0); // flowpackage
2007 }
2008 if(fUseFP) { // flowpackage
2009 fVZEevent->SetNumberOfRPs(rfp);
2010 if(fDebug>0) printf("Inserted tracks in FlowEventVZE %d ==> %.1f\n",fVZEevent->NumberOfTracks(),fQVZEC+fQVZEA);
24373b38 2011 }
24373b38 2012}
2013//=======================================================================
2014void AliAnalysisTaskFlowStrange::AddTPCRFPSpy(TList *me) {
2015 TH1D *tH1D;
2016 tH1D = new TH1D("PT", "PT", 50,0,5); me->Add(tH1D);
2017 tH1D = new TH1D("PHI", "PHI", 90,0,TMath::TwoPi()); me->Add(tH1D);
2018 tH1D = new TH1D("ETA", "ETA", 40,-1,+1); me->Add(tH1D);
2019 tH1D = new TH1D("TPCS", "TPC Signal", 100,0,500); me->Add(tH1D);
2020 tH1D = new TH1D("IPXY", "IPXY", 100,-2,+2); me->Add(tH1D);
2021 tH1D = new TH1D("IPZ", "IPZ", 100,-2,+2); me->Add(tH1D);
2022 // TPC
2023 tH1D = new TH1D("TPCNCLS", "NCLS", 170,-0.5,+169.5); me->Add(tH1D);
2024 tH1D = new TH1D("TPCSHCL", "NSCLS / NCLS", 100,0,1); me->Add(tH1D);
2025 tH1D = new TH1D("TPCFICL", "NCLS1I / NCLS",100,0,1); me->Add(tH1D);
2026 tH1D = new TH1D("TPCXRNF", "XROW / NFCLS", 100,0,1.5); me->Add(tH1D);
2027 tH1D = new TH1D("TPCRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
2028 // ITS
2029 tH1D = new TH1D("ITSNCLS", "NCLS", 7,-0.5,+6.5); me->Add(tH1D);
2030 tH1D = new TH1D("ITSRCHI", "CHI2 / NCLS", 50,0,5); me->Add(tH1D);
2031
2032}
2033//=======================================================================
2034Bool_t AliAnalysisTaskFlowStrange::PassesRFPTPCCuts(AliESDtrack *track, Double_t aodchi2cls, Float_t aodipxy, Float_t aodipz) {
2035 if(track->GetKinkIndex(0)>0) return kFALSE;
2036 if( (track->GetStatus()&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2037 Double_t pt = track->Pt();
2038 Double_t phi = track->Phi();
2039 Double_t eta = track->Eta();
2040 Double_t tpcs = track->GetTPCsignal();
2041 Float_t ipxy, ipz;
2042 track->GetImpactParameters(ipxy,ipz);
2043 Int_t cls = track->GetTPCclusters(0);
2044 Double_t xrows, findcls, chi2;
2045 findcls = track->GetTPCNclsF();
2046 xrows = track->GetTPCCrossedRows();
2047 chi2 = track->GetTPCchi2();
2048 Double_t rchi2 = chi2/cls;
2049 if(!fReadESD) {
2050 rchi2 = aodchi2cls;
2051 ipxy = aodipxy;
2052 ipz = aodipz;
2053 }
2054 Double_t xrnfcls = xrows/findcls;
2055 Double_t scls, cls1i, itschi2;
2056 Int_t itscls;
2057 cls1i = track->GetTPCNclsIter1();
2058 scls = track->GetTPCnclsS();
2059 itscls = track->GetITSclusters(0);
2060 itschi2 = track->GetITSchi2();
2061 Double_t shcl = scls/cls;
2062 Double_t ficl = cls1i/cls;
2063 Double_t itsrchi2 = itscls/itschi2;
2064 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PT"))->Fill( pt );
2065 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("PHI"))->Fill( phi );
2066 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ETA"))->Fill( eta );
2067 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCS"))->Fill( tpcs );
2068 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPXY"))->Fill( ipxy );
2069 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("IPZ"))->Fill( ipz );
2070 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCNCLS"))->Fill( cls );
2071 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCSHCL"))->Fill( shcl );
2072 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCFICL"))->Fill( ficl );
2073 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2074 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2075 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSNCLS"))->Fill( itscls );
2076 ((TH1D*)((TList*)fList->FindObject("TPCRFPall"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2077 if(pt<fRFPminPt) return kFALSE; //0.2
2078 if(pt>fRFPmaxPt) return kFALSE; //5.0
6c0d2d67 2079 if(eta<fRFPCminEta||(eta>fRFPCmaxEta&&eta<fRFPAminEta)||eta>fRFPAmaxEta) return kFALSE; // -0.8 0.0 0.0 +0.8
24373b38 2080 if(tpcs<fRFPTPCsignal) return kFALSE; //10.0
2081 if( TMath::Sqrt(ipxy*ipxy/fRFPmaxIPxy/fRFPmaxIPxy+ipz*ipz/fRFPmaxIPz/fRFPmaxIPz)>1 ) return kFALSE; //2.4 3.2
2082 if(cls<fRFPTPCncls) return kFALSE; //70
2083 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PT"))->Fill( pt );
2084 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("PHI"))->Fill( phi );
2085 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ETA"))->Fill( eta );
2086 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCS"))->Fill( tpcs );
2087 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPXY"))->Fill( ipxy );
2088 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("IPZ"))->Fill( ipz );
2089 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCNCLS"))->Fill( cls );
2090 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCSHCL"))->Fill( shcl );
2091 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCFICL"))->Fill( ficl );
2092 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCXRNF"))->Fill( xrnfcls );
2093 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("TPCRCHI"))->Fill( rchi2 );
2094 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSNCLS"))->Fill( itscls );
2095 ((TH1D*)((TList*)fList->FindObject("TPCRFPsel"))->FindObject("ITSRCHI"))->Fill( itsrchi2 );
2096 return kTRUE;
2097}
2098//=======================================================================
6c0d2d67 2099void AliAnalysisTaskFlowStrange::MakeQTPC(AliVEvent *tevent) {
24373b38 2100 AliESDEvent *tESD = (AliESDEvent*) (tevent);
2101 AliAODEvent *tAOD = (AliAODEvent*) (tevent);
2102 if(fReadESD) {
2103 if(!tESD) return;
6c0d2d67 2104 MakeQTPC(tESD);
24373b38 2105 } else {
2106 if(!tAOD) return;
6c0d2d67 2107 MakeQTPC(tAOD);
24373b38 2108 }
2109}
2110//=======================================================================
6c0d2d67 2111void AliAnalysisTaskFlowStrange::MakeQTPC(AliAODEvent *tAOD) {
24373b38 2112 //=>cleaning
6c0d2d67 2113 if(fUseFP) fTPCevent->ClearFast(); // flowpackage
2114 fQTPCACos=fQTPCASin=fQTPCA=fQTPC2hCos=0;
2115 fQTPCCCos=fQTPCCSin=fQTPCC=fQTPC2hSin=0;
2116 fQTPCA_nTracks = 0;
2117 fQTPCC_nTracks = 0;
24373b38 2118 Int_t rfp=0;
2119 Double_t eta, phi, w;
2120 //=>aod stuff
2121 AliAODVertex *vtx = tAOD->GetPrimaryVertex();
2122 Double_t pos[3],cov[6];
2123 vtx->GetXYZ(pos);
2124 vtx->GetCovarianceMatrix(cov);
2125 const AliESDVertex vESD(pos,cov,100.,100);
2126 AliAODTrack *track;
2127 //=>looping
2128 Int_t rawN = tAOD->GetNumberOfTracks();
2129 for(Int_t id=0; id!=rawN; ++id) {
2130 track = tAOD->GetTrack(id);
2131 //=>cuts
2132 if(!track->TestFilterBit(fRFPFilterBit)) continue;
2133 if( fExcludeTPCEdges )
2134 if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tAOD->GetMagneticField() ) ) continue;
2135 AliESDtrack etrack( track );
2136 etrack.SetTPCClusterMap( track->GetTPCClusterMap() );
2137 etrack.SetTPCSharedMap( track->GetTPCSharedMap() );
2138 etrack.SetTPCPointsF( track->GetTPCNclsF() );
2139 Float_t ip[2];
2140 etrack.GetDZ(pos[0], pos[1], pos[2], tAOD->GetMagneticField(), ip);
2141 if(!PassesRFPTPCCuts(&etrack,track->Chi2perNDF(),ip[0],ip[1])) continue;
2142 //=>collecting info
2143 phi = track->Phi();
2144 eta = track->Eta();
2145 w = 1;
6c0d2d67 2146 if(eta<0) {
2147 fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2148 fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2149 fQTPCC += w;
2150 fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
24373b38 2151 } else {
6c0d2d67 2152 fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2153 fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2154 fQTPCA += w;
2155 fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
24373b38 2156 }
6c0d2d67 2157 fQTPC2hCos += w*TMath::Cos(2.0*fHarmonic*phi);
2158 fQTPC2hSin += w*TMath::Sin(2.0*fHarmonic*phi);
24373b38 2159 rfp++;
6c0d2d67 2160 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2161 if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flow package
2162 }
2163 if(fUseFP) {
2164 fTPCevent->SetNumberOfRPs(rfp);
2165 if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
24373b38 2166 }
24373b38 2167}
2168//=======================================================================
6c0d2d67 2169void AliAnalysisTaskFlowStrange::MakeQTPC(AliESDEvent *tESD) {
24373b38 2170 //=>cleaning
6c0d2d67 2171 if(fUseFP) fTPCevent->ClearFast(); // flow package
2172 fQTPCACos=fQTPCASin=fQTPCA=0;
2173 fQTPCCCos=fQTPCCSin=fQTPCC=0;
2174 fQTPCA_nTracks = 0;
2175 fQTPCC_nTracks = 0;
24373b38 2176 Int_t rfp=0;
2177 Double_t eta, phi, w;
2178 //=>looping
2179 AliESDtrack *track;
2180 Int_t rawN = tESD->GetNumberOfTracks();
2181 for(Int_t id=0; id!=rawN; ++id) {
2182 track = tESD->GetTrack(id);
2183 //=>cuts
2184 if( fExcludeTPCEdges )
2185 if( IsAtTPCEdge( track->Phi(), track->Pt(), track->Charge(), tESD->GetMagneticField() ) ) continue;
2186 if(!PassesFilterBit(track)) continue;
2187 if(!PassesRFPTPCCuts(track)) continue;
2188 //=>collecting info
2189 phi = track->Phi();
2190 eta = track->Eta();
2191 w = 1;
24373b38 2192 if(eta<0) {
6c0d2d67 2193 fQTPCCCos += w*TMath::Cos(fHarmonic*phi);
2194 fQTPCCSin += w*TMath::Sin(fHarmonic*phi);
2195 fQTPCC += w;
2196 fQTPCC_fID[fQTPCC_nTracks++] = track->GetID();
24373b38 2197 } else {
6c0d2d67 2198 fQTPCACos += w*TMath::Cos(fHarmonic*phi);
2199 fQTPCASin += w*TMath::Sin(fHarmonic*phi);
2200 fQTPCA += w;
2201 fQTPCA_fID[fQTPCA_nTracks++] = track->GetID();
24373b38 2202 }
24373b38 2203 rfp++;
6c0d2d67 2204 ((TH2D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("TPCAllPhiEta"))->Fill( phi, eta, w );
2205 if(fUseFP) PushBackFlowTrack(fTPCevent,track->Pt(),phi,eta,w,track->GetID()); // flowpackage
2206 }
2207 if(fUseFP) {
2208 fTPCevent->SetNumberOfRPs(rfp);
2209 if(fDebug) printf("Inserted tracks in FlowEventTPC %d ==> %.1f\n",fTPCevent->NumberOfTracks(),fQTPCA+fQTPCC);
24373b38 2210 }
24373b38 2211}
2212//=======================================================================
2213void AliAnalysisTaskFlowStrange::AddMCParticleSpy(TList *me) {
2214 TH1D *tH1D;
2215 TH2D *tH2D;
512ced40 2216 TProfile *tPro;
24373b38 2217 tH1D = new TH1D("Pt", "Pt", 100,0.0,20); me->Add(tH1D);
2218 tH1D = new TH1D("Phi", "Phi", 100,0,TMath::TwoPi()); me->Add(tH1D);
2219 tH1D = new TH1D("Eta", "Eta", 100,-1,+1); me->Add(tH1D);
2220 tH1D = new TH1D("Rad2", "Rad2", 1000,0,+100); me->Add(tH1D);
2221 tH2D = new TH2D("Dphi", "phi-MCEP;pt;dphi",100,0,20, 72,0,TMath::Pi()); me->Add(tH2D);
512ced40 2222 tPro = new TProfile("Cos2dphi","Cos2dphi",100,0,20); me->Add(tPro);
24373b38 2223 return;
2224}
2225//=======================================================================
2226void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, AliAODMCParticle *p) {
2227 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2228 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2229 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2230 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Xv()*p->Xv() +
2231 p->Yv()*p->Yv() ) );
2232 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Dphi" ))->Fill( p->Pt(), GetMCDPHI(p->Phi()) );
512ced40 2233 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphi" ))->Fill( p->Pt(), TMath::Cos( 2*GetMCDPHI(p->Phi()) ), 1 );
24373b38 2234 return;
2235}
2236//=======================================================================
2237Double_t AliAnalysisTaskFlowStrange::GetMCDPHI(Double_t phi) {
2238 Double_t dDPHI = phi - fMCEP;
6c0d2d67 2239 //if( dDPHI < 0 ) dDPHI += TMath::TwoPi();
2240 //if( dDPHI > TMath::Pi() ) dDPHI = TMath::TwoPi()-dDPHI;
24373b38 2241 return dDPHI;
2242}
2243//=======================================================================
2244void AliAnalysisTaskFlowStrange::FillMCParticleSpy(TString listName, TParticle *p) {
2245 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Pt" ))->Fill( p->Pt() );
2246 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Eta" ))->Fill( p->Eta() );
2247 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Phi" ))->Fill( p->Phi() );
2248 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("Rad2" ))->Fill( TMath::Sqrt( p->Vx()*p->Vx() +
2249 p->Vy()*p->Vy() ) );
2250 return;
2251}
2252//=======================================================================
512ced40
RAB		 2253void AliAnalysisTaskFlowStrange::AddCandidatesSpy(TList *me,Bool_t res) {
2254 TH1D *tH1D;
24373b38 2255 TH2D *tH2D;
512ced40
RAB		 2256 TProfile *tPro;
2257 TProfile2D *tPro2;
24373b38 2258 tH2D = new TH2D("PhiEta", "PhiEta;Phi;Eta", 100,0,TMath::TwoPi(),100,-1,+1); me->Add(tH2D);
2259 tH2D = new TH2D("PtRAP", "PtRAP;Pt;Y", 100,0,20,100,-2.0,+2.0); me->Add(tH2D);
2260 tH2D = new TH2D("PtDCA", "PtDCA;Pt;DCA", 100,0,20,100,0,10); me->Add(tH2D);
2261 tH2D = new TH2D("PtCTP", "PtCTP;Pt;CTP", 100,0,20,100,-1,+1); me->Add(tH2D);
24373b38 2262 tH2D = new TH2D("PtD0D0", "PtD0D0;Pt;D0D0", 100,0,20,100,-5,+5); me->Add(tH2D);
2263 tH2D = new TH2D("PtRad2", "PtRad2;Pt;RadXY", 100,0,20,100,0,+50); me->Add(tH2D);
2264 tH2D = new TH2D("PtDL", "PtDL;Pt;DL", 100,0,20,100,0,+50); me->Add(tH2D);
2265 tH2D = new TH2D("PtMASS", "PtMASS;Pt;MASS", 100,0,20,fMassBins,fMinMass,fMaxMass); me->Add(tH2D);
2266 tH2D = new TH2D("APPOS", "APPOS;alphaPOS;QtPOS",100,-2,+2,100,0,0.3); me->Add(tH2D);
2267 tH2D = new TH2D("D0PD0N", "D0PD0N;D0P;D0N", 200,-10,+10,200,-10,+10); me->Add(tH2D);
2268 tH2D = new TH2D("XPOSXNEG","XPOSXNEG;XPOS;XNEG", 200,-50,+50,200,-50,+50); me->Add(tH2D);
512ced40
RAB		 2269
2270 if(fReadMC) {
2271 if(res) {
2272 tH1D = new TH1D("MCOrigin", "MCOrigin;Rad2",1000,0,50); me->Add(tH1D);
2273 tH2D = new TH2D("PHIRes","PHIRes;PHI;MC-DAT", 72, 0, TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2274 tH2D = new TH2D("ETARes","ETARes;ETA;MC-DAT", 16,-0.8, +0.8,100,-0.2,+0.2); me->Add(tH2D);
2275 tH2D = new TH2D("PTRes", "PTRes;Pt;MC-DAT", 100, 0, 20,100,-0.4,+0.4); me->Add(tH2D);
2276 tH2D = new TH2D("RXYRes","RXYRes;RXY;MC-DAT",100, 0, 50,100,-4.0,+4.0); me->Add(tH2D);
2277 }
2278 tH2D = new TH2D("PTDPHIMC","PtDPHIMC;Pt;PHI-MCEP",100,0,20,72,0,TMath::Pi()); me->Add(tH2D);
2279 tPro = new TProfile("Cos2dphiMC", "Cos2dphiMC",100,0,20); me->Add(tPro);
2280 tPro = new TProfile("Cos2dphiMCPK","Cos2dphiMC PK",100,0,20); me->Add(tPro);
2281 tPro2=new TProfile2D("C2DPHIMCMASS","C2DPHIMCMASS",100,0,20,fMassBins,fMinMass,fMaxMass); me->Add(tPro2);
2282 }
24373b38 2283 return;
2284}
2285//=======================================================================
512ced40 2286void AliAnalysisTaskFlowStrange::FillCandidateSpy(TString listName, Bool_t fillRes) {
24373b38 2287 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PhiEta"))->Fill( fDecayPhi, fDecayEta );
2288 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRAP" ))->Fill( fDecayPt, fDecayRapidity );
2289 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDCA" ))->Fill( fDecayPt, fDecayDCAdaughters );
2290 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtCTP" ))->Fill( fDecayPt, fDecayCosinePointingAngleXY );
2291 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtD0D0"))->Fill( fDecayPt, fDecayProductIPXY );
2292 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtRad2"))->Fill( fDecayPt, fDecayRadXY );
2293 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtDL" ))->Fill( fDecayPt, fDecayDecayLength );
2294 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PtMASS"))->Fill( fDecayPt, fDecayMass );
2295 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("APPOS" ))->Fill( fDecayAlpha, fDecayQt );
512ced40
RAB		 2296 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("D0PD0N"))->Fill( fDecayIPpos, fDecayIPneg );
2297 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("XPOSXNEG"))->Fill( fDecayXpos, fDecayXneg );
2298 if(fReadMC) {
2299 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTDPHIMC" ))->Fill( fDecayPt, GetMCDPHI( fDecayPhi ) );
2300 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphiMC" ))->Fill( fDecayPt, TMath::Cos( 2*GetMCDPHI(fDecayPhi) ), 1 );
2301 if( fDecayMass>fMinMassX && fDecayMass<fMaxMassX ) {
2302 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("Cos2dphiMCPK" ))->Fill( fDecayPt, TMath::Cos( 2*GetMCDPHI(fDecayPhi) ), 1 );
2303 }
2304 ((TProfile2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("C2DPHIMCMASS" ))->Fill(fDecayPt,fDecayMass,TMath::Cos(2*GetMCDPHI(fDecayPhi)), 1 );
2305 if(fillRes) {
2306 ((TH1D*)((TList*)fList->FindObject(listName.Data()))->FindObject("MCOrigin"))->Fill( fDecayMatchOrigin );
2307 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDecayPhi, fDecayMatchPhi-fDecayPhi );
2308 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDecayEta, fDecayMatchEta-fDecayEta );
2309 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDecayPt, fDecayMatchPt-fDecayPt );
2310 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("RXYRes"))->Fill( fDecayRadXY, fDecayMatchRadXY-fDecayRadXY );
2311 }
2312 }
24373b38 2313}
2314//=======================================================================
2315Bool_t AliAnalysisTaskFlowStrange::AcceptCandidate() {
2316 if(fDecayEta<fDecayMinEta) return kFALSE;
2317 if(fDecayEta>fDecayMaxEta) return kFALSE;
2318 if(fDecayPt<fDecayMinPt) return kFALSE;
2319 if(fDecayProductIPXY>fDecayMaxProductIPXY) return kFALSE;
2320 if(fDecayDCAdaughters>fDecayMaxDCAdaughters) return kFALSE;
2321 if(fDecayCosinePointingAngleXY<fDecayMinCosinePointingAngleXY) return kFALSE;
2322 if(fDecayRadXY<fDecayMinRadXY) return kFALSE;
24373b38 2323 if(TMath::Abs(fDecayRapidity)>fDecayMaxRapidity) return kFALSE;
2324 if(fSpecie==0) {
2325 if(fDecayAPCutPie) {
2326 if(fDecayQt/TMath::Abs(fDecayAlpha)<fDecayMinQt) return kFALSE;
2327 } else {
2328 if(fDecayQt<fDecayMinQt) return kFALSE;
2329 }
512ced40
RAB		 2330 if(fDecayDecayLength>fDecayMaxDecayLength*2.6842) return kFALSE;
2331 } else {
2332 if(fDecayDecayLength>fDecayMaxDecayLength*7.89) return kFALSE;
24373b38 2333 }
2334 if(fSpecie==1) if(fDecayAlpha>0) return kFALSE;
2335 if(fSpecie==2) if(fDecayAlpha<0) return kFALSE;
2336 return kTRUE;
2337}
2338//=======================================================================
6c0d2d67 2339void AliAnalysisTaskFlowStrange::AddTrackSpy(TList *me,Bool_t res) {
24373b38 2340 TH2D *tH2D;
2341 tH2D = new TH2D("PHIETA", "PHIETA;PHI;ETA", 100,0,TMath::TwoPi(),100,-2,2); me->Add(tH2D);
2342 tH2D = new TH2D("IPXYIPZ", "IPXYIPZ;IPXY;IPZ", 1000,-20,+20,1000,-20,+20); me->Add(tH2D);
2343 tH2D = new TH2D("PTTPCNCLS", "PTTPCNCLS;PT;NCLS", 100,0,20,170,0,170); me->Add(tH2D);
512ced40
RAB		 2344 tH2D = new TH2D("PTITSLAY", "PTITSLAY;PT;ITSLAYER", 100,0,20,6,-0.5,+5.5);me->Add(tH2D);
2345 tH2D = new TH2D("PTITSTPCrefit","PTITSTPCrefit;PT", 100,0,20,2,-0.5,+1.5);me->Add(tH2D);
2346 tH2D->GetYaxis()->SetBinLabel(1,"ITS refit");
2347 tH2D->GetYaxis()->SetBinLabel(2,"TPC refit");
2348
24373b38 2349 tH2D = new TH2D("POSTPCNCLCHI2","POSTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8); me->Add(tH2D);
2350 tH2D = new TH2D("POSTPCNFCLNXR","POSTPCNFCLNXR;NFCLS;NXR", 170,0,170,170,0,170); me->Add(tH2D);
2351 tH2D = new TH2D("POSTPCNCLNFCL","POSTPCNCLNFCL;NCLS;NFCLS", 170,0,170,170,0,170); me->Add(tH2D);
2352 tH2D = new TH2D("POSTPCNCLNSCL","POSTPCNCLNSCL;NCLS;NSCLS", 170,0,170,170,0,170); me->Add(tH2D);
512ced40 2353
24373b38 2354 tH2D = new TH2D("NEGTPCNCLCHI2","NEGTPCNCLCHI2;NCLS;CHI2/NCLS", 170,0,170,100,0,8); me->Add(tH2D);
2355 tH2D = new TH2D("NEGTPCNFCLNXR","NEGTPCNFCLNXR;NFCLS;NXR", 170,0,170,170,0,170); me->Add(tH2D);
2356 tH2D = new TH2D("NEGTPCNCLNFCL","NEGTPCNCLNFCL;NCLS;NFCLS", 170,0,170,170,0,170); me->Add(tH2D);
2357 tH2D = new TH2D("NEGTPCNCLNSCL","NEGTPCNCLNSCL;NCLS;NSCLS", 170,0,170,170,0,170); me->Add(tH2D);
512ced40 2358
6c0d2d67 2359 if(fReadMC) {
2360 TProfile *tPro;
2361 tPro = new TProfile("COSNDPHIMC","COSNDPHIMC",100,0,20); me->Add(tPro);
2362 }
2363
512ced40
RAB		 2364 if(res) {
2365 tH2D = new TH2D("PHIRes", "PHIRes;PHI;MC-DAT", 72, 0, TMath::TwoPi(),100,-0.12,+0.12); me->Add(tH2D);
2366 tH2D = new TH2D("ETARes", "ETARes;ETA;MC-DAT", 16,-0.8, +0.8,100,-0.2,+0.2); me->Add(tH2D);
2367 tH2D = new TH2D("PTRes", "PTRes;Pt;MC-DAT", 100, 0, 20,100,-0.4,+0.4); me->Add(tH2D);
2368 }
24373b38 2369}
2370//=======================================================================
512ced40 2371void AliAnalysisTaskFlowStrange::FillTrackSpy(TString listName,Bool_t res) {
24373b38 2372 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PHIETA" ))->Fill( fDaughterPhi, fDaughterEta );
2373 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "IPXYIPZ" ))->Fill( fDaughterImpactParameterXY, fDaughterImpactParameterZ );
2374 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTTPCNCLS" ))->Fill( fDaughterPt, fDaughterNClsTPC );
512ced40
RAB		 2375 if( TESTBIT(fDaughterITScm,0) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 0 );
2376 if( TESTBIT(fDaughterITScm,1) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 1 );
2377 if( TESTBIT(fDaughterITScm,2) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 2 );
2378 if( TESTBIT(fDaughterITScm,3) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 3 );
2379 if( TESTBIT(fDaughterITScm,4) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 4 );
2380 if( TESTBIT(fDaughterITScm,5) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSLAY" ))->Fill( fDaughterPt, 5 );
2381 if( (fDaughterStatus&AliESDtrack::kITSrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 0 );
2382 if( (fDaughterStatus&AliESDtrack::kTPCrefit) ) ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( "PTITSTPCrefit" ))->Fill( fDaughterPt, 1 );
2383
24373b38 2384 TString ch="NEG";
2385 if(fDaughterCharge>0) ch="POS";
2386 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLCHI2",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterChi2PerNClsTPC );
2387 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNFCLNXR",ch.Data()) ))->Fill( fDaughterNFClsTPC, fDaughterXRows );
2388 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNFCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNFClsTPC );
2389 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject( Form("%sTPCNCLNSCL",ch.Data()) ))->Fill( fDaughterNClsTPC, fDaughterNSClsTPC );
512ced40 2390
6c0d2d67 2391 if(fReadMC) {
2392 Double_t cosn = TMath::Cos( fHarmonic*GetMCDPHI(fDaughterPhi) );
2393 ((TProfile*)((TList*)fList->FindObject(listName.Data()))->FindObject("COSNDPHIMC" ))->Fill( fDaughterPt, cosn, 1 );
2394 }
512ced40 2395 if(res) {
6c0d2d67 2396 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PHIRes"))->Fill( fDaughterPhi, fDaughterMatchPhi-fDaughterAtSecPhi );
2397 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("ETARes"))->Fill( fDaughterEta, fDaughterMatchEta-fDaughterAtSecEta );
2398 ((TH2D*)((TList*)fList->FindObject(listName.Data()))->FindObject("PTRes"))->Fill( fDaughterPt, fDaughterMatchPt-fDaughterAtSecPt );
512ced40 2399 }
24373b38 2400}
2401//=======================================================================
2402void AliAnalysisTaskFlowStrange::LoadTrack(AliESDtrack *myTrack, Double_t aodChi2NDF) {
2403 fDaughterCharge = myTrack->Charge();
2404 fDaughterXRows = myTrack->GetTPCCrossedRows();
2405 fDaughterNFClsTPC = myTrack->GetTPCNclsF();
2406 fDaughterNSClsTPC = myTrack->GetTPCnclsS();
2407 fDaughterNClsTPC = myTrack->GetTPCclusters(0);
2408 if(fReadESD) {
2409 if(fDaughterNClsTPC>0) fDaughterChi2PerNClsTPC = myTrack->GetTPCchi2()/fDaughterNClsTPC;
2410 } else {
2411 fDaughterChi2PerNClsTPC = aodChi2NDF;
2412 }
2413 myTrack->GetImpactParameters(fDaughterImpactParameterXY,fDaughterImpactParameterZ);
2414 fDaughterStatus = myTrack->GetStatus();
512ced40 2415 fDaughterITScm = myTrack->GetITSClusterMap();
24373b38 2416 fDaughterPhi = myTrack->Phi();
2417 fDaughterEta = myTrack->Eta();
2418 fDaughterPt = myTrack->Pt();
2419 fDaughterKinkIndex = myTrack->GetKinkIndex(0);
2420}
2421//=======================================================================
2422Bool_t AliAnalysisTaskFlowStrange::AcceptDaughter() {
2423 if(fDaughterKinkIndex>0) return kFALSE;
2424 if( (fDaughterStatus&AliESDtrack::kTPCrefit)==0 ) return kFALSE;
2425 if(fDaughterNFClsTPC<1) return kFALSE;
2426 if(fDaughterPt<fDaughterMinPt) return kFALSE;
2427 if(fDaughterEta<fDaughterMinEta) return kFALSE;
2428 if(fDaughterEta>fDaughterMaxEta) return kFALSE;
2429 if(fDaughterNClsTPC<fDaughterMinNClsTPC) return kFALSE;
2430 if(fDaughterXRows<fDaughterMinXRows) return kFALSE;
2431 if(fDaughterChi2PerNClsTPC>fDaughterMaxChi2PerNClsTPC) return kFALSE;
2432 if(TMath::Abs(fDaughterImpactParameterXY)<fDaughterMinImpactParameterXY) return kFALSE;
2433 if(fDaughterXRows<fDaughterMinXRowsOverNClsFTPC*fDaughterNFClsTPC) return kFALSE;
512ced40
RAB		 2434 for(Int_t lay=0; lay!=6; ++lay)
2435 if(fDaughterITSConfig[lay]>-0.5) {
2436 if(fDaughterITSConfig[lay]) {
2437 if(!TESTBIT(fDaughterITScm,lay)) return kFALSE;
2438 } else {
2439 if(TESTBIT(fDaughterITScm,lay)) return kFALSE;
2440 }
2441 }
24373b38 2442 return kTRUE;
2443}
2444//=======================================================================
2445Double_t AliAnalysisTaskFlowStrange::GetWDist(const AliVVertex* v0, const AliVVertex* v1) {
2446 // calculate sqrt of weighted distance to other vertex
2447 if (!v0 || !v1) {
2448 printf("One of vertices is not valid\n");
2449 return 0;
2450 }
2451 static TMatrixDSym vVb(3);
2452 double dist = -1;
2453 double dx = v0->GetX()-v1->GetX();
2454 double dy = v0->GetY()-v1->GetY();
2455 double dz = v0->GetZ()-v1->GetZ();
2456 double cov0[6],cov1[6];
2457 v0->GetCovarianceMatrix(cov0);
2458 v1->GetCovarianceMatrix(cov1);
2459 vVb(0,0) = cov0[0]+cov1[0];
2460 vVb(1,1) = cov0[2]+cov1[2];
2461 vVb(2,2) = cov0[5]+cov1[5];
2462 vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2463 vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2464 vVb.InvertFast();
2465 if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
2466 dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2467 + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2468 return dist>0 ? TMath::Sqrt(dist) : -1;
2469}
2470//=======================================================================
2471Bool_t AliAnalysisTaskFlowStrange::plpMV(const AliVEvent *event) {
2472 // check for multi-vertexer pile-up
2473 const AliAODEvent *aod = (const AliAODEvent*)event;
2474 const AliESDEvent *esd = (const AliESDEvent*)event;
2475 //
2476 const int kMinPlpContrib = 5;
2477 const double kMaxPlpChi2 = 5.0;
2478 const double kMinWDist = 15;
2479 //
2480 if (!aod && !esd) {
2481 printf("Event is neither of AOD nor ESD\n");
2482 exit(1);
2483 }
2484 //
2485 const AliVVertex* vtPrm = 0;
2486 const AliVVertex* vtPlp = 0;
2487 int nPlp = 0;
2488 //
2489 if (aod) {
2490 if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
2491 vtPrm = aod->GetPrimaryVertex();
2492 if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
2493 }
2494 else {
2495 if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
2496 vtPrm = esd->GetPrimaryVertexTracks();
2497 if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
2498 }
2499 //int bcPrim = vtPrm->GetBC();
2500 //
2501 for (int ipl=0;ipl<nPlp;ipl++) {
2502 vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
2503 //
2504 if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2505 if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2506 // int bcPlp = vtPlp->GetBC();
2507 // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
2508 //
2509 double wDst = GetWDist(vtPrm,vtPlp);
2510 if (wDst<kMinWDist) continue;
2511 //
2512 return kTRUE; // pile-up: well separated vertices
2513 }
2514 //
2515 return kFALSE;
2516 //
2517}
2518//=======================================================================
2519void AliAnalysisTaskFlowStrange::MakeFilterBits() {
2520 //FilterBit 1
2521 fFB1 = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
2522 //FilterBit1024
2523 fFB1024 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
2524 fFB1024->SetMinNCrossedRowsTPC(120);
2525 fFB1024->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
2526 fFB1024->SetMaxChi2PerClusterITS(36);
2527 fFB1024->SetMaxFractionSharedTPCClusters(0.4);
2528 fFB1024->SetMaxChi2TPCConstrainedGlobal(36);
2529 fFB1024->SetEtaRange(-0.9,0.9);
2530 fFB1024->SetPtRange(0.15, 1e10);
2531}
2532//=======================================================================
2533Bool_t AliAnalysisTaskFlowStrange::PassesFilterBit(AliESDtrack *track) {
2534 Bool_t ret=kFALSE;
2535 switch(fRFPFilterBit) {
2536 case(1024):
2537 ret = fFB1024->AcceptTrack(track);
2538 break;
2539 default:
2540 ret = fFB1->AcceptTrack(track);
2541 }
2542 return ret;
2543}
2544//=======================================================================
2545void AliAnalysisTaskFlowStrange::LoadVZEROResponse() {
2546 if(fVZEResponse) {
2547 TString run = fVZEResponse->GetTitle();
2548 if( run.Atoi() == fRunNumber ) return;
2549 fVZEResponse = NULL;
2550 }
2551 //==>loading
2552 fVZEResponse = dynamic_cast<TH2D*> (fVZEload->FindObject( Form("%d",fRunNumber) ));
2553 if(fVZEResponse) {
2554 printf("New VZE calibration: run %d || %s -> Entries %.0f\n",fRunNumber, fVZEResponse->GetTitle(),fVZEResponse->GetEntries());
2555 } else {
6c0d2d67 2556 printf("VZE calibration: requested but not found!!!\n");
24373b38 2557 }
2558}
2559//=======================================================================
2560void AliAnalysisTaskFlowStrange::AddVZEQA() {
2561 fVZEQA = new TList();
2562 fVZEQA->SetName( Form("VZEQA%d",fRunNumber) );
2563 fVZEQA->SetOwner();
2564 if(fQAlevel>0) {
2565 TProfile2D *prof = new TProfile2D("LINP","LINP;VZEcell;VZEmult;SPDtrkl", 64,0,64,500,0,700,0,10000); fVZEQA->Add( prof );
2566 prof = new TProfile2D("MULP","MULP;VZEcell;CENTR;VZEmult", 64,0,64,100,0,100,0,10000); fVZEQA->Add( prof );
2567 TH3D *tH3D = new TH3D("EQU","EQU;VZEeqmult;VZEmult",100,0,700,100,0,700,64,0,64); fVZEQA->Add( tH3D );
2568 }
2569 fList->Add(fVZEQA);
2570}
2571//=======================================================================
2572void AliAnalysisTaskFlowStrange::SaveVZEROQA() {
2573 AliAODEvent *event = dynamic_cast<AliAODEvent*> (InputEvent());
2574 if(!event) return;
2575 AliVVZERO *vzero = event->GetVZEROData();
2576 AliAODTracklets *tracklets = event->GetTracklets();
2577 if(!vzero) return;
2578 if(!tracklets) return;
2579 Double_t mult, eqmult;
2580 Int_t trkl;
2581 for(int id=0; id!=64; ++id) {
2582 trkl = tracklets->GetNumberOfTracklets();
2583 mult = vzero->GetMultiplicity(id);
2584 eqmult = event->GetVZEROEqMultiplicity(id);
2585 ((TProfile2D*) fVZEQA->FindObject( "LINP" ))->Fill(id,mult,trkl,1);
2586 ((TProfile2D*) fVZEQA->FindObject( "MULP" ))->Fill(id,fThisCent,mult,1);
2587 ((TH3D*) fVZEQA->FindObject("EQU"))->Fill(eqmult,mult,id);
2588 }
2589}
2590//=======================================================================
2591void AliAnalysisTaskFlowStrange::AddVZEROResponse() {
2592 fVZEResponse = NULL;
2593 AliVEvent *event = InputEvent();
2594 if(!event) return;
2595 Int_t thisrun = event->GetRunNumber();
6c0d2d67 2596 fVZEResponse = new TH2D( Form("%d",thisrun), Form("%d;cell;CC",thisrun), 64,0,64, 110, -10, 100);
24373b38 2597 fList->Add(fVZEResponse);
2598}
2599//=======================================================================
2600void AliAnalysisTaskFlowStrange::SaveVZEROResponse() {
2601 if(!fVZEResponse) return;
2602 AliVEvent *event = InputEvent();
2603 if(!event) return;
2604 Double_t w;
2605 for(int id=0; id!=64; ++id) {
2606 w = event->GetVZEROEqMultiplicity(id);
2607 fVZEResponse->Fill(id,fThisCent,w);
2608 }
2609}
2610//=======================================================================
6c0d2d67 2611Int_t AliAnalysisTaskFlowStrange::RefMult(AliAODEvent *tAOD, Int_t fb) {
2612 Int_t found = 0;
2613 for(int i=0; i!=tAOD->GetNumberOfTracks(); ++i) {
2614 AliAODTrack *t = tAOD->GetTrack( i );
2615 if(!t) continue;
2616 if( !t->TestFilterBit(fb) ) continue;
2617 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2618 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2619 if( t->GetTPCNcls()<70 ) continue;
2620 //if( t->GetTPCsignal()<10.0 ) continue;
2621 if( t->Chi2perNDF()<0.2 ) continue;
2622 ++found;
2623 }
2624 return found;
2625}
2626//=======================================================================
24373b38 2627Int_t AliAnalysisTaskFlowStrange::RefMultTPC() {
2628 AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2629 if(!ev) return -1;
2630 Int_t found = 0;
2631 for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2632 AliAODTrack *t = ev->GetTrack( i );
2633 if(!t) continue;
2634 if( !t->TestFilterBit(1) ) continue;
2635 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2636 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2637 if( t->GetTPCNcls()<70 ) continue;
2638 if( t->GetTPCsignal()<10.0 ) continue;
2639 if( t->Chi2perNDF()<0.2 ) continue;
2640 ++found;
2641 }
2642 return found;
2643}
2644//=======================================================================
2645Int_t AliAnalysisTaskFlowStrange::RefMultGlobal() {
2646 AliAODEvent *ev = dynamic_cast<AliAODEvent*> (InputEvent());
2647 if(!ev) return -1;
2648 Int_t found = 0;
2649 for(int i=0; i!=ev->GetNumberOfTracks(); ++i) {
2650 AliAODTrack *t = ev->GetTrack( i );
2651 if(!t) continue;
2652 if( !t->TestFilterBit(16) ) continue;
2653 if( t->Eta()<-0.8 || t->Eta()>+0.8 ) continue;
2654 if( t->Pt()<0.2 || t->Pt()>5.0 ) continue;
2655 if( t->GetTPCNcls()<70 ) continue;
2656 if( t->GetTPCsignal()<10.0 ) continue;
2657 if( t->Chi2perNDF()<0.1 ) continue;
2658 Double_t b[3], bcov[3];
2659 if( !t->PropagateToDCA(ev->GetPrimaryVertex(),ev->GetMagneticField(),100,b,bcov) ) continue;
2660 if( b[0]>+0.3 || b[0]<-0.3 || b[1]>+0.3 || b[1]<-0.3) continue;
2661 ++found;
2662 }
2663 return found;
2664}
2665//=======================================================================
6c0d2d67 2666void AliAnalysisTaskFlowStrange::ResetContainers() {
2667 if(fUseFP) {
2668 fTPCevent->ClearFast();
2669 fVZEevent->ClearFast();
2670 }
2671}
2672//=======================================================================
2673void AliAnalysisTaskFlowStrange::AddTrackVn(TList *tList) {
2674 TProfile *tProfile;
2675 TH1D *tH1D;
2676 // ptbins
2677 Int_t npt = 29;
2678 Double_t ptbin[30] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2679 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.3, 2.6, 3.0, 3.5,
2680 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10., 12., 16., 20.};
2681 // vze
2682 tProfile = new TProfile("SP_uVZEA","u x Q_{VZEA}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2683 tProfile = new TProfile("SP_uVZEC","u x Q_{VZEC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2684 tProfile = new TProfile("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2685 tProfile = new TProfile("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2686 tProfile = new TProfile("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2687 // error
2688 tProfile = new TProfile("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2689 tProfile = new TProfile("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2690 tProfile = new TProfile("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2691 // tpc
2692 tProfile = new TProfile("SP_uTPCA","u x Q_{TPCA}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2693 tProfile = new TProfile("SP_uTPCC","u x Q_{TPCC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2694 tProfile = new TProfile("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2695 // error
2696 tProfile = new TProfile("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2697 tProfile = new TProfile("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2698 // control
2699 tH1D = new TH1D("QC_HistPt_P","HistPt_P",npt,ptbin); tList->Add( tH1D );
2700 tH1D = new TH1D("QC_HistPt_Q","HistPt_Q",npt,ptbin); tList->Add( tH1D );
2701 // qc
2702 tProfile = new TProfile("QC_C2","QC_C2",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2703 tProfile = new TProfile("QC_C4","QC_C4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2704 tProfile = new TProfile("QC_DC2","QC_DC2",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2705 tProfile = new TProfile("QC_DC4","QC_DC4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2706 tProfile = new TProfile("QC_C2C4","QC_C2C4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2707 tProfile = new TProfile("QC_C2DC2","QC_C2DC2",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2708 tProfile = new TProfile("QC_C2DC4","QC_C2DC4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2709 tProfile = new TProfile("QC_C4DC2","QC_C4DC2",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2710 tProfile = new TProfile("QC_C4DC4","QC_C4DC4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2711 tProfile = new TProfile("QC_DC2DC4","QC_DC2DC4",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2712 // qc transient
2713 tProfile = new TProfile("QC_pCos","QC_pCos",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2714 tProfile = new TProfile("QC_pSin","QC_pSin",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2715 tProfile = new TProfile("QC_qCos","QC_qCos",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2716 tProfile = new TProfile("QC_qSin","QC_qSin",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2717 tProfile = new TProfile("QC_q2hCos","QC_q2hCos",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2718 tProfile = new TProfile("QC_q2hSin","QC_q2hSin",npt,ptbin,-3,+3,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2719 // measurements
2720 tH1D = new TH1D("SP_vnVZEA","SP_vnVZEA",npt,ptbin); tList->Add( tH1D );
2721 tH1D = new TH1D("SP_vnVZEC","SP_vnVZEC",npt,ptbin); tList->Add( tH1D );
2722 tH1D = new TH1D("SP_vnTPCA","SP_vnTPCA",npt,ptbin); tList->Add( tH1D );
2723 tH1D = new TH1D("SP_vnTPCC","SP_vnTPCC",npt,ptbin); tList->Add( tH1D );
2724 tH1D = new TH1D("QC_Cum2","QC_Cum2",npt,ptbin); tList->Add( tH1D );
2725 tH1D = new TH1D("QC_Cum4","QC_Cum4",npt,ptbin); tList->Add( tH1D );
2726 tH1D = new TH1D("QC_DCum2","QC_DCum2",npt,ptbin); tList->Add( tH1D );
2727 tH1D = new TH1D("QC_DCum4","QC_DCum4",npt,ptbin); tList->Add( tH1D );
2728 tH1D = new TH1D("SP_vnVZEGA","SP_vnVZEGA",npt,ptbin); tList->Add( tH1D );
2729 tH1D = new TH1D("SP_vnVZEWA","SP_vnVZEWA",npt,ptbin); tList->Add( tH1D );
2730 tH1D = new TH1D("SP_vnTPCAA","SP_vnTPCAA",npt,ptbin); tList->Add( tH1D );
2731 tH1D = new TH1D("QC_vn2","QC_vn2",npt,ptbin); tList->Add( tH1D );
2732 tH1D = new TH1D("QC_vn4","QC_vn4",npt,ptbin); tList->Add( tH1D );
2733 if(fReadMC) {
2734 TH2D *tH2D;
2735 tProfile = new TProfile("MC_COSNDPHI","MC_COSNDPHI",npt,ptbin,-3,+3); tList->Add( tProfile );
2736 tH2D = new TH2D("MC_COSNDPHI_uQVZEA","MC_COSNDPHI_uQVZEA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2737 tH2D = new TH2D("MC_COSNDPHI_uQVZEC","MC_COSNDPHI_uQVZEC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2738 tH2D = new TH2D("MC_COSNDPHI_uQTPCA","MC_COSNDPHI_uQTPCA",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2739 tH2D = new TH2D("MC_COSNDPHI_uQTPCC","MC_COSNDPHI_uQTPCC",100,-1,+1,100,-0.3,+0.3); tList->Add( tH2D );
2740 }
2741}
2742//=======================================================================
2743void AliAnalysisTaskFlowStrange::AddDecayVn(TList *tList) {
2744 TProfile2D *tProfile;
2745 TH2D *tH2D;
2746 // ptbins
2747 Int_t npt = 29;
2748 Double_t ptbin[30] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2749 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.3, 2.6, 3.0, 3.5,
2750 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10., 12., 16., 20.};
2751 Int_t nms = 14;
2752 Double_t msbin[15] = { 0.400, 0.420, 0.440, 0.460, 0.480,
2753 0.490, 0.495, 0.500, 0.505, 0.510,
2754 0.520, 0.540, 0.560, 0.580, 0.600};
2755 // vze
2756 tProfile = new TProfile2D("SP_uVZEA","u x Q_{VZEA}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2757 tProfile = new TProfile2D("SP_uVZEC","u x Q_{VZEC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2758 tProfile = new TProfile2D("SP_VZEAVZEC","Q_{VZEA} x Q_{VZEC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2759 tProfile = new TProfile2D("SP_VZEATPC","Q_{VZEA} x Q_{TPC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2760 tProfile = new TProfile2D("SP_VZECTPC","Q_{VZEC} x Q_{TPC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2761 // error
2762 tProfile = new TProfile2D("SP_uVZEAuVZEC","u x Q_{VZEA} . u x Q_{VZEC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2763 tProfile = new TProfile2D("SP_uVZEAVZEAVZEC","u x Q_{VZEA} . Q_{VZEA} x Q_{VZEC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2764 tProfile = new TProfile2D("SP_uVZECVZEAVZEC","u x Q_{VZEC} . Q_{VZEA} x Q_{VZEC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2765 // tpc
2766 tProfile = new TProfile2D("SP_uTPCA","u x Q_{TPCA}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2767 tProfile = new TProfile2D("SP_uTPCC","u x Q_{TPCC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2768 tProfile = new TProfile2D("SP_TPCATPCC","Q_{TPCA} x Q_{TPCC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2769 // error
2770 tProfile = new TProfile2D("SP_uTPCATPCATPCC","u x Q_{TPCA} . Q_{TPCA} x Q_{TPCC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2771 tProfile = new TProfile2D("SP_uTPCCTPCATPCC","u x Q_{TPCC} . Q_{TPCA} x Q_{TPCC}",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2772 // control
2773 tH2D = new TH2D("QC_HistPt_P","HistPt_P",npt,ptbin,nms,msbin); tList->Add( tH2D );
2774 tH2D = new TH2D("QC_HistPt_Q","HistPt_Q",npt,ptbin,nms,msbin); tList->Add( tH2D );
2775 // qc
2776 tProfile = new TProfile2D("QC_C2","QC_C2",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2777 tProfile = new TProfile2D("QC_C4","QC_C4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2778 tProfile = new TProfile2D("QC_DC2","QC_DC2",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2779 tProfile = new TProfile2D("QC_DC4","QC_DC4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2780 tProfile = new TProfile2D("QC_C2C4","QC_C2C4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2781 tProfile = new TProfile2D("QC_C2DC2","QC_C2DC2",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2782 tProfile = new TProfile2D("QC_C2DC4","QC_C2DC4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2783 tProfile = new TProfile2D("QC_C4DC2","QC_C4DC2",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2784 tProfile = new TProfile2D("QC_C4DC4","QC_C4DC4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2785 tProfile = new TProfile2D("QC_DC2DC4","QC_DC2DC4",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2786 // qc transient
2787 tProfile = new TProfile2D("QC_pCos","QC_pCos",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2788 tProfile = new TProfile2D("QC_pSin","QC_pSin",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2789 tProfile = new TProfile2D("QC_qCos","QC_qCos",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2790 tProfile = new TProfile2D("QC_qSin","QC_qSin",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2791 tProfile = new TProfile2D("QC_q2hCos","QC_q2hCos",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2792 tProfile = new TProfile2D("QC_q2hSin","QC_q2hSin",npt,ptbin,nms,msbin,"s"); tList->Add( tProfile ); tProfile->Sumw2();
2793 // measurements
2794 tH2D = new TH2D("SP_vnVZEA","SP_vnVZEA",npt,ptbin,nms,msbin); tList->Add( tH2D );
2795 tH2D = new TH2D("SP_vnVZEC","SP_vnVZEC",npt,ptbin,nms,msbin); tList->Add( tH2D );
2796 tH2D = new TH2D("SP_vnTPCA","SP_vnTPCA",npt,ptbin,nms,msbin); tList->Add( tH2D );
2797 tH2D = new TH2D("SP_vnTPCC","SP_vnTPCC",npt,ptbin,nms,msbin); tList->Add( tH2D );
2798 tH2D = new TH2D("QC_Cum2","QC_Cum2",npt,ptbin,nms,msbin); tList->Add( tH2D );
2799 tH2D = new TH2D("QC_Cum4","QC_Cum4",npt,ptbin,nms,msbin); tList->Add( tH2D );
2800 tH2D = new TH2D("QC_DCum2","QC_DCum2",npt,ptbin,nms,msbin); tList->Add( tH2D );
2801 tH2D = new TH2D("QC_DCum4","QC_DCum4",npt,ptbin,nms,msbin); tList->Add( tH2D );
2802 tH2D = new TH2D("SP_vnVZEGA","SP_vnVZEGA",npt,ptbin,nms,msbin); tList->Add( tH2D );
2803 tH2D = new TH2D("SP_vnVZEWA","SP_vnVZEWA",npt,ptbin,nms,msbin); tList->Add( tH2D );
2804 tH2D = new TH2D("SP_vnTPCAA","SP_vnTPCAA",npt,ptbin,nms,msbin); tList->Add( tH2D );
2805 tH2D = new TH2D("QC_vn2","QC_vn2",npt,ptbin,nms,msbin); tList->Add( tH2D );
2806 tH2D = new TH2D("QC_vn4","QC_vn4",npt,ptbin,nms,msbin); tList->Add( tH2D );
2807 if(fReadMC) {
2808 tProfile = new TProfile2D("MC_COSNDPHI","MC_COSNDPHI",npt,ptbin,nms,msbin); tList->Add( tProfile );
2809 }
2810}
2811//=======================================================================
2812void AliAnalysisTaskFlowStrange::QCStoreTrackVn(TString name) {
2813 // getting transients
2814 TProfile *pCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ));
2815 TProfile *pSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ));
2816 TProfile *qCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ));
2817 TProfile *qSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ));
2818 TProfile *q2hCos = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ));
2819 TProfile *q2hSin = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ));
2820 if(fDebug) {
2821 printf("**QCStoreTrackVn( %s )\n",name.Data());
2822 printf(" Read %.0f entries in p set and %.0f entries in q set\n",pCos->GetEntries(),qCos->GetEntries());
2823 printf(" pCos[5] %.16f | pSin[5] %.16f \n", pCos->GetBinContent(5), pSin->GetBinContent(5));
2824 printf(" qCos[5] %.16f | qSin[5] %.16f \n", qCos->GetBinContent(5), qSin->GetBinContent(5));
2825 printf(" q2hCos[5] %.16f | q2hSin[5] %.16f \n", q2hCos->GetBinContent(5), q2hSin->GetBinContent(5));
2826 }
2827 // filling {2} and {4} correlator
2828 TProfile *c2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
2829 TProfile *c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
2830 TProfile *dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
2831 TProfile *dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
2832 TProfile *c2c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
2833 TProfile *c2dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
2834 TProfile *c2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
2835 TProfile *c4dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
2836 TProfile *c4dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
2837 TProfile *dc2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
2838 double tpc_qcos = fQTPCACos + fQTPCCCos;
2839 double tpc_qsin = fQTPCASin + fQTPCCSin;
2840 double tpc_qsqr = tpc_qcos*tpc_qcos + tpc_qsin*tpc_qsin;
2841 double tpc_q2hsqr = fQTPC2hCos*fQTPC2hCos + fQTPC2hSin*fQTPC2hSin;
2842 double tpc_qmul = fQTPCA + fQTPCC;
2843 int n = c2->GetNbinsX();
2844 for(int i=1; i!=n+1; ++i) {
2845 double mp = pCos->GetBinEntries( i );
2846 if( mp<1 ) { if(fDebug>2) printf(" bin %d:: mp (%.16f) < 1!\n",i,mp); continue; }
2847 double mm1 = tpc_qmul*(tpc_qmul-1);
2848 if( mm1<1e-100 ) { if(fDebug>2) printf(" bin %d:: mm1<1e-100!\n",i); continue; }
2849 double mq = qCos->GetBinEntries( i );
2850 double mpmmq = mp*tpc_qmul-mq;
2851 if( mpmmq<1e-100 ) { if(fDebug>2) printf(" bin %d:: mpmmq<1e-100!\n",i); continue; }
2852 double pt = c2->GetBinCenter( i );
2853 double pcos = pCos->GetBinContent( i )*mp;
2854 double psin = pSin->GetBinContent( i )*mp;
2855 double qcos = qCos->GetBinContent( i )*mq;
2856 double qsin = qSin->GetBinContent( i )*mq;
2857 double q2hcos = q2hCos->GetBinContent( i )*mq;
2858 double q2hsin = q2hSin->GetBinContent( i )*mq;
2859 double pQ = pcos*tpc_qcos+psin*tpc_qsin;
2860 double q2nQnQn = (qcos*tpc_qcos + qsin*tpc_qsin)*tpc_qcos + (qsin*tpc_qcos-qcos*tpc_qsin)*tpc_qsin;
2861 double pnQnQ2n = (pcos*tpc_qcos - psin*tpc_qsin)*fQTPC2hCos + (psin*tpc_qcos+pcos*tpc_qsin)*fQTPC2hSin;
2862 double tC2 = (tpc_qsqr-tpc_qmul)/mm1;
2863 double tDC2 = (pQ-mq)/mpmmq;
2864 c2->Fill( pt, tC2, mm1 );
2865 dc2->Fill( pt, tDC2, mpmmq );
2866 c2dc2->Fill( pt, tC2*tDC2, mm1*mpmmq );
2867 double mm1m2m3 = tpc_qmul*(tpc_qmul-1)*(tpc_qmul-2)*(tpc_qmul-3);
2868 if(mm1m2m3<1e-100) continue;
2869 double mpm3mqm1m2 = (mp*tpc_qmul-3*mq)*(tpc_qmul-1)*(tpc_qmul-2);
2870 if(mpm3mqm1m2<1e-100) continue;
2871 double req2hqnqn = fQTPC2hCos*(tpc_qcos*tpc_qcos-tpc_qsin*tpc_qsin)+2*fQTPC2hSin*tpc_qcos*tpc_qsin;
2872 double tC4 = (tpc_qsqr*tpc_qsqr + tpc_q2hsqr - 2*req2hqnqn - 2*(2*tpc_qsqr*(tpc_qmul-2)-tpc_qmul*(tpc_qmul-3)))/mm1m2m3;
2873 double tDC4 = pQ*tpc_qsqr -q2nQnQn -pnQnQ2n -2*(tpc_qmul-1)*pQ -2*mq*(tpc_qsqr-tpc_qmul+3) +6*(qcos*tpc_qcos+qsin*tpc_qsin) +(q2hcos*fQTPC2hCos+q2hsin*fQTPC2hSin);
2874 tDC4 /= mpm3mqm1m2;
2875 c4->Fill( pt, tC4, mm1m2m3 );
2876 dc4->Fill( pt, tDC4, mpm3mqm1m2 );
2877 c2c4->Fill( pt, tC2*tC4, mm1*mm1m2m3 );
2878 c2dc4->Fill( pt, tC2*tDC4, mm1*mpm3mqm1m2 );
2879 c4dc2->Fill( pt, tC4*tDC2, mm1m2m3*mpmmq );
2880 c4dc4->Fill( pt, tC4*tDC4, mm1m2m3*mpm3mqm1m2 );
2881 dc2dc4->Fill( pt, tDC2*tDC4, mpmmq*mpm3mqm1m2 );
2882 }
2883 // clean for next event
2884 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Reset();
2885 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Reset();
2886 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Reset();
2887 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Reset();
2888 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Reset();
2889 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Reset();
2890}
2891//=======================================================================
2892void AliAnalysisTaskFlowStrange::QCStoreDecayVn(TString name) {
2893 // getting transients
2894 TProfile2D *pCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ));
2895 TProfile2D *pSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ));
2896 TProfile2D *qCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ));
2897 TProfile2D *qSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ));
2898 TProfile2D *q2hCos = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ));
2899 TProfile2D *q2hSin = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ));
2900 if(fDebug) {
2901 printf("**QCStoreTrackVn( %s )\n",name.Data());
2902 printf(" Read %.0f entries in p set and %.0f entries in q set\n",pCos->GetEntries(),qCos->GetEntries());
2903 printf(" pCos[5][5] %.16f | pSin[5][5] %.16f \n", pCos->GetBinContent(5,5), pSin->GetBinContent(5,5));
2904 printf(" qCos[5][5] %.16f | qSin[5][5] %.16f \n", qCos->GetBinContent(5,5), qSin->GetBinContent(5,5));
2905 printf(" q2hCos[5][5] %.16f | q2hSin[5][5] %.16f \n", q2hCos->GetBinContent(5,5), q2hSin->GetBinContent(5,5));
2906 }
2907 // filling {2} and {4} correlator
2908 TProfile2D *c2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
2909 TProfile2D *c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
2910 TProfile2D *dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
2911 TProfile2D *dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
2912 TProfile2D *c2c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
2913 TProfile2D *c2dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
2914 TProfile2D *c2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
2915 TProfile2D *c4dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
2916 TProfile2D *c4dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
2917 TProfile2D *dc2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
2918 double tpc_qcos = fQTPCACos + fQTPCCCos;
2919 double tpc_qsin = fQTPCASin + fQTPCCSin;
2920 double tpc_qsqr = tpc_qcos*tpc_qcos + tpc_qsin*tpc_qsin;
2921 double tpc_q2hsqr = fQTPC2hCos*fQTPC2hCos + fQTPC2hSin*fQTPC2hSin;
2922 double tpc_qmul = fQTPCA + fQTPCC;
2923 int n = c2->GetNbinsX();
2924 int m = c2->GetNbinsY();
2925 for(int i=1; i!=n+1; ++i) {
2926 double pt = c2->GetXaxis()->GetBinCenter( i );
2927 for(int j=1; j!=m+1; ++j) {
2928 double ms = c2->GetYaxis()->GetBinCenter( j );
2929 int k = pCos->GetBin(i,j);
2930 double mp = pCos->GetBinEntries( k );
2931 if( mp<1 ) { if(fDebug>2) printf(" bin %d,%d:: mp (%.16f) < 1!\n",i,j,mp); continue; }
2932 double mm1 = tpc_qmul*(tpc_qmul-1);
2933 if( mm1<1e-100 ) { if(fDebug>2) printf(" bin %d,%d:: mm1<1e-100!\n",i,j); continue; }
2934 double mq = qCos->GetBinEntries( k );
2935 double mpmmq = mp*tpc_qmul-mq;
2936 if( mpmmq<1e-100 ) { if(fDebug>2) printf(" bin %d,%d:: mpmmq<1e-100!\n",i,j); continue; }
2937 double pcos = pCos->GetBinContent( i,j )*mp;
2938 double psin = pSin->GetBinContent( i,j )*mp;
2939 double qcos = qCos->GetBinContent( i,j )*mq;
2940 double qsin = qSin->GetBinContent( i,j )*mq;
2941 double q2hcos = q2hCos->GetBinContent( i,j )*mq;
2942 double q2hsin = q2hSin->GetBinContent( i,j )*mq;
2943 double pQ = pcos*tpc_qcos+psin*tpc_qsin;
2944 double q2nQnQn = (qcos*tpc_qcos + qsin*tpc_qsin)*tpc_qcos + (qsin*tpc_qcos-qcos*tpc_qsin)*tpc_qsin;
2945 double pnQnQ2n = (pcos*tpc_qcos - psin*tpc_qsin)*fQTPC2hCos + (psin*tpc_qcos+pcos*tpc_qsin)*fQTPC2hSin;
2946 double tC2 = (tpc_qsqr-tpc_qmul)/mm1;
2947 double tDC2 = (pQ-mq)/mpmmq;
2948 c2->Fill( pt, ms, tC2, mm1 );
2949 dc2->Fill( pt, ms, tDC2, mpmmq );
2950 c2dc2->Fill( pt, ms, tC2*tDC2, mm1*mpmmq );
2951 double mm1m2m3 = tpc_qmul*(tpc_qmul-1)*(tpc_qmul-2)*(tpc_qmul-3);
2952 if(mm1m2m3<1e-100) continue;
2953 double mpm3mqm1m2 = (mp*tpc_qmul-3*mq)*(tpc_qmul-1)*(tpc_qmul-2);
2954 if(mpm3mqm1m2<1e-100) continue;
2955 double req2hqnqn = fQTPC2hCos*(tpc_qcos*tpc_qcos-tpc_qsin*tpc_qsin)+2*fQTPC2hSin*tpc_qcos*tpc_qsin;
2956 double tC4 = (tpc_qsqr*tpc_qsqr + tpc_q2hsqr - 2*req2hqnqn - 2*(2*tpc_qsqr*(tpc_qmul-2)-tpc_qmul*(tpc_qmul-3)))/mm1m2m3;
2957 double tDC4 = pQ*tpc_qsqr -q2nQnQn -pnQnQ2n -2*(tpc_qmul-1)*pQ -2*mq*(tpc_qsqr-tpc_qmul+3) +6*(qcos*tpc_qcos+qsin*tpc_qsin) +(q2hcos*fQTPC2hCos+q2hsin*fQTPC2hSin);
2958 tDC4 /= mpm3mqm1m2;
2959 c4->Fill( pt, ms, tC4, mm1m2m3 );
2960 dc4->Fill( pt, ms, tDC4, mpm3mqm1m2 );
2961 c2c4->Fill( pt, ms, tC2*tC4, mm1*mm1m2m3 );
2962 c2dc4->Fill( pt, ms, tC2*tDC4, mm1*mpm3mqm1m2 );
2963 c4dc2->Fill( pt, ms, tC4*tDC2, mm1m2m3*mpmmq );
2964 c4dc4->Fill( pt, ms, tC4*tDC4, mm1m2m3*mpm3mqm1m2 );
2965 dc2dc4->Fill( pt, ms, tDC2*tDC4, mpmmq*mpm3mqm1m2 );
2966 }
2967 }
2968 // clean for next event
2969 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Reset();
2970 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Reset();
2971 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Reset();
2972 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Reset();
2973 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Reset();
2974 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Reset();
2975}
2976//=======================================================================
2977void AliAnalysisTaskFlowStrange::FillTrackVn(TString name,Double_t pt,Double_t phi,Double_t eta,Int_t fid) {
2978 // reading vze qm
2979 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
2980 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
2981 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
2982 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
2983 // reading tpc qm
2984 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
2985 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
2986 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
2987 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
2988 Double_t qtpc = fQTPCA+fQTPCC;
2989 Double_t tpc_qmcos = (fQTPCACos+fQTPCCCos)/qtpc;
2990 Double_t tpc_qmsin = (fQTPCASin+fQTPCCSin)/qtpc;
2991 // computing u
2992 Double_t cosn = TMath::Cos(fHarmonic*phi);
2993 Double_t sinn = TMath::Sin(fHarmonic*phi);
2994 Double_t cos2n = TMath::Cos(2.0*fHarmonic*phi);
2995 Double_t sin2n = TMath::Sin(2.0*fHarmonic*phi);
2996 // Scalar Product
2997 Double_t uQ, uQa, uQc, qaqc;
2998 // filling flow with vze
2999 qaqc = (vzea_qmcos*vzec_qmcos + vzea_qmsin*vzec_qmsin);
3000 uQa = (cosn*vzea_qmcos + sinn*vzea_qmsin);
3001 uQc = (cosn*vzec_qmcos + sinn*vzec_qmsin);
3002 Double_t cosmc = TMath::Cos( fHarmonic*GetMCDPHI(phi) );
3003 if(fReadMC) {
3004 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQVZEA" ))->Fill( cosmc,uQa );
3005 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQVZEC" ))->Fill( cosmc,uQc );
3006 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI" ))->Fill( pt,cosmc );
3007 }
3008 Double_t qaqt = (vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin);
3009 Double_t qcqt = (vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin);
3010 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" ))->Fill( pt,uQa,fQVZEA );
3011 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" ))->Fill( pt,uQc,fQVZEC );
3012 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" ))->Fill( pt,qaqc,fQVZEA*fQVZEC );
3013 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" ))->Fill( pt,qaqt,fQVZEA*qtpc );
3014 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" ))->Fill( pt,qcqt,fQVZEC*qtpc );
3015 // error vze
3016 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" ))->Fill( pt,uQa*uQc,fQVZEA*fQVZEC );
3017 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" ))->Fill( pt,uQa*qaqc,fQVZEA*fQVZEA*fQVZEC );
3018 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" ))->Fill( pt,uQc*qaqc,fQVZEC*fQVZEA*fQVZEC );
3019 // filling flow with tpc
3020 qaqc = (tpca_qmcos*tpcc_qmcos + tpca_qmsin*tpcc_qmsin);
3021 if(eta<0) {
3022 uQ = (cosn*tpca_qmcos + sinn*tpca_qmsin);
3023 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" ))->Fill( pt,uQ,fQTPCA );
3024 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" ))->Fill( pt,uQ*qaqc,fQTPCA*fQTPCA*fQTPCC );
3025 if(fReadMC) {
3026 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQTPCA" ))->Fill( cosmc,uQ );
3027 }
3028 } else {
3029 uQ = (cosn*tpcc_qmcos + sinn*tpcc_qmsin);
3030 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" ))->Fill( pt,uQ,fQTPCC );
3031 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" ))->Fill( pt,uQ*qaqc,fQTPCC*fQTPCA*fQTPCC );
3032 if(fReadMC) {
3033 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI_uQTPCC" ))->Fill( cosmc,uQ );
24373b38 3034 }
3035 }
6c0d2d67 3036 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" ))->Fill( pt,qaqc,fQTPCA*fQTPCC );
3037 // QC
3038 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_P" ))->Fill( pt );
3039 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Fill( pt, cosn, 1.0 );
3040 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Fill( pt, sinn, 1.0 );
3041 if(InQTPC(fid)) {
3042 ((TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_Q" ))->Fill( pt );
3043 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Fill( pt, cosn, 1.0 );
3044 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Fill( pt, sinn, 1.0 );
3045 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Fill( pt, cos2n, 1.0 );
3046 ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Fill( pt, sin2n, 1.0 );
3047 }
3048}
3049//=======================================================================
3050void AliAnalysisTaskFlowStrange::FillDecayVn(TString name,Double_t ms,Double_t pt,Double_t phi,Double_t eta,Int_t fid1,Int_t fid2) {
3051 // reading vze qm
3052 Double_t vzec_qmcos = fQVZECCos/fQVZEC;
3053 Double_t vzec_qmsin = fQVZECSin/fQVZEC;
3054 Double_t vzea_qmcos = fQVZEACos/fQVZEA;
3055 Double_t vzea_qmsin = fQVZEASin/fQVZEA;
3056 // reading tpc qm
3057 Double_t tpcc_qmcos = fQTPCCCos/fQTPCC;
3058 Double_t tpcc_qmsin = fQTPCCSin/fQTPCC;
3059 Double_t tpca_qmcos = fQTPCACos/fQTPCA;
3060 Double_t tpca_qmsin = fQTPCASin/fQTPCA;
3061 Double_t qtpc = fQTPCA+fQTPCC;
3062 Double_t tpc_qmcos = (fQTPCACos+fQTPCCCos)/qtpc;
3063 Double_t tpc_qmsin = (fQTPCASin+fQTPCCSin)/qtpc;
3064 // computing u
3065 Double_t cosn = TMath::Cos(fHarmonic*phi);
3066 Double_t sinn = TMath::Sin(fHarmonic*phi);
3067 Double_t cos2n = TMath::Cos(2.0*fHarmonic*phi);
3068 Double_t sin2n = TMath::Sin(2.0*fHarmonic*phi);
3069 // Scalar Product
3070 Double_t uQ, uQa, uQc, qaqc;
3071 // filling flow with vze
3072 qaqc = (vzea_qmcos*vzec_qmcos + vzea_qmsin*vzec_qmsin);
3073 uQa = (cosn*vzea_qmcos + sinn*vzea_qmsin);
3074 uQc = (cosn*vzec_qmcos + sinn*vzec_qmsin);
3075 Double_t cosmc = TMath::Cos( fHarmonic*GetMCDPHI(phi) );
3076 if(fReadMC) {
3077 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "MC_COSNDPHI" ))->Fill( pt,ms,cosmc );
3078 }
3079 Double_t qaqt = (vzea_qmcos*tpc_qmcos + vzea_qmsin*tpc_qmsin);
3080 Double_t qcqt = (vzec_qmcos*tpc_qmcos + vzec_qmsin*tpc_qmsin);
3081 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" ))->Fill( pt,ms,uQa,fQVZEA );
3082 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" ))->Fill( pt,ms,uQc,fQVZEC );
3083 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" ))->Fill( pt,ms,qaqc,fQVZEA*fQVZEC );
3084 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" ))->Fill( pt,ms,qaqt,fQVZEA*qtpc );
3085 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" ))->Fill( pt,ms,qcqt,fQVZEC*qtpc );
3086 // error vze
3087 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" ))->Fill( pt,ms,uQa*uQc,fQVZEA*fQVZEC );
3088 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" ))->Fill( pt,ms,uQa*qaqc,fQVZEA*fQVZEA*fQVZEC );
3089 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" ))->Fill( pt,ms,uQc*qaqc,fQVZEC*fQVZEA*fQVZEC );
3090 // filling flow with tpc
3091 qaqc = (tpca_qmcos*tpcc_qmcos + tpca_qmsin*tpcc_qmsin);
3092 if(eta<0) {
3093 uQ = (cosn*tpca_qmcos + sinn*tpca_qmsin);
3094 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" ))->Fill( pt,ms,uQ,fQTPCA );
3095 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" ))->Fill( pt,ms,uQ*qaqc,fQTPCA*fQTPCA*fQTPCC );
3096 } else {
3097 uQ = (cosn*tpcc_qmcos + sinn*tpcc_qmsin);
3098 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" ))->Fill( pt,ms,uQ,fQTPCC );
3099 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" ))->Fill( pt,ms,uQ*qaqc,fQTPCC*fQTPCA*fQTPCC );
3100 }
3101 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" ))->Fill( pt,ms,qaqc,fQTPCA*fQTPCC );
3102 // QC
3103 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_P" ))->Fill( pt,ms );
3104 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pCos" ))->Fill( pt, ms, cosn, 1.0 );
3105 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_pSin" ))->Fill( pt, ms, sinn, 1.0 );
3106 if(InQTPC(fid1)||InQTPC(fid2)) {
3107 ((TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_HistPt_Q" ))->Fill( pt,ms );
3108 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qCos" ))->Fill( pt, ms, cosn, 1.0 );
3109 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_qSin" ))->Fill( pt, ms, sinn, 1.0 );
3110 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hCos" ))->Fill( pt, ms, cos2n, 1.0 );
3111 ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_q2hSin" ))->Fill( pt, ms, sin2n, 1.0 );
3112 }
3113}
3114//=======================================================================
3115Bool_t AliAnalysisTaskFlowStrange::InQTPC(Int_t id) {
3116 Bool_t ret = kFALSE;
3117 for(int i=0; i!=fQTPCA_nTracks; ++i)
3118 if(fQTPCA_fID[i]==id) {
3119 ret=kTRUE;
3120 break;
3121 }
3122 if(ret) return kTRUE;
3123 for(int i=0; i!=fQTPCC_nTracks; ++i)
3124 if(fQTPCC_fID[i]==id) {
3125 ret=kTRUE;
3126 break;
3127 }
3128 return ret;
24373b38 3129}
3130//=======================================================================
6c0d2d67 3131void AliAnalysisTaskFlowStrange::ComputeTrackVn(TString name) {
3132 TProfile *uQa, *uQc, *qaqc, *uQaqaqc, *uQcqaqc;
3133 TArrayD *pasww, *pbsww, *pcsww;
3134 //ScalarProducr TPC
3135 printf("<<%s>> SP TPC\n",name.Data());
3136 uQa = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" );
3137 uQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" );
3138 qaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" );
3139 uQaqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" );
3140 uQcqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" );
3141 pasww = uQa->GetBinSumw2();
3142 pbsww = uQc->GetBinSumw2();
3143 pcsww = qaqc->GetBinSumw2();
3144 //
3145 TH1D *sptpca = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCA" );
3146 TH1D *sptpcc = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCC" );
3147 TH1D *sptpcaa = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCAA" );
3148 //
3149 for(Int_t i=1; i!=sptpcaa->GetNbinsX()+1; ++i) {
3150 sptpcaa->SetBinContent(i,0);
3151 sptpcaa->SetBinError(i,0);
3152 sptpca->SetBinContent(i,0);
3153 sptpca->SetBinError(i,0);
3154 sptpcc->SetBinContent(i,0);
3155 sptpcc->SetBinError(i,0);
3156 double a = uQa->GetBinContent(i);
3157 double b = uQc->GetBinContent(i);
3158 double c = qaqc->GetBinContent(i);
3159 //if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3160 //if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3161 if(c<1e-100) continue;
3162 // nominal sptpca
3163 double vna = a/TMath::Sqrt(c);
3164 sptpca->SetBinContent(i,vna);
3165 // nominal sptpcc
3166 double vnc = b/TMath::Sqrt(c);
3167 sptpcc->SetBinContent(i,vnc);
3168 // nominal sptpc
3169 double vn = (vna + vnc)/2.0;
3170 sptpcaa->SetBinContent(i,vn);
3171 // errors
3172 double asw = uQa->GetBinEntries(i);
3173 double bsw = uQc->GetBinEntries(i);
3174 double csw = qaqc->GetBinEntries(i);
3175 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3176 double asww = pasww->At(i);
3177 double bsww = pbsww->At(i);
3178 double csww = pcsww->At(i);
3179 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3180 if((1<1e-100+asww/asw/asw)||(1<1e-100+bsww/bsw/bsw)||(1<1e-100+csww/csw/csw)) continue;
3181 if(TMath::AreEqualAbs(asww,asw*asw,1e-100)||
3182 TMath::AreEqualAbs(bsww,bsw*bsw,1e-100)||
3183 TMath::AreEqualAbs(csww,csw*csw,1e-100)) continue;
3184 double ac = uQaqaqc->GetBinContent(i);
3185 double bc = uQcqaqc->GetBinContent(i);
3186 double acsw = uQaqaqc->GetBinEntries(i);
3187 double bcsw = uQcqaqc->GetBinEntries(i);
3188 double ea = uQa->GetBinError(i)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3189 double eb = uQc->GetBinError(i)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3190 double ec = qaqc->GetBinError(i)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3191 //printf("%d >> ea^2 %.16f |||| asww %.16f | asw %.16f | 1-asww/asw/asw %.16f \n", i,ea, asww, asw, 1-asww/asw/asw);
3192 //printf("%d >> eb^2 %.16f |||| bsww %.16f | bsw %.16f | 1-bsww/bsw/bsw %.16f \n", i,eb, bsww, bsw, 1-bsww/bsw/bsw);
3193 //printf("%d >> ec^2 %.16f |||| csww %.16f | csw %.16f | 1-csww/csw/csw %.16f \n", i,ec, csww, csw, 1-csww/csw/csw);
3194 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3195 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3196 double evna = 1.0/TMath::Abs(c) * ( ea*ea + vna*vna/TMath::Abs(c)/4.0*ec*ec - vna/TMath::Sqrt(c)*eac );
3197 double evnc = 1.0/TMath::Abs(c) * ( eb*eb + vnc*vnc/TMath::Abs(c)/4.0*ec*ec - vnc/TMath::Sqrt(c)*ebc );
3198 //printf("%d >> evna^2 %.16f |||| ea %.16f | ec %.16f | eac %.16f | c %.16f\n", i,evna, ea, ec, eac,c);
3199 //printf("%d >> evnc^2 %.16f |||| eb %.16f | ec %.16f | ebc %.16f | c %.16f\n", i,evnc, eb, ec, ebc,c);
3200 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3201 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3202 sptpca->SetBinError(i,evna);
3203 sptpcc->SetBinError(i,evnc);
3204 sptpcaa->SetBinError(i,TMath::Sqrt(evna*evna+evnc*evnc)/2.0);
3205 }
3206 //ScalarProduct VZE
3207 printf("<<%s>> SP VZE\n",name.Data());
3208 double cvzea2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->GetBinContent( 1 );
3209 double cvzec2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->GetBinContent( 1 );
3210 if( TMath::AreEqualAbs(cvzea2+cvzec2,0,1e-100) ) return;
3211 uQa = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" );
3212 uQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" );
3213 qaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" );
3214 uQaqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" );
3215 uQcqaqc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" );
3216 pasww = uQa->GetBinSumw2();
3217 pbsww = uQc->GetBinSumw2();
3218 pcsww = qaqc->GetBinSumw2();
3219 //
3220 TProfile *qaqt = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" );
3221 TProfile *qcqt = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" );
3222 TProfile *uQauQc = (TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" );
3223 //
3224 TH1D *spvzea = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEA" );
3225 TH1D *spvzec = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEC" );
3226 TH1D *spvzega = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEGA" );
3227 TH1D *spvzewa = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEWA" );
3228 for(Int_t i=1; i!=spvzewa->GetNbinsX()+1; ++i) {
3229 spvzega->SetBinContent(i,0);
3230 spvzega->SetBinError(i,0);
3231 spvzewa->SetBinContent(i,0);
3232 spvzewa->SetBinError(i,0);
3233 spvzea->SetBinContent(i,0);
3234 spvzea->SetBinError(i,0);
3235 spvzec->SetBinContent(i,0);
3236 spvzec->SetBinError(i,0);
3237 double a = uQa->GetBinContent(i);
3238 double b = uQc->GetBinContent(i);
3239 double c = qaqc->GetBinContent(i);
3240 double at = qaqt->GetBinContent(i);
3241 double bt = qcqt->GetBinContent(i);
3242 if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3243 if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3244 if(TMath::AreEqualAbs(c,0,1e-100)) continue;
3245 if(TMath::AreEqualAbs(at,0,1e-100)) continue;
3246 if(TMath::AreEqualAbs(bt,0,1e-100)) continue;
3247 // nominal spvzea
3248 double aa = c*at/bt;
3249 if(aa<1e-100) continue;
3250 double vna = a/TMath::Sqrt(aa);
3251 spvzea->SetBinContent(i,vna);
3252 // nominal spvzec
3253 double bb = c*bt/at;
3254 if(bb<1e-100) continue;
3255 double vnc = b/TMath::Sqrt(bb);
3256 spvzec->SetBinContent(i,vnc);
3257 //nominal spvzewa
3258 double vnwa = (cvzea2*vna + cvzec2*vnc) / (cvzea2+cvzec2);
3259 spvzewa->SetBinContent(i,vnwa);
3260 // nominal spvzega
3261 double vnga = a*b/c;
3262 if(vnga<1e-100) continue;
3263 vnga = TMath::Sqrt(vnga);
3264 spvzega->SetBinContent(i,vnga);
3265 // errors
3266 double asw = uQa->GetBinEntries(i);
3267 double bsw = uQc->GetBinEntries(i);
3268 double csw = qaqc->GetBinEntries(i);
3269 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3270 double asww = pasww->At(i);
3271 double bsww = pbsww->At(i);
3272 double csww = pcsww->At(i);
3273 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3274 if((1<asww/asw/asw)||(1<bsww/bsw/bsw)||(1<csww/csw/csw)) continue;
3275 double ab = uQauQc->GetBinContent(i);
3276 double ac = uQaqaqc->GetBinContent(i);
3277 double bc = uQcqaqc->GetBinContent(i);
3278 double absw = uQauQc->GetBinEntries(i);
3279 double acsw = uQaqaqc->GetBinEntries(i);
3280 double bcsw = uQcqaqc->GetBinEntries(i);
3281 if(TMath::AreEqualAbs(1,absw/asw/bsw,1e-100)||TMath::AreEqualAbs(1,bcsw/bsw/csw,1e-100)||TMath::AreEqualAbs(1,acsw/asw/csw,1e-100)) continue;
3282 double ea = uQa->GetBinError(i)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3283 double eb = uQc->GetBinError(i)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3284 double ec = qaqc->GetBinError(i)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3285 double eab = (ab-a*b)/(1-absw/asw/bsw)*absw/asw/bsw;
3286 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3287 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3288 double nc, nec, neac, nebc;
3289 nc = c*at/bt;
3290 nec = ec*at/bt;
3291 neac = eac*at/bt;
3292 nebc = ebc*at/bt;
3293 double evna = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*ea*ea + a*a/4.0*nec*nec - a*TMath::Abs(nc)*neac*neac );
3294 nc = c*bt/at;
3295 nec = ec*bt/at;
3296 neac = eac*bt/at;
3297 nebc = ebc*bt/at;
3298 double evnc = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*eb*eb + b*b/4.0*nec*nec - b*TMath::Abs(nc)*nebc*nebc );
3299 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3300 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3301 spvzea->SetBinError(i,evna);
3302 spvzec->SetBinError(i,evnc);
3303 double evnwa = TMath::Sqrt( cvzea2*cvzea2*evna*evna + cvzec2*cvzec2*evnc*evnc )/(cvzea2+cvzec2);
3304 spvzewa->SetBinError(i,evnwa);
3305 double evnga = 0.25/c/c * ( TMath::Abs(b*c/a)*ea*ea + TMath::Abs(a*c/b)*eb*eb + TMath::Abs(a*b/c)*ec*ec + 2*c*eab - 2*a*ebc - 2*b*eac );
3306 if(evnga>1e-100) evnga = TMath::Sqrt(evnga); else evnga=0;
3307 spvzega->SetBinError(i,evnga);
3308 }
3309 printf("<<%s>> QC TPC\n",name.Data());
3310 //Qcumulants
3311 TH1D *resC2 = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum2" );
3312 TH1D *resC4 = (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum4" );
3313 TH1D *resDC2= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum2" );
3314 TH1D *resDC4= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum4" );
3315 TH1D *resvn2= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn2" );
3316 TH1D *resvn4= (TH1D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn4" );
3317 //correlators
3318 TProfile *c2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
3319 TProfile *c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
3320 TProfile *dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
3321 TProfile *dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
3322 TProfile *c2c4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
3323 TProfile *c2dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
3324 TProfile *c2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
3325 TProfile *c4dc2 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
3326 TProfile *c4dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
3327 TProfile *dc2dc4 = ((TProfile*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
3328 TArrayD *c2sww = c2->GetBinSumw2();
3329 TArrayD *c4sww = c4->GetBinSumw2();
3330 TArrayD *dc2sww= dc2->GetBinSumw2();
3331 TArrayD *dc4sww= dc4->GetBinSumw2();
3332 for(Int_t i=1; i!=resvn2->GetNbinsX()+1; ++i) {
3333 // cn{2}
3334 double v_c2sw = c2->GetBinEntries(i);
3335 if(v_c2sw<1e-100) continue;
3336 double v_c2sww = c2sww->At(i);
3337 double v_c2 = c2->GetBinContent(i);
3338 double e_c2 = TMath::Sqrt(v_c2sww)/v_c2sw*c2->GetBinError(i);
3339 double cum2 = v_c2;
3340 double ecum2= e_c2;
3341 resC2->SetBinContent(i, cum2 );
3342 resC2->SetBinError(i, ecum2 );
3343 // cn{4}
3344 double v_c4sw = c4->GetBinEntries(i);
3345 if(v_c4sw<1e-100) continue;
3346 double v_c4sww = c4sww->At(i);
3347 double v_c4 = c4->GetBinContent(i);
3348 double e_c4 = TMath::Sqrt(v_c4sww)/v_c4sw*c4->GetBinError(i);
3349 double v_c2c4 = c2c4->GetBinContent(i);
3350 double v_c2c4sw = c2c4->GetBinEntries(i);
3351 if(TMath::AreEqualAbs(v_c2c4sw/v_c2sw/v_c4sw,1,1e-100)) continue;
3352 double covc2c4 = v_c2c4sw/v_c2sw/v_c4sw*(v_c2c4 - v_c2*v_c4)/(1-v_c2c4sw/v_c2sw/v_c4sw);
3353 double cum4 = v_c4 - 2*v_c2*v_c2;
3354 double ecum4= 16.0*v_c2*v_c2*e_c2*e_c2 + e_c4*e_c4 - 8.0*v_c2*covc2c4;
3355 if(ecum4<1e-100) continue;
3356 ecum4 = TMath::Sqrt( ecum4 );
3357 resC4->SetBinContent(i, cum4 );
3358 resC4->SetBinError(i, ecum4 );
3359 // dn{2}
3360 double v_dc2sw = dc2->GetBinEntries(i);
3361 if(v_dc2sw<1) continue;
3362 double v_dc2 = dc2->GetBinContent(i);
3363 double v_dc2sww = dc2sww->At(i);
3364 double e_dc2 = TMath::Sqrt(v_dc2sww)/v_dc2sw*dc2->GetBinError(i);
3365 double dcum2 = v_dc2;
3366 double edcum2= e_dc2;
3367 resDC2->SetBinContent(i, dcum2 );
3368 resDC2->SetBinError(i, edcum2 );
3369 // v2{2}
3370 if(v_c2<1e-100) continue;
3371 double dv22 = v_dc2/TMath::Sqrt(v_c2);
3372 double v_c2dc2 = c2dc2->GetBinContent(i);
3373 double v_c2dc2sw = c2dc2->GetBinEntries(i);
3374 if(TMath::AreEqualAbs(v_c2dc2sw/v_c2sw/v_dc2sw,1,1e-100)) continue;
3375 double covc2dc2 = v_c2dc2sw/v_c2sw/v_dc2sw*(v_c2dc2 - v_c2*v_dc2)/(1-v_c2dc2sw/v_c2sw/v_dc2sw);
3376 double edv22 = 0.25/v_c2/v_c2/v_c2*(v_dc2*v_dc2*e_c2*e_c2 + 4*v_c2*v_c2*e_dc2*e_dc2 - 4*v_c2*v_dc2*covc2dc2);
3377 //printf("%d >> dv22 %.16f || edv22^2 %.16f |||| v_c2dc2 %.16f | v_c2dc2sw %.16f | covc2dc2 %.16f \n", i,dv22,edv22,v_c2dc2,v_c2dc2sw,covc2dc2);
3378 if(edv22<1e-100) continue;
3379 edv22 = TMath::Sqrt(edv22);
3380 resvn2->SetBinContent(i,dv22);
3381 resvn2->SetBinError(i,edv22);
3382 // dn{4}
3383 double v_dc4sw = dc4->GetBinEntries(i);
3384 if(v_dc4sw<1) continue;
3385 double v_dc4 = dc4->GetBinContent(i);
3386 double v_dc4sww = dc4sww->At(i);
3387 double e_dc4 = TMath::Sqrt(v_dc4sww)/v_dc4sw*dc4->GetBinError(i);
3388 double dcum4 = v_dc4 - 2*v_c2*v_dc2;
3389 double v_c2dc4 = c2dc4->GetBinContent(i);
3390 double v_c2dc4sw = c2dc4->GetBinEntries(i);
3391 if(TMath::AreEqualAbs(v_c2dc4sw/v_c2sw/v_dc4sw,1,1e-100)) continue;
3392 double covc2dc4 = v_c2dc4sw/v_c2sw/v_dc4sw*(v_c2dc4 - v_c2*v_dc4)/(1-v_c2dc4sw/v_c2sw/v_dc4sw);
3393 double v_dc2dc4 = dc2dc4->GetBinContent(i);
3394 double v_dc2dc4sw = dc2dc4->GetBinEntries(i);
3395 if(TMath::AreEqualAbs(v_dc2dc4sw/v_dc2sw/v_dc4sw,1,1e-100)) continue;
3396 double covdc2dc4 = v_dc2dc4sw/v_dc2sw/v_dc4sw*(v_dc2dc4 - v_dc2*v_dc4)/(1-v_dc2dc4sw/v_dc2sw/v_dc4sw);
3397 double edcum4= ( +4.0*v_dc2*v_dc2*e_c2*e_c2
3398 +4.0*v_c2*v_c2*e_dc2*e_dc2
3399 +e_dc4*e_dc4
3400 +8.0*v_c2*v_dc2*covc2dc2
3401 -4.0*v_dc2*covc2dc4
3402 -4.0*v_c2*covdc2dc4 );
3403 if(edcum4<1e-100) continue;
3404 edcum4 = TMath::Sqrt(edcum4);
3405 resDC4->SetBinContent(i, dcum4 );
3406 resDC4->SetBinError(i, edcum4 );
3407 // v2{4}
3408 if(cum4>1e-100) continue;
3409 double dv24 = -dcum4/TMath::Power(-cum4,0.75);
3410 double dterm1 = 2*v_c2*v_c2*v_dc2 - 3*v_c2*v_dc4 + 2*v_c4*v_dc2;
3411 double dterm2 = 9.0/16.0*dcum4*dcum4;
3412 double dterm3 = 4.0*v_c2*v_c2*cum4*cum4;
3413 double dterm4 = cum4*cum4;
3414 double dterm5 = -3.0/2.0*dcum4*dterm1;
3415 double dterm6 = -4.0*v_c2*cum4*dterm1;
3416 double dterm7 = -2.0*cum4*dterm1;
3417 double dterm8 = 3.0*v_c2*cum4*dcum4;
3418 double dterm9 = 3.0/2.0*cum4*dcum4;
3419 double dterm10= 4*v_c2*cum4*cum4;
3420 double v_c4dc2 = c4dc2->GetBinContent(i);
3421 double v_c4dc2sw = c4dc2->GetBinEntries(i);
3422 if(TMath::AreEqualAbs(v_c4dc2sw/v_c4sw/v_dc2sw,1,1e-100)) continue;
3423 double covc4dc2 = v_c4dc2sw/v_c4sw/v_dc2sw*(v_c4dc2 - v_c4*v_dc2)/(1-v_c4dc2sw/v_c4sw/v_dc2sw);
3424 double v_c4dc4 = c4dc4->GetBinContent(i);
3425 double v_c4dc4sw = c4dc4->GetBinEntries(i);
3426 if(TMath::AreEqualAbs(v_c4dc4sw/v_c4sw/v_dc4sw,1,1e-100)) continue;
3427 double covc4dc4 = v_c4dc4sw/v_c4sw/v_dc4sw*(v_c4dc4 - v_c4*v_dc4)/(1-v_c4dc4sw/v_c4sw/v_dc4sw);
3428 double edv24= 1.0/TMath::Power(-cum4,3.5)*(+dterm1*dterm1*e_c2*e_c2
3429 +dterm2*e_c4*e_c4
3430 +dterm3*e_dc2*e_dc2
3431 +dterm4*e_dc4*e_dc4
3432 -dterm5*covc2c4
3433 -dterm6*covc2dc2
3434 +dterm7*covc2dc4
3435 +dterm8*covc4dc2
3436 -dterm9*covc4dc4
3437 -dterm10*covdc2dc4);
3438 if(edv24<1e-100) continue;
3439 edv24 = TMath::Sqrt(edv24);
3440 resvn4->SetBinContent(i,dv24);
3441 resvn4->SetBinError(i,edv24);
3442 }
3443}
3444//=======================================================================
3445void AliAnalysisTaskFlowStrange::ComputeDecayVn(TString name) {
3446 TProfile2D *uQa, *uQc, *qaqc, *uQaqaqc, *uQcqaqc;
3447 TArrayD *pasww, *pbsww, *pcsww;
3448 //ScalarProducr TPC
3449 printf("<<%s>> SP TPC\n",name.Data());
3450 uQa = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCA" );
3451 uQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCC" );
3452 qaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_TPCATPCC" );
3453 uQaqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCATPCATPCC" );
3454 uQcqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uTPCCTPCATPCC" );
3455 pasww = uQa->GetBinSumw2();
3456 pbsww = uQc->GetBinSumw2();
3457 pcsww = qaqc->GetBinSumw2();
3458 //
3459 TH2D *sptpca = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCA" );
3460 TH2D *sptpcc = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCC" );
3461 TH2D *sptpcaa = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnTPCAA" );
3462 //
3463 for(Int_t i=1; i!=sptpcaa->GetNbinsX()+1; ++i) {
3464 for(Int_t j=1; j!=sptpcaa->GetNbinsY()+1; ++j) {
3465 sptpcaa->SetBinContent(i,j,0);
3466 sptpcaa->SetBinError(i,j,0);
3467 sptpca->SetBinContent(i,j,0);
3468 sptpca->SetBinError(i,j,0);
3469 sptpcc->SetBinContent(i,j,0);
3470 sptpcc->SetBinError(i,j,0);
3471 double a = uQa->GetBinContent(i,j);
3472 double b = uQc->GetBinContent(i,j);
3473 double c = qaqc->GetBinContent(i,j);
3474 //if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3475 //if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3476 if(c<1e-100) continue;
3477 // nominal sptpca
3478 double vna = a/TMath::Sqrt(c);
3479 sptpca->SetBinContent(i,vna);
3480 // nominal sptpcc
3481 double vnc = b/TMath::Sqrt(c);
3482 sptpcc->SetBinContent(i,vnc);
3483 // nominal sptpc
3484 double vn = (vna + vnc)/2.0;
3485 sptpcaa->SetBinContent(i,vn);
3486 // errors
3487 int k = sptpcaa->GetBin(i,j);
3488 double asw = uQa->GetBinEntries(k);
3489 double bsw = uQc->GetBinEntries(k);
3490 double csw = qaqc->GetBinEntries(k);
3491 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3492 double asww = pasww->At(k);
3493 double bsww = pbsww->At(k);
3494 double csww = pcsww->At(k);
3495 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3496 if((1<1e-100+asww/asw/asw)||(1<1e-100+bsww/bsw/bsw)||(1<1e-100+csww/csw/csw)) continue;
3497 if(TMath::AreEqualAbs(asww,asw*asw,1e-100)||
3498 TMath::AreEqualAbs(bsww,bsw*bsw,1e-100)||
3499 TMath::AreEqualAbs(csww,csw*csw,1e-100)) continue;
3500 double ac = uQaqaqc->GetBinContent(i,j);
3501 double bc = uQcqaqc->GetBinContent(i,j);
3502 double acsw = uQaqaqc->GetBinEntries(k);
3503 double bcsw = uQcqaqc->GetBinEntries(k);
3504 double ea = uQa->GetBinError(i,j)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3505 double eb = uQc->GetBinError(i,j)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3506 double ec = qaqc->GetBinError(i,j)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3507 //printf("%d >> ea^2 %.16f |||| asww %.16f | asw %.16f | 1-asww/asw/asw %.16f \n", i,ea, asww, asw, 1-asww/asw/asw);
3508 //printf("%d >> eb^2 %.16f |||| bsww %.16f | bsw %.16f | 1-bsww/bsw/bsw %.16f \n", i,eb, bsww, bsw, 1-bsww/bsw/bsw);
3509 //printf("%d >> ec^2 %.16f |||| csww %.16f | csw %.16f | 1-csww/csw/csw %.16f \n", i,ec, csww, csw, 1-csww/csw/csw);
3510 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3511 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3512 double evna = 1.0/TMath::Abs(c) * ( ea*ea + vna*vna/TMath::Abs(c)/4.0*ec*ec - vna/TMath::Sqrt(c)*eac );
3513 double evnc = 1.0/TMath::Abs(c) * ( eb*eb + vnc*vnc/TMath::Abs(c)/4.0*ec*ec - vnc/TMath::Sqrt(c)*ebc );
3514 //printf("%d >> evna^2 %.16f |||| ea %.16f | ec %.16f | eac %.16f | c %.16f\n", i,evna, ea, ec, eac,c);
3515 //printf("%d >> evnc^2 %.16f |||| eb %.16f | ec %.16f | ebc %.16f | c %.16f\n", i,evnc, eb, ec, ebc,c);
3516 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3517 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3518 sptpca->SetBinError(i,j,evna);
3519 sptpcc->SetBinError(i,j,evnc);
3520 sptpcaa->SetBinError(i,j,TMath::Sqrt(evna*evna+evnc*evnc)/2.0);
3521 }
3522 }
3523 //ScalarProduct VZE
3524 printf("<<%s>> SP VZE\n",name.Data());
3525 double cvzea2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEA"))->GetBinContent( 1 );
3526 double cvzec2 = ((TH1D*)((TList*)fList->FindObject("MakeQSpy"))->FindObject("ChiSquaredVZEC"))->GetBinContent( 1 );
3527 if( TMath::AreEqualAbs(cvzea2+cvzec2,0,1e-100) ) return;
3528 uQa = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEA" );
3529 uQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEC" );
3530 qaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEAVZEC" );
3531 uQaqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAVZEAVZEC" );
3532 uQcqaqc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZECVZEAVZEC" );
3533 pasww = uQa->GetBinSumw2();
3534 pbsww = uQc->GetBinSumw2();
3535 pcsww = qaqc->GetBinSumw2();
3536 //
3537 TProfile2D *qaqt = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZEATPC" );
3538 TProfile2D *qcqt = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_VZECTPC" );
3539 TProfile2D *uQauQc = (TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_uVZEAuVZEC" );
3540 //
3541 TH2D *spvzea = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEA" );
3542 TH2D *spvzec = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEC" );
3543 TH2D *spvzega = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEGA" );
3544 TH2D *spvzewa = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "SP_vnVZEWA" );
3545 for(Int_t i=1; i!=spvzewa->GetNbinsX()+1; ++i) {
3546 for(Int_t j=1; j!=spvzewa->GetNbinsX()+1; ++j) {
3547 spvzega->SetBinContent(i,j,0);
3548 spvzega->SetBinError(i,j,0);
3549 spvzewa->SetBinContent(i,j,0);
3550 spvzewa->SetBinError(i,j,0);
3551 spvzea->SetBinContent(i,j,0);
3552 spvzea->SetBinError(i,j,0);
3553 spvzec->SetBinContent(i,j,0);
3554 spvzec->SetBinError(i,j,0);
3555 double a = uQa->GetBinContent(i,j);
3556 double b = uQc->GetBinContent(i,j);
3557 double c = qaqc->GetBinContent(i,j);
3558 double at = qaqt->GetBinContent(i,j);
3559 double bt = qcqt->GetBinContent(i,j);
3560 if(TMath::AreEqualAbs(a,0,1e-100)) continue;
3561 if(TMath::AreEqualAbs(b,0,1e-100)) continue;
3562 if(TMath::AreEqualAbs(c,0,1e-100)) continue;
3563 if(TMath::AreEqualAbs(at,0,1e-100)) continue;
3564 if(TMath::AreEqualAbs(bt,0,1e-100)) continue;
3565 // nominal spvzea
3566 double aa = c*at/bt;
3567 if(aa<1e-100) continue;
3568 double vna = a/TMath::Sqrt(aa);
3569 spvzea->SetBinContent(i,vna);
3570 // nominal spvzec
3571 double bb = c*bt/at;
3572 if(bb<1e-100) continue;
3573 double vnc = b/TMath::Sqrt(bb);
3574 spvzec->SetBinContent(i,j,vnc);
3575 //nominal spvzewa
3576 double vnwa = (cvzea2*vna + cvzec2*vnc) / (cvzea2+cvzec2);
3577 spvzewa->SetBinContent(i,vnwa);
3578 // nominal spvzega
3579 double vnga = a*b/c;
3580 if(vnga<1e-100) continue;
3581 vnga = TMath::Sqrt(vnga);
3582 spvzega->SetBinContent(i,j,vnga);
3583 // errors
3584 int k = spvzea->GetBin(i,j);
3585 double asw = uQa->GetBinEntries(k);
3586 double bsw = uQc->GetBinEntries(k);
3587 double csw = qaqc->GetBinEntries(k);
3588 if(asw<1e-100||bsw<1e-100||csw<1e-100) continue;
3589 double asww = pasww->At(k);
3590 double bsww = pbsww->At(k);
3591 double csww = pcsww->At(k);
3592 if(asww<1e-100||bsww<1e-100||csww<1e-100) continue;
3593 if((1<asww/asw/asw)||(1<bsww/bsw/bsw)||(1<csww/csw/csw)) continue;
3594 double ab = uQauQc->GetBinContent(i,j);
3595 double ac = uQaqaqc->GetBinContent(i,j);
3596 double bc = uQcqaqc->GetBinContent(i,j);
3597 double absw = uQauQc->GetBinEntries(k);
3598 double acsw = uQaqaqc->GetBinEntries(k);
3599 double bcsw = uQcqaqc->GetBinEntries(k);
3600 if(TMath::AreEqualAbs(1,absw/asw/bsw,1e-100)||TMath::AreEqualAbs(1,bcsw/bsw/csw,1e-100)||TMath::AreEqualAbs(1,acsw/asw/csw,1e-100)) continue;
3601 double ea = uQa->GetBinError(i,j)*TMath::Sqrt(asww)/asw/TMath::Sqrt(1-asww/asw/asw);
3602 double eb = uQc->GetBinError(i,j)*TMath::Sqrt(bsww)/bsw/TMath::Sqrt(1-bsww/bsw/bsw);
3603 double ec = qaqc->GetBinError(i,j)*TMath::Sqrt(csww)/csw/TMath::Sqrt(1-csww/csw/csw);
3604 double eab = (ab-a*b)/(1-absw/asw/bsw)*absw/asw/bsw;
3605 double ebc = (bc-b*c)/(1-bcsw/bsw/csw)*bcsw/bsw/csw;
3606 double eac = (ac-a*c)/(1-acsw/asw/csw)*acsw/asw/csw;
3607 double nc, nec, neac, nebc;
3608 nc = c*at/bt;
3609 nec = ec*at/bt;
3610 neac = eac*at/bt;
3611 nebc = ebc*at/bt;
3612 double evna = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*ea*ea + a*a/4.0*nec*nec - a*TMath::Abs(nc)*neac*neac );
3613 nc = c*bt/at;
3614 nec = ec*bt/at;
3615 neac = eac*bt/at;
3616 nebc = ebc*bt/at;
3617 double evnc = 1.0/TMath::Abs(nc*nc*nc) * ( nc*nc*eb*eb + b*b/4.0*nec*nec - b*TMath::Abs(nc)*nebc*nebc );
3618 if(evna>1e-100) evna = TMath::Sqrt(evna); else evna=0;
3619 if(evnc>1e-100) evnc = TMath::Sqrt(evnc); else evnc=0;
3620 spvzea->SetBinError(i,j,evna);
3621 spvzec->SetBinError(i,j,evnc);
3622 double evnwa = TMath::Sqrt( cvzea2*cvzea2*evna*evna + cvzec2*cvzec2*evnc*evnc )/(cvzea2+cvzec2);
3623 spvzewa->SetBinError(i,j,evnwa);
3624 double evnga = 0.25/c/c * ( TMath::Abs(b*c/a)*ea*ea + TMath::Abs(a*c/b)*eb*eb + TMath::Abs(a*b/c)*ec*ec + 2*c*eab - 2*a*ebc - 2*b*eac );
3625 if(evnga>1e-100) evnga = TMath::Sqrt(evnga); else evnga=0;
3626 spvzega->SetBinError(i,j,evnga);
3627 }
3628 }
3629 printf("<<%s>> QC TPC\n",name.Data());
3630 //Qcumulants
3631 TH2D *resC2 = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum2" );
3632 TH2D *resC4 = (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_Cum4" );
3633 TH2D *resDC2= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum2" );
3634 TH2D *resDC4= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DCum4" );
3635 TH2D *resvn2= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn2" );
3636 TH2D *resvn4= (TH2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_vn4" );
3637 //correlators
3638 TProfile2D *c2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2" ));
3639 TProfile2D *c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4" ));
3640 TProfile2D *dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2" ));
3641 TProfile2D *dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC4" ));
3642 TProfile2D *c2c4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2C4" ));
3643 TProfile2D *c2dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC2" ));
3644 TProfile2D *c2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C2DC4" ));
3645 TProfile2D *c4dc2 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC2" ));
3646 TProfile2D *c4dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_C4DC4" ));
3647 TProfile2D *dc2dc4 = ((TProfile2D*)((TList*)fList->FindObject(name.Data()))->FindObject( "QC_DC2DC4" ));
3648 TArrayD *c2sww = c2->GetBinSumw2();
3649 TArrayD *c4sww = c4->GetBinSumw2();
3650 TArrayD *dc2sww= dc2->GetBinSumw2();
3651 TArrayD *dc4sww= dc4->GetBinSumw2();
3652 for(Int_t i=1; i!=resvn2->GetNbinsX()+1; ++i) {
3653 for(Int_t j=1; j!=resvn2->GetNbinsY()+1; ++j) {
3654 // cn{2}
3655 int k = c2->GetBin(i,j);
3656 double v_c2sw = c2->GetBinEntries(k);
3657 if(v_c2sw<1e-100) continue;
3658 double v_c2sww = c2sww->At(k);
3659 double v_c2 = c2->GetBinContent(i,j);
3660 double e_c2 = TMath::Sqrt(v_c2sww)/v_c2sw*c2->GetBinError(i,j);
3661 double cum2 = v_c2;
3662 double ecum2= e_c2;
3663 resC2->SetBinContent(i,j, cum2 );
3664 resC2->SetBinError(i,j, ecum2 );
3665 // cn{4}
3666 double v_c4sw = c4->GetBinEntries(k);
3667 if(v_c4sw<1e-100) continue;
3668 double v_c4sww = c4sww->At(k);
3669 double v_c4 = c4->GetBinContent(i,j);
3670 double e_c4 = TMath::Sqrt(v_c4sww)/v_c4sw*c4->GetBinError(i,j);
3671 double v_c2c4 = c2c4->GetBinContent(i,j);
3672 double v_c2c4sw = c2c4->GetBinEntries(k);
3673 if(TMath::AreEqualAbs(v_c2c4sw/v_c2sw/v_c4sw,1,1e-100)) continue;
3674 double covc2c4 = v_c2c4sw/v_c2sw/v_c4sw*(v_c2c4 - v_c2*v_c4)/(1-v_c2c4sw/v_c2sw/v_c4sw);
3675 double cum4 = v_c4 - 2*v_c2*v_c2;
3676 double ecum4= 16.0*v_c2*v_c2*e_c2*e_c2 + e_c4*e_c4 - 8.0*v_c2*covc2c4;
3677 if(ecum4<1e-100) continue;
3678 ecum4 = TMath::Sqrt( ecum4 );
3679 resC4->SetBinContent(i,j, cum4 );
3680 resC4->SetBinError(i,j, ecum4 );
3681 // dn{2}
3682 double v_dc2sw = dc2->GetBinEntries(k);
3683 if(v_dc2sw<1) continue;
3684 double v_dc2 = dc2->GetBinContent(i,j);
3685 double v_dc2sww = dc2sww->At(k);
3686 double e_dc2 = TMath::Sqrt(v_dc2sww)/v_dc2sw*dc2->GetBinError(i,j);
3687 double dcum2 = v_dc2;
3688 double edcum2= e_dc2;
3689 resDC2->SetBinContent(i,j, dcum2 );
3690 resDC2->SetBinError(i,j, edcum2 );
3691 // v2{2}
3692 if(v_c2<1e-100) continue;
3693 double dv22 = v_dc2/TMath::Sqrt(v_c2);
3694 double v_c2dc2 = c2dc2->GetBinContent(i,j);
3695 double v_c2dc2sw = c2dc2->GetBinEntries(k);
3696 if(TMath::AreEqualAbs(v_c2dc2sw/v_c2sw/v_dc2sw,1,1e-100)) continue;
3697 double covc2dc2 = v_c2dc2sw/v_c2sw/v_dc2sw*(v_c2dc2 - v_c2*v_dc2)/(1-v_c2dc2sw/v_c2sw/v_dc2sw);
3698 double edv22 = 0.25/v_c2/v_c2/v_c2*(v_dc2*v_dc2*e_c2*e_c2 + 4*v_c2*v_c2*e_dc2*e_dc2 - 4*v_c2*v_dc2*covc2dc2);
3699 //printf("%d >> dv22 %.16f || edv22^2 %.16f |||| v_c2dc2 %.16f | v_c2dc2sw %.16f | covc2dc2 %.16f \n", i,dv22,edv22,v_c2dc2,v_c2dc2sw,covc2dc2);
3700 if(edv22<1e-100) continue;
3701 edv22 = TMath::Sqrt(edv22);
3702 resvn2->SetBinContent(i,j,dv22);
3703 resvn2->SetBinError(i,j,edv22);
3704 // dn{4}
3705 double v_dc4sw = dc4->GetBinEntries(k);
3706 if(v_dc4sw<1) continue;
3707 double v_dc4 = dc4->GetBinContent(i,j);
3708 double v_dc4sww = dc4sww->At(k);
3709 double e_dc4 = TMath::Sqrt(v_dc4sww)/v_dc4sw*dc4->GetBinError(i,j);
3710 double dcum4 = v_dc4 - 2*v_c2*v_dc2;
3711 double v_c2dc4 = c2dc4->GetBinContent(i,j);
3712 double v_c2dc4sw = c2dc4->GetBinEntries(k);
3713 if(TMath::AreEqualAbs(v_c2dc4sw/v_c2sw/v_dc4sw,1,1e-100)) continue;
3714 double covc2dc4 = v_c2dc4sw/v_c2sw/v_dc4sw*(v_c2dc4 - v_c2*v_dc4)/(1-v_c2dc4sw/v_c2sw/v_dc4sw);
3715 double v_dc2dc4 = dc2dc4->GetBinContent(i,j);
3716 double v_dc2dc4sw = dc2dc4->GetBinEntries(k);
3717 if(TMath::AreEqualAbs(v_dc2dc4sw/v_dc2sw/v_dc4sw,1,1e-100)) continue;
3718 double covdc2dc4 = v_dc2dc4sw/v_dc2sw/v_dc4sw*(v_dc2dc4 - v_dc2*v_dc4)/(1-v_dc2dc4sw/v_dc2sw/v_dc4sw);
3719 double edcum4= ( +4.0*v_dc2*v_dc2*e_c2*e_c2
3720 +4.0*v_c2*v_c2*e_dc2*e_dc2
3721 +e_dc4*e_dc4
3722 +8.0*v_c2*v_dc2*covc2dc2
3723 -4.0*v_dc2*covc2dc4
3724 -4.0*v_c2*covdc2dc4 );
3725 if(edcum4<1e-100) continue;
3726 edcum4 = TMath::Sqrt(edcum4);
3727 resDC4->SetBinContent(i,j, dcum4 );
3728 resDC4->SetBinError(i,j, edcum4 );
3729 // v2{4}
3730 if(cum4>1e-100) continue;
3731 double dv24 = -dcum4/TMath::Power(-cum4,0.75);
3732 double dterm1 = 2*v_c2*v_c2*v_dc2 - 3*v_c2*v_dc4 + 2*v_c4*v_dc2;
3733 double dterm2 = 9.0/16.0*dcum4*dcum4;
3734 double dterm3 = 4.0*v_c2*v_c2*cum4*cum4;
3735 double dterm4 = cum4*cum4;
3736 double dterm5 = -3.0/2.0*dcum4*dterm1;
3737 double dterm6 = -4.0*v_c2*cum4*dterm1;
3738 double dterm7 = -2.0*cum4*dterm1;
3739 double dterm8 = 3.0*v_c2*cum4*dcum4;
3740 double dterm9 = 3.0/2.0*cum4*dcum4;
3741 double dterm10= 4*v_c2*cum4*cum4;
3742 double v_c4dc2 = c4dc2->GetBinContent(i,j);
3743 double v_c4dc2sw = c4dc2->GetBinEntries(k);
3744 if(TMath::AreEqualAbs(v_c4dc2sw/v_c4sw/v_dc2sw,1,1e-100)) continue;
3745 double covc4dc2 = v_c4dc2sw/v_c4sw/v_dc2sw*(v_c4dc2 - v_c4*v_dc2)/(1-v_c4dc2sw/v_c4sw/v_dc2sw);
3746 double v_c4dc4 = c4dc4->GetBinContent(i,j);
3747 double v_c4dc4sw = c4dc4->GetBinEntries(k);
3748 if(TMath::AreEqualAbs(v_c4dc4sw/v_c4sw/v_dc4sw,1,1e-100)) continue;
3749 double covc4dc4 = v_c4dc4sw/v_c4sw/v_dc4sw*(v_c4dc4 - v_c4*v_dc4)/(1-v_c4dc4sw/v_c4sw/v_dc4sw);
3750 double edv24= 1.0/TMath::Power(-cum4,3.5)*(+dterm1*dterm1*e_c2*e_c2
3751 +dterm2*e_c4*e_c4
3752 +dterm3*e_dc2*e_dc2
3753 +dterm4*e_dc4*e_dc4
3754 -dterm5*covc2c4
3755 -dterm6*covc2dc2
3756 +dterm7*covc2dc4
3757 +dterm8*covc4dc2
3758 -dterm9*covc4dc4
3759 -dterm10*covdc2dc4);
3760 if(edv24<1e-100) continue;
3761 edv24 = TMath::Sqrt(edv24);
3762 resvn4->SetBinContent(i,j,dv24);
3763 resvn4->SetBinError(i,j,edv24);
3764 }
3765 }
24373b38 3766}