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