]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fixed coverity errors
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalQCSP.cxx
CommitLineData
914042c2 1/**************************************************************************
96387d94 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
914042c2 14 **************************************************************************/
15
16
17////////////////////////////////////////////////////////////////////////
18// //
9d63a2f6 19// Task for Heavy Flavour Electron Flow with TPC plus EMCal //
914042c2 20// Non-Photonic Electron identified with Invariant mass //
21// analysis methos in function SelectPhotonicElectron //
96387d94 22// //
914042c2 23// //
24// Author: Andrea Dubla (Utrecht University) //
96387d94 25// //
914042c2 26// //
27////////////////////////////////////////////////////////////////////////
28
29#include "TChain.h"
30#include "TTree.h"
31#include "TH2F.h"
32#include "TMath.h"
33#include "TCanvas.h"
34#include "THnSparse.h"
35#include "TLorentzVector.h"
36#include "TString.h"
37#include "TFile.h"
38#include "AliAnalysisTask.h"
39#include "AliAnalysisManager.h"
40#include "AliESDEvent.h"
41#include "AliESDHandler.h"
42#include "AliAODEvent.h"
43#include "AliAODHandler.h"
a70c9e97 44#include "AliAnalysisTaskFlowTPCEMCalQCSP.h"
914042c2 45#include "TGeoGlobalMagField.h"
46#include "AliLog.h"
47#include "AliAnalysisTaskSE.h"
48#include "TRefArray.h"
49#include "TVector.h"
50#include "AliESDInputHandler.h"
51#include "AliESDpid.h"
52#include "AliAODInputHandler.h"
53#include "AliAODPid.h"
54#include "AliESDtrackCuts.h"
55#include "AliPhysicsSelection.h"
56#include "AliCentralitySelectionTask.h"
57#include "AliESDCaloCluster.h"
58#include "AliAODCaloCluster.h"
59#include "AliESDCaloTrigger.h"
60#include "AliEMCALRecoUtils.h"
61#include "AliEMCALGeometry.h"
62#include "AliGeomManager.h"
63#include "stdio.h"
64#include "TGeoManager.h"
65#include "iostream"
66#include "fstream"
67#include "AliEMCALTrack.h"
68//#include "AliEMCALTracker.h"
69#include "AliMagF.h"
70#include "AliKFParticle.h"
71#include "AliKFVertex.h"
72#include "AliPID.h"
73#include "AliPIDResponse.h"
74#include "AliHFEcontainer.h"
75#include "AliHFEcuts.h"
76#include "AliHFEpid.h"
77#include "AliHFEpidBase.h"
78#include "AliHFEpidQAmanager.h"
79#include "AliHFEtools.h"
80#include "AliCFContainer.h"
81#include "AliCFManager.h"
82#include "AliKFParticle.h"
83#include "AliKFVertex.h"
84#include "AliCentrality.h"
85#include "AliVEvent.h"
86#include "AliStack.h"
87#include "AliMCEvent.h"
88#include "TProfile.h"
89#include "AliFlowCandidateTrack.h"
90#include "AliFlowTrackCuts.h"
91#include "AliFlowEventSimple.h"
92#include "AliFlowCommonConstants.h"
93#include "AliFlowEvent.h"
94#include "TVector3.h"
95#include "TRandom2.h"
96#include "AliESDVZERO.h"
97#include "AliAODVZERO.h"
98#include "AliPID.h"
99#include "AliPIDResponse.h"
100#include "AliFlowTrack.h"
101#include "AliAnalysisTaskVnV0.h"
20791315 102#include "AliSelectNonHFE.h"
914042c2 103
9d63a2f6 104
914042c2 105class AliFlowTrackCuts;
106
9d63a2f6 107using namespace std;
914042c2 108
a70c9e97 109ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
914042c2 110//________________________________________________________________________
96387d94 111AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
112: AliAnalysisTaskSE(name)
9d63a2f6 113,fDebug(0)
114,fAOD(0)
115,fGeom(0)
116,fOutputList(0)
117,fCuts(0)
118,fIdentifiedAsOutInz(kFALSE)
119,fPassTheEventCut(kFALSE)
120,fCFM(0)
121,fPID(0)
122,fPIDqa(0)
123,fCutsRP(0) // track cuts for reference particles
124,fNullCuts(0) // dummy cuts for flow event tracks
125,fFlowEvent(0) //! flow events (one for each inv mass band)
126,fkCentralityMethod(0)
127,fCentrality(0)
128,fCentralityMin(0)
129,fCentralityMax(0)
130,fInvmassCut(0)
040cc6fa 131,fpTCut(0)
9d63a2f6 132,fTrigger(0)
133,fPhi(0)
134,fEta(0)
135,fVZEROA(0)
136,fVZEROC(0)
137,fTPCM(0)
138,fNoEvents(0)
139,fTrkEovPBef(0)
140//,fdEdxBef(0)
141,fInclusiveElecPt(0)
142,fTPCnsigma(0)
143,fTPCnsigmaAft(0)
144,fCentralityPass(0)
145,fCentralityNoPass(0)
146,fInvmassLS1(0)
147,fInvmassULS1(0)
148,fPhotoElecPt(0)
149,fSemiInclElecPt(0)
150,fULSElecPt(0)
151,fLSElecPt(0)
152,fminTPC(-1)
153,fmaxTPC(3)
154,fminEovP(0.8)
155,fmaxEovP(1.2)
156,fminM20(0.03)
157,fmaxM20(0.3)
158,fminM02(0.03)
159,fmaxM02(0.5)
160,fDispersion(1)
161,fMultCorAfterCuts(0)
162,fMultvsCentr(0)
163,fSubEventDPhiv2(0)
164,EPVzA(0)
165,EPVzC(0)
166,EPTPC(0)
167,fV2Phi(0)
168,fSparseElectronHadron(0)
169,fvertex(0)
170,fMultCorBeforeCuts(0)
171,fSideBandsFlow(kFALSE)
172,fPhiminusPsi(kFALSE)
173,fFlowEventCont(0) //! flow events (one for each inv mass band)
174,fpurity(kFALSE)
175,fSparseElectronpurity(0)
20791315 176,fOpeningAngleLS(0)
177,fOpeningAngleULS(0)
178,fNonHFE(new AliSelectNonHFE)
179,fDCA(0)
180,fOpeningAngleCut(0)
181,fOP_angle(0)
8e8d9d2a 182,fAssoTPCCluster(0)
183,fAssoITSRefit(0)
ff36a82c 184,fMultCut(0)
185,fMultCorAfterCentrBeforeCuts(0)
186,fMultCorAfterVZTRKComp(0)
187,fCentralityBeforePileup(0)
188,fCentralityAfterVZTRK(0)
9f508565 189,EPVz(0)
190,EPTPCp(0)
191,EPTPCn(0)
192,fSubEventDPhiv2new(0)
193,fV2Phivzerotot(0)
194,fHistCentrDistr(0x0)
195,fCentralityNoPassForFlattening(0)
196,fInvmassLS1highpt(0)
197,fInvmassULS1highpt(0)
914042c2 198{
96387d94 199 //Named constructor
200
201 fPID = new AliHFEpid("hfePid");
202 // Define input and output slots here
203 // Input slot #0 works with a TChain
204 DefineInput(0, TChain::Class());
205 // Output slot #0 id reserved by the base class for AOD
206 // Output slot #1 writes into a TH1 container
207 // DefineOutput(1, TH1I::Class());
208 DefineOutput(1, TList::Class());
209 DefineOutput(2, AliFlowEventSimple::Class());
210 if(fSideBandsFlow){
211 DefineOutput(3, AliFlowEventSimple::Class());
212 }
213 // DefineOutput(3, TTree::Class());
914042c2 214}
215
216//________________________________________________________________________
96387d94 217AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
218: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
9d63a2f6 219,fDebug(0)
220,fAOD(0)
221,fGeom(0)
222,fOutputList(0)
223,fCuts(0)
224,fIdentifiedAsOutInz(kFALSE)
225,fPassTheEventCut(kFALSE)
226,fCFM(0)
227,fPID(0)
228,fPIDqa(0)
229,fCutsRP(0) // track cuts for reference particles
230,fNullCuts(0) // dummy cuts for flow event tracks
231,fFlowEvent(0) //! flow events (one for each inv mass band)
232,fkCentralityMethod(0)
233,fCentrality(0)
234,fCentralityMin(0)
235,fCentralityMax(0)
236,fInvmassCut(0)
040cc6fa 237,fpTCut(0)
9d63a2f6 238,fTrigger(0)
239,fPhi(0)
240,fEta(0)
241,fVZEROA(0)
242,fVZEROC(0)
243,fTPCM(0)
244,fNoEvents(0)
245,fTrkEovPBef(0)
246//,fdEdxBef(0)
247,fInclusiveElecPt(0)
248,fTPCnsigma(0)
249,fTPCnsigmaAft(0)
250,fCentralityPass(0)
251,fCentralityNoPass(0)
252,fInvmassLS1(0)
253,fInvmassULS1(0)
254,fPhotoElecPt(0)
255,fSemiInclElecPt(0)
256,fULSElecPt(0)
257,fLSElecPt(0)
258,fminTPC(-1)
259,fmaxTPC(3)
260,fminEovP(0.8)
261,fmaxEovP(1.2)
262,fminM20(0.03)
263,fmaxM20(0.3)
264,fminM02(0.03)
265,fmaxM02(0.5)
266,fDispersion(1)
267,fMultCorAfterCuts(0)
268,fMultvsCentr(0)
269,fSubEventDPhiv2(0)
270,EPVzA(0)
271,EPVzC(0)
272,EPTPC(0)
273,fV2Phi(0)
274,fSparseElectronHadron(0)
275,fvertex(0)
276,fMultCorBeforeCuts(0)
277,fSideBandsFlow(kFALSE)
278,fPhiminusPsi(kFALSE)
279,fFlowEventCont(0) //! flow events (one for each inv mass band)
280,fpurity(kFALSE)
281,fSparseElectronpurity(0)
20791315 282,fOpeningAngleLS(0)
283,fOpeningAngleULS(0)
284,fNonHFE(new AliSelectNonHFE)
285,fDCA(0)
286,fOpeningAngleCut(0)
287,fOP_angle(0)
8e8d9d2a 288,fAssoTPCCluster(0)
289,fAssoITSRefit(0)
ff36a82c 290,fMultCut(0)
291,fMultCorAfterCentrBeforeCuts(0)
292,fMultCorAfterVZTRKComp(0)
293,fCentralityBeforePileup(0)
294,fCentralityAfterVZTRK(0)
9f508565 295,EPVz(0)
296,EPTPCp(0)
297,EPTPCn(0)
298,fSubEventDPhiv2new(0)
299,fV2Phivzerotot(0)
300,fHistCentrDistr(0x0)
301,fCentralityNoPassForFlattening(0)
302,fInvmassLS1highpt(0)
303,fInvmassULS1highpt(0)
914042c2 304{
96387d94 305 //Default constructor
306 fPID = new AliHFEpid("hfePid");
307 // Constructor
308 // Define input and output slots here
309 // Input slot #0 works with a TChain
310 DefineInput(0, TChain::Class());
311 // Output slot #0 id reserved by the base class for AOD
312 // Output slot #1 writes into a TH1 container
313 // DefineOutput(1, TH1I::Class());
314 DefineOutput(1, TList::Class());
315 DefineOutput(2, AliFlowEventSimple::Class());
9d63a2f6 316 // DefineOutput(3, TTree::Class());
96387d94 317 if(fSideBandsFlow){
318 DefineOutput(3, AliFlowEventSimple::Class());
319 }
9d63a2f6 320 //DefineOutput(3, TTree::Class());
914042c2 321}
322//_________________________________________
323
a70c9e97 324AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
914042c2 325{
96387d94 326 //Destructor
327
328 delete fOutputList;
329 delete fGeom;
330 delete fPID;
331 delete fCFM;
332 delete fPIDqa;
333 if (fDCA) delete fNonHFE;
334 if (fOutputList) delete fOutputList;
335 if (fFlowEvent) delete fFlowEvent;
336 if (fFlowEventCont) delete fFlowEventCont;
337
914042c2 338}
339//_________________________________________
340
a70c9e97 341void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
914042c2 342{
96387d94 343 //Main loop
344 //Called for each event
345
346 // create pointer to event
347
348 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
349 if (!fAOD)
350 {
351 printf("ERROR: fAOD not available\n");
352 return;
353 }
354
355 if(!fCuts)
356 {
357 AliError("HFE cuts not available");
358 return;
359 }
360
361 if(!fPID->IsInitialized())
362 {
363 // Initialize PID with the given run number
364 AliWarning("PID not initialised, get from Run no");
365 fPID->InitializePID(fAOD->GetRunNumber());
366 }
914042c2 367
9f508565 368 // cout << "kTrigger == " << fTrigger <<endl;
914042c2 369
370 if(fTrigger==0){
96387d94 371 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
914042c2 372 }
373 if(fTrigger==1){
96387d94 374 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
914042c2 375 }
376 if(fTrigger==2){
96387d94 377 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
914042c2 378 }
379 if(fTrigger==3){
96387d94 380 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
914042c2 381 }
9d63a2f6 382
96387d94 383
384 //---------------CENTRALITY AND EVENT SELECTION-----------------------
385 Int_t fNOtrks = fAOD->GetNumberOfTracks();
9d63a2f6 386 Float_t vtxz = -999;
387 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
388 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
389 TString vtxTtl = trkVtx->GetTitle();
390 if (!vtxTtl.Contains("VertexerTracks"))return;
391 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
392 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
393 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
394 vtxz = trkVtx->GetZ();
395 if(TMath::Abs(vtxz)>10)return;
96387d94 396
397 // Event cut
9d63a2f6 398 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
399 if(fNOtrks<2) return;
400
ff36a82c 401 Bool_t pass = kFALSE; //to select centrality and pile up protection
9d63a2f6 402 CheckCentrality(fAOD,pass);
403 if(!pass)return;
ff36a82c 404 fvertex->Fill(vtxz);
a70c9e97 405
96387d94 406 fNoEvents->Fill(0);
407 PlotVZeroMultiplcities(fAOD);
408
409 SetNullCuts(fAOD);
410 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
411
9d63a2f6 412 if(fSideBandsFlow){
96387d94 413 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont); //Calculate event plane Qvector and EP resolution for inclusive
9d63a2f6 414 }
415
96387d94 416 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
417 if(!pidResponse)
418 {
419 AliDebug(1, "Using default PID Response");
420 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
914042c2 421 }
914042c2 422
96387d94 423 fPID->SetPIDResponse(pidResponse);
424 fCFM->SetRecEventInfo(fAOD);
425
426 // Look for kink mother
427 Int_t numberofvertices = fAOD->GetNumberOfVertices();
428 Double_t listofmotherkink[numberofvertices];
429 Int_t numberofmotherkink = 0;
430 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
431 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
432 if(!aodvertex) continue;
433 if(aodvertex->GetType()==AliAODVertex::kKink) {
434 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
435 if(!mother) continue;
436 Int_t idmother = mother->GetID();
437 listofmotherkink[numberofmotherkink] = idmother;
438 //printf("ID %d\n",idmother);
439 numberofmotherkink++;
440 }
441 }
442
9f508565 443
96387d94 444 //=============================================V0EP from Alex======================================================================
9d63a2f6 445 Double_t qxEPa = 0, qyEPa = 0;
446 Double_t qxEPc = 0, qyEPc = 0;
9f508565 447 Double_t qxEP = 0, qyEP = 0;
9d63a2f6 448
449 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
450 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
9f508565 451 Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
9d63a2f6 452
453
454 Double_t Qx2 = 0, Qy2 = 0;
9f508565 455 Double_t Qx2p = 0, Qy2p = 0;
456 Double_t Qx2n = 0, Qy2n = 0;
9d63a2f6 457
458 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
459
460 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
461
462 if (!aodTrack)
463 continue;
464
465 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
466 continue;
467
468 if (!aodTrack->TestFilterBit(128))
469 continue;
470
9f508565 471
472 if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
473
474 Qx2p += TMath::Cos(2*aodTrack->Phi());
475 Qy2p += TMath::Sin(2*aodTrack->Phi());
476 }
477 if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
478
479 Qx2n += TMath::Cos(2*aodTrack->Phi());
480 Qy2n += TMath::Sin(2*aodTrack->Phi());
481 }
482
483
9d63a2f6 484 Qx2 += TMath::Cos(2*aodTrack->Phi());
485 Qy2 += TMath::Sin(2*aodTrack->Phi());
9f508565 486
487
488
489
9d63a2f6 490 }
491
492 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
9f508565 493 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
494 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
9d63a2f6 495
496 EPVzA->Fill(evPlAngV0A);
497 EPVzC->Fill(evPlAngV0C);
498 EPTPC->Fill(evPlAngTPC);
499
9f508565 500 EPTPCn->Fill(evPlAngTPCn);
501 EPTPCp->Fill(evPlAngTPCp);
502 EPVz->Fill(evPlAngV0);
503
504
505
506
9d63a2f6 507 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
508 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
509 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
9f508565 510
511
512 fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
513 fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
514 fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
515
516
517
518
96387d94 519 //====================================================================================================================
9d63a2f6 520
521
96387d94 522 AliAODTrack *track = NULL;
523
524 // Track loop
525 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
526 {
527 track = fAOD->GetTrack(iTracks);
528 if (!track)
914042c2 529 {
96387d94 530 printf("ERROR: Could not receive track %d\n", iTracks);
531 continue;
914042c2 532 }
96387d94 533
534 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
535 //----------hfe begin---------
536 if(track->Eta()<-0.7 || track->Eta()>0.7) continue; //eta cuts on candidates
537
538 // RecKine: ITSTPC cuts
539 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
540
541 // Reject kink mother
542 Bool_t kinkmotherpass = kTRUE;
543 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
544 if(track->GetID() == listofmotherkink[kinkmother]) {
545 kinkmotherpass = kFALSE;
546 continue;
547 }
548 }
549 if(!kinkmotherpass) continue;
550
551 // RecPrim
552 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
553 // HFEcuts: ITS layers cuts
554 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
555 // HFE cuts: TPC PID cleanup
556 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
557
558 Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
559 // Track extrapolation
560 Int_t fClsId = track->GetEMCALcluster();
561 if(fClsId < 0) continue;
562 AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
563 if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
564
565 pt = track->Pt(); //pt track after cuts
566 if(pt<fpTCut) continue;
567 fClsE = cluster->E();
568 p = track->P();
569 // dEdx = track->GetTPCsignal();
570 fEovP = fClsE/p;
571 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
572 Double_t m20 =cluster->GetM20();
573 Double_t m02 =cluster->GetM02();
574 Double_t disp=cluster->GetDispersion();
575 if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
576 fTPCnsigma->Fill(p,fTPCnSigma);
577 // fdEdxBef->Fill(p,dEdx);
578 Double_t eta = track->Eta();
579 Double_t phi = track->Phi();
580 //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
581 //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
582 Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
583 if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
584 Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
585 if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
586
587 if(fPhiminusPsi){
588 Double_t valueElh[8] = {
589 pt,
590 fEovP,
591 fTPCnSigma,
592 m20,
593 m02,
594 disp,
595 dPhi_aeh,
596 dPhi_ceh};
597 fSparseElectronHadron->Fill(valueElh);
598 }
599 //----------------------------------------------------------------------------------------------------------
600 //---------------------------From here usual electron selection---------------------------------------------
601 //----------------------------------------------------------------------------------------------------------
602 if(m20 < fminM20 || m20 > fmaxM20) continue;
603 if(m02 < fminM02 || m02 > fmaxM02) continue;
604 if(disp > fDispersion ) continue;
605 //---------------------------------for purity---------------------------------------------------------------
606 if(fpurity){
607 Double_t valuepurity[3] = {
608 pt,
609 fEovP,
610 fTPCnSigma};
9f508565 611 fSparseElectronpurity->Fill(valuepurity);
96387d94 612 }
613 //----------------------------------------------------------------------------------------------------------
614 //----------------------------------------------------------------------------------------------------------
615 if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
616 //===============================Flow Event for Contamination=============================================
617 if(fSideBandsFlow){
618 if(fEovP>0 && fEovP<0.6){
619 AliFlowTrack *sTrackCont = new AliFlowTrack();
620 sTrackCont->Set(track);
621 sTrackCont->SetID(track->GetID());
622 sTrackCont->SetForRPSelection(kFALSE);
623 sTrackCont->SetForPOISelection(kTRUE);
624 sTrackCont->SetMass(2637);
625 for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
626 {
627 // cout << " no of rps " << iRPs << endl;
628 AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
629 if (!iRPCont) continue;
630 if (!iRPCont->InRPSelection()) continue;
631 if( sTrackCont->GetID() == iRPCont->GetID())
632 {
633 if(fDebug) printf(" was in RP set");
634 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
635 iRPCont->SetForRPSelection(kFALSE);
636 fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
637 }
638 } //end of for loop on RPs
639 fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
640 }
641 }
642 //==========================================================================================================
643 //===============================From here eovP cut is used fro QC, SP and EPV0=============================
644 if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
645 //==========================================================================================================
646 //============================Event Plane Method with V0====================================================
647 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
648 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
649 Double_t v2Phi[3] = {
650 v2PhiV0A,
651 v2PhiV0C,
652 pt};
653 fV2Phi->Fill(v2Phi);
9f508565 654
655 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
656 Double_t v2PhiV0tot[2] = {
657 v2PhiVz,
658 pt};
659 fV2Phivzerotot->Fill(v2PhiV0tot);
96387d94 660 //=========================================================================================================
661 fTPCnsigmaAft->Fill(p,fTPCnSigma);
662 fInclusiveElecPt->Fill(pt);
663 fPhi->Fill(phi);
664 fEta->Fill(eta);
665 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
666 AliFlowTrack *sTrack = new AliFlowTrack();
667 sTrack->Set(track);
668 sTrack->SetID(track->GetID());
669 sTrack->SetForRPSelection(kFALSE);
670 sTrack->SetForPOISelection(kTRUE);
671 sTrack->SetMass(263732);
672 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
673 {
674 // cout << " no of rps " << iRPs << endl;
675 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
676 if (!iRP) continue;
677 if (!iRP->InRPSelection()) continue;
678 if( sTrack->GetID() == iRP->GetID())
679 {
680 if(fDebug) printf(" was in RP set");
681 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
682 iRP->SetForRPSelection(kFALSE);
683 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
684 }
685 } //end of for loop on RPs
686 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
687
688
689
690 if(fDCA){
691 //----------------------Selection of Photonic Electrons DCA-----------------------------
692 fNonHFE = new AliSelectNonHFE();
693 fNonHFE->SetAODanalysis(kTRUE);
694 fNonHFE->SetInvariantMassCut(fInvmassCut);
695 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
696 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
697 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
698 fNonHFE->SetAlgorithm("DCA"); //KF
699 fNonHFE->SetPIDresponse(pidResponse);
700 fNonHFE->SetTrackCuts(-3,3);
701
702 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
703 fNonHFE->SetHistAngle(fOpeningAngleULS);
704 //fNonHFE->SetHistDCABack(fDCABack);
705 //fNonHFE->SetHistDCA(fDCA);
706 fNonHFE->SetHistMassBack(fInvmassLS1);
707 fNonHFE->SetHistMass(fInvmassULS1);
708
709 fNonHFE->FindNonHFE(iTracks,track,fAOD);
710
711 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
712 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
713 // Bool_t fUlsIsPartner = kFALSE;
714 // Bool_t fLsIsPartner = kFALSE;
715 if(fNonHFE->IsULS()){
716 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
717 fULSElecPt->Fill(track->Pt());
718 }
719 }
720
721 if(fNonHFE->IsLS()){
722 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
723 fLSElecPt->Fill(track->Pt());
724 }
725 }
726 }
727
728 if(!fDCA){
729 //----------------------Selection of Photonic Electrons KFParticle-----------------------------
730 Bool_t fFlagPhotonicElec = kFALSE;
731 SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
732 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
733 // Semi inclusive electron
734 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
735 }
736
737
738
739 }//end loop on track
740
741 PostData(1, fOutputList);
742 PostData(2, fFlowEvent);
9d63a2f6 743 if(fSideBandsFlow){
96387d94 744 PostData(3, fFlowEventCont);
9d63a2f6 745 }
96387d94 746
747 //----------hfe end---------
914042c2 748}
749//_________________________________________
a70c9e97 750void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
914042c2 751{
96387d94 752 //Identify non-heavy flavour electrons using Invariant mass method KF
753
754 Bool_t flagPhotonicElec = kFALSE;
755
756 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
757 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
758 if (!trackAsso) {
759 printf("ERROR: Could not receive track %d\n", jTracks);
760 continue;
761 }
762 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
763 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
764 // if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit) || (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
765
766 if(fAssoITSRefit){
767 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
768 }
769
770 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
771
772 if(jTracks == itrack) continue;
773 Double_t ptAsso=-999., nsigma=-999.0;
774 Double_t mass=-999., width = -999;
775 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
776 Double_t openingAngle = -999.;
9f508565 777 Double_t ptcutonmasshighpt = track->Pt();
778
96387d94 779 ptAsso = trackAsso->Pt();
780 Short_t chargeAsso = trackAsso->Charge();
781 Short_t charge = track->Charge();
9f508565 782
96387d94 783 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
784
785 //80
786 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
787 if(nsigma < -3 || nsigma > 3) continue;
788 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
789 if(ptAsso <0.3) continue;
790
791 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
792 if(charge>0) fPDGe1 = -11;
793 if(chargeAsso>0) fPDGe2 = -11;
794
795 if(charge == chargeAsso) fFlagLS = kTRUE;
796 if(charge != chargeAsso) fFlagULS = kTRUE;
797
798 AliKFParticle::SetField(fAOD->GetMagneticField());
799 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
800 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
801 AliKFParticle recg(ge1, ge2);
802
803 if(recg.GetNDF()<1) continue;
804 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
805 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
806 recg.GetMass(mass,width);
807
808 openingAngle = ge1.GetAngle(ge2);
809 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
810 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
811 if(fOP_angle)if(openingAngle > fOpeningAngleCut) continue;
812
813
814 if(fFlagLS) fInvmassLS1->Fill(mass);
815 if(fFlagULS) fInvmassULS1->Fill(mass);
816
9f508565 817 if(ptcutonmasshighpt >= 8.){
818 if(fFlagLS) fInvmassLS1highpt->Fill(mass);
819 if(fFlagULS) fInvmassULS1highpt->Fill(mass);
820 }
821
822
96387d94 823 if(mass<fInvmassCut){
824 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
825 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
826 }
827
828 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
829 flagPhotonicElec = kTRUE;
830 }
831 }//track loop
832 fFlagPhotonicElec = flagPhotonicElec;
914042c2 833}
8e8d9d2a 834//__________________________________________________________________________________
a70c9e97 835void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
914042c2 836{
96387d94 837 //Create histograms
838
839 //----------hfe initialising begin---------
840 fNullCuts = new AliFlowTrackCuts("null_cuts");
841
842 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
843 cc->SetNbinsMult(10000);
844 cc->SetMultMin(0);
845 cc->SetMultMax(10000);
846
847 cc->SetNbinsPt(100);
848 cc->SetPtMin(0);
849 cc->SetPtMax(50);
850
851 cc->SetNbinsPhi(180);
852 cc->SetPhiMin(0.0);
853 cc->SetPhiMax(TMath::TwoPi());
854
855 cc->SetNbinsEta(30);
856 cc->SetEtaMin(-7.0);
857 cc->SetEtaMax(+7.0);
858
859 cc->SetNbinsQ(500);
860 cc->SetQMin(0.0);
861 cc->SetQMax(3.0);
862
863 //--------Initialize PID
864 fPID->SetHasMCData(kFALSE);
865 if(!fPID->GetNumberOfPIDdetectors())
866 {
867 fPID->AddDetector("TPC", 0);
868 fPID->AddDetector("EMCAL", 1);
869 }
870
871 fPID->SortDetectors();
872 fPIDqa = new AliHFEpidQAmanager();
873 fPIDqa->Initialize(fPID);
874
875 //--------Initialize correction Framework and Cuts
876 fCFM = new AliCFManager;
877 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
878 fCFM->SetNStepParticle(kNcutSteps);
879 for(Int_t istep = 0; istep < kNcutSteps; istep++)
880 fCFM->SetParticleCutsList(istep, NULL);
881
882 if(!fCuts){
883 AliWarning("Cuts not available. Default cuts will be used");
884 fCuts = new AliHFEcuts;
885 fCuts->CreateStandardCuts();
886 }
887
888 fCuts->SetAOD();
889 fCuts->Initialize(fCFM);
890 //----------hfe initialising end--------
891 //---------Output Tlist
892 fOutputList = new TList();
893 fOutputList->SetOwner();
894 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
895
896 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
897 fOutputList->Add(fNoEvents);
898
899 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
900 fOutputList->Add(fTPCnsigma);
901
902 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
903 fOutputList->Add(fTPCnsigmaAft);
904
905 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
906 fOutputList->Add(fTrkEovPBef);
907
908 // fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
909 // fOutputList->Add(fdEdxBef);
910
911 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
912 fOutputList->Add(fInclusiveElecPt);
913
914 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
915 fOutputList->Add(fPhotoElecPt);
ff36a82c 916
96387d94 917 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
918 fOutputList->Add(fSemiInclElecPt);
919
920 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
921 fOutputList->Add(fULSElecPt);
922
914042c2 923 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
924 fOutputList->Add(fLSElecPt);
96387d94 925
926 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
927 fOutputList->Add(fInvmassLS1);
928
929 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
930 fOutputList->Add(fInvmassULS1);
9f508565 931
932 fInvmassLS1highpt = new TH1F("fInvmassLS1highpt", "Inv mass of LS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
933 fOutputList->Add(fInvmassLS1highpt);
934
935 fInvmassULS1highpt = new TH1F("fInvmassULS1highpt", "Inv mass of ULS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
936 fOutputList->Add(fInvmassULS1highpt);
96387d94 937
938 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
939 fOutputList->Add(fCentralityPass);
940
941 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
942 fOutputList->Add(fCentralityNoPass);
943
9f508565 944 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
945 fOutputList->Add(fCentralityNoPassForFlattening);
946
ff36a82c 947 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
948 fOutputList->Add(fCentralityBeforePileup);
9f508565 949
ff36a82c 950 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
951 fOutputList->Add(fCentralityAfterVZTRK);
9f508565 952
96387d94 953 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
954 fOutputList->Add(fPhi);
955
956 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
957 fOutputList->Add(fEta);
958
959 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
960 fOutputList->Add(fVZEROA);
961
962 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
963 fOutputList->Add(fVZEROC);
964
965 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
966 fOutputList->Add(fTPCM);
967
968 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
969 fOutputList->Add(fvertex);
970
971 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
972 fOutputList->Add(fMultCorBeforeCuts);
973
974 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
975 fOutputList->Add(fMultCorAfterCuts);
976
ff36a82c 977 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
978 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
979
980 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
981 fOutputList->Add(fMultCorAfterVZTRKComp);
982
96387d94 983 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
984 fOutputList->Add(fMultvsCentr);
9d63a2f6 985
20791315 986 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
987 fOutputList->Add(fOpeningAngleLS);
988
989 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
990 fOutputList->Add(fOpeningAngleULS);
991
96387d94 992
993 //----------------------------------------------------------------------------
9d63a2f6 994 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
995 fOutputList->Add(EPVzA);
996 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
997 fOutputList->Add(EPVzC);
998 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
999 fOutputList->Add(EPTPC);
9f508565 1000
1001
1002 EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1003 fOutputList->Add(EPVz);
1004 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1005 fOutputList->Add(EPTPCp);
1006 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1007 fOutputList->Add(EPTPCn);
1008
96387d94 1009 //----------------------------------------------------------------------------
9d63a2f6 1010 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1011 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1012 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1013 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1014 fOutputList->Add(fSubEventDPhiv2);
9f508565 1015
1016
1017
1018 fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1019 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1020 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1021 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1022 fOutputList->Add(fSubEventDPhiv2new);
1023 //================================Event Plane with VZERO A & C=====================
1024 const Int_t nPtBins = 12;
1025 Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10, 13};
9d63a2f6 1026 // v2A, v2C, pt
1027 Int_t bins[3] = { 50, 50, nPtBins};
1028 Double_t xmin[3] = { -1., -1., 0};
9f508565 1029 Double_t xmax[3] = { 1., 1., 13};
9d63a2f6 1030 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1031 // Set bin limits for axes which are not standard binned
1032 fV2Phi->SetBinEdges(2, binsPt);
1033 // set axes titles
1034 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1035 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1036 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1037 fOutputList->Add(fV2Phi);
9f508565 1038
1039
1040
1041 //================================Event Plane with VZERO=====================
1042 // const Int_t nPtBins = 10;
1043 // Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
1044 // v2, pt
1045 Int_t binsV[2] = { 50, nPtBins};
1046 Double_t xminV[2] = { -1., 0};
1047 Double_t xmaxV[2] = { 1., 13};
1048 fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1049 // Set bin limits for axes which are not standard binned
1050 fV2Phivzerotot->SetBinEdges(1, binsPt);
1051 // set axes titles
1052 fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1053 fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1054 fOutputList->Add(fV2Phivzerotot);
1055
1056
1057
96387d94 1058 //----------------------------------------------------------------------------
9f508565 1059
9d63a2f6 1060 if(fPhiminusPsi){
96387d94 1061 Int_t binsvElectH[8]={ 600, 200, 200 ,100, 100, 100, 10, 10}; //pt, E/p,TPCnSigma,M20,M02,Disp Phi-psiV0A ,Phi-PsiV0C,eta (commented)
1062 Double_t xminvElectH[8]={0, 0, -10 , 0, 0, 0, 0, 0};
1063 Double_t xmaxvElectH[8]={20, 2, 10 , 2, 2, 2, TMath::Pi(), TMath::Pi()};
1064 fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
1065 fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1066 fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
1067 fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
1068 fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
1069 fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
1070 fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
1071 fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
1072 fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
1073 fOutputList->Add(fSparseElectronHadron);
9d63a2f6 1074 }
96387d94 1075 //----------------------------------------------------------------------------
1076 //----------------------------------------------------------------------------
9d63a2f6 1077 if(fpurity){
96387d94 1078 Int_t binsvpurity[3]={ 600,200, 200}; //pt, E/p,TPCnSigma
1079 Double_t xminvpurity[3]={0, 0, -10};
ff36a82c 1080 Double_t xmaxvpurity[3]={30, 2, 10};
96387d94 1081 fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
1082 fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1083 fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
1084 fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
1085 fOutputList->Add(fSparseElectronpurity);
9d63a2f6 1086 }
96387d94 1087 //----------------------------------------------------------------------------
1088
1089 PostData(1,fOutputList);
1090 // create and post flowevent
1091 fFlowEvent = new AliFlowEvent(10000);
1092 PostData(2, fFlowEvent);
9d63a2f6 1093
1094 if(fSideBandsFlow){
96387d94 1095 fFlowEventCont = new AliFlowEvent(10000);
1096 PostData(3, fFlowEventCont);
9d63a2f6 1097 }
96387d94 1098}
914042c2 1099
1100//________________________________________________________________________
a70c9e97 1101void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
914042c2 1102{
96387d94 1103 // Info("Terminate");
1104 AliAnalysisTaskSE::Terminate();
914042c2 1105}
1106//_____________________________________________________________________________
a70c9e97 1107template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
914042c2 1108{
96387d94 1109 // QA multiplicity plots
1110 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1111 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
914042c2 1112}
1113//_____________________________________________________________________________
a70c9e97 1114template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
914042c2 1115{
96387d94 1116 //Set null cuts
914042c2 1117 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1118 fCutsRP->SetEvent(event, MCEvent());
1119 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1120 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1121 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1122 fNullCuts->SetEvent(event, MCEvent());
1123}
1124//_____________________________________________________________________________
96387d94 1125void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
914042c2 1126{
96387d94 1127 //Prepare flow events
1128 FlowEv->ClearFast();
1129 FlowEv->Fill(fCutsRP, fNullCuts);
1130 FlowEv->SetReferenceMultiplicity(iMulti);
1131 FlowEv->DefineDeadZone(0, 0, 0, 0);
1132 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
914042c2 1133}
1134//_____________________________________________________________________________
a70c9e97 1135Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
914042c2 1136{
96387d94 1137 // Check single track cuts for a given cut step
1138 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1139 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1140 return kTRUE;
914042c2 1141}
1142//_________________________________________
a70c9e97 1143void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
914042c2 1144{
9f508565 1145 //============================Multiplicity TPV vs Global===============================================================================
9d63a2f6 1146 const Int_t nGoodTracks = event->GetNumberOfTracks();
9d63a2f6 1147 Float_t multTPC(0.); // tpc mult estimate
1148 Float_t multGlob(0.); // global multiplicity
1149 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1150 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1151 if (!trackAOD) continue;
1152 if (!(trackAOD->TestFilterBit(1))) continue;
1153 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
1154 multTPC++;
1155 }
1156 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1157 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1158 if (!trackAOD) continue;
1159 if (!(trackAOD->TestFilterBit(16))) continue;
1160 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
1161 Double_t b[2] = {-99., -99.};
1162 Double_t bCov[3] = {-99., -99., -99.};
1163 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1164 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1165 multGlob++;
9f508565 1166 } //track loop
1167 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1168 //============================================================================================================================
ff36a82c 1169 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1170 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1171 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
9f508565 1172 // cout << "--------------Centrality evaluated-------------------------"<<endl;
ff36a82c 1173 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1174 {
1175 fCentralityNoPass->Fill(fCentrality);
1176 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1177 centralitypass = kFALSE;
1178 }else
1179 {
1180 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1181 centralitypass = kTRUE;
1182 }
1183 if (centralitypass){
9f508565 1184 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1185 fCentralityBeforePileup->Fill(fCentrality);
1186 }//...after centrality selectrion
1187 //============================================================================================================================
1188 //to remove the bias introduced by multeplicity outliers---------------------
ff36a82c 1189 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1190 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1191 if (TMath::Abs(centv0 - centTrk) > 5.0){
1192 centralitypass = kFALSE;
1193 fCentralityNoPass->Fill(fCentrality);
1194 }
1195 if (centralitypass){
9f508565 1196 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1197 fCentralityAfterVZTRK->Fill(fCentrality);
ff36a82c 1198 }//...after centrality selectrion
9f508565 1199 //============================================================================================================================
1200 if(fMultCut){
1201 if(fTrigger==1){
1202 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1203 // cout <<" Trigger ==" <<fTrigger<< endl;
1204 centralitypass = kFALSE;
1205 fCentralityNoPass->Fill(fCentrality);
1206 }//2011 Semicentral
1207 }
1208 if(fTrigger==0){
1209 if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1210 // cout <<" Trigger ==" <<fTrigger<< endl;
1211 centralitypass = kFALSE;
1212 fCentralityNoPass->Fill(fCentrality);
1213 }//2011
1214 }//2011 Central
1215 }
1216 //=================================All cuts are passed==================++++==================================================
1217 //=================================Now Centrality flattening for central trigger==================++++==================================================
edaa818d 1218 if(fTrigger==0){
9f508565 1219 if(!IsEventSelectedForCentrFlattening(fCentrality)){
edaa818d 1220 centralitypass = kFALSE;
9f508565 1221 fCentralityNoPassForFlattening->Fill(fCentrality);
1222 }
1223 }
1224 //==============================fill histo after all cuts==============================++++==================================================
9d63a2f6 1225 if(centralitypass){
96387d94 1226 fCentralityPass->Fill(fCentrality);
1227 fMultCorAfterCuts->Fill(multGlob, multTPC);
1228 fMultvsCentr->Fill(fCentrality, multTPC);
a70c9e97 1229 }
914042c2 1230}
1231//_____________________________________________________________________________
a70c9e97 1232void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
914042c2 1233{
96387d94 1234 // Set a centrality range ]min, max] and define the method to use for centrality selection
1235 fCentralityMin = CentralityMin;
1236 fCentralityMax = CentralityMax;
1237 fkCentralityMethod = CentralityMethod;
914042c2 1238}
1239//_____________________________________________________________________________
a70c9e97 1240void AliAnalysisTaskFlowTPCEMCalQCSP::SetIDCuts(Double_t minTPC, Double_t maxTPC, Double_t minEovP, Double_t maxEovP, Double_t minM20, Double_t maxM20, Double_t minM02, Double_t maxM02, Double_t Dispersion)
914042c2 1241{
1242 //Set ID cuts
1243 fminTPC = minTPC;
1244 fmaxTPC = maxTPC;
1245 fminEovP = minEovP;
1246 fmaxEovP = maxEovP;
1247 fminM20 = minM20;
1248 fmaxM20 = maxM20;
1249 fminM02 = minM02;
1250 fmaxM02 = maxM02;
1251 fDispersion = Dispersion;
1252}
1253//_____________________________________________________________________________
9f508565 1254//_____________________________________________________________________________
1255void AliAnalysisTaskFlowTPCEMCalQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1256 // set the histo for centrality flattening
1257 // the centrality is flatten in the range minCentr,maxCentr
1258 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1259 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1260 // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
1261 // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
1262
1263 if(maxCentr<minCentr){
1264 AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1265 }
1266
1267 if(fHistCentrDistr)delete fHistCentrDistr;
1268 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1269 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1270 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1271 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1272 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1273 Double_t ref=0.,bincont=0.,binrefwidth=1.;
1274 Int_t binref=0;
1275 if(TMath::Abs(centrRef)<0.0001){
1276 binref=fHistCentrDistr->GetMinimumBin();
1277 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1278 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1279 }
1280 else if(centrRef>0.){
1281 binref=h->FindBin(centrRef);
1282 if(binref<1||binref>h->GetNbinsX()){
1283 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1284 }
1285 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1286 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1287 }
1288 else{
1289 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1290 binref=fHistCentrDistr->GetMaximumBin();
1291 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1292 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1293 }
1294
1295 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1296 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1297 bincont=h->GetBinContent(j);
1298 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1299 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1300 }
1301 else{
1302 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1303 }
1304 }
1305
1306 fHistCentrDistr->SetBinContent(0,switchTRand);
1307 return;
1308
1309}
1310
1311//-------------------------------------------------
1312Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1313 //
1314 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1315 // Can be faster if it was required that fHistCentrDistr covers
1316 // exactly the desired centrality range (e.g. part of the lines below should be done during the
1317 // setting of the histo) and TH1::SetMinimum called
1318 //
1319
1320 if(!fHistCentrDistr) return kTRUE;
1321 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1322 // if(maxbin>fHistCentrDistr->GetNbinsX()){
1323 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1324 // }
1325
1326 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1327 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1328 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1329
1330 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1331 if(gRandom->Uniform(1.)<bincont)return kTRUE;
1332 return kFALSE;
1333 }
1334
1335 if(centDigits*100.<bincont)return kTRUE;
1336 return kFALSE;
1337
1338}
1339//---------------------------------------------------------------------------
1340
1341
1342
1343
1344