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