]> 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
Merge branch 'feature-movesplit'
[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"
e7195f3b 77#include "AliHFEpidTPC.h"
914042c2 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"
20791315 103#include "AliSelectNonHFE.h"
914042c2 104
9d63a2f6 105
914042c2 106class AliFlowTrackCuts;
107
9d63a2f6 108using namespace std;
914042c2 109
a70c9e97 110ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
914042c2 111//________________________________________________________________________
96387d94 112AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
113: AliAnalysisTaskSE(name)
9d63a2f6 114,fDebug(0)
115,fAOD(0)
e7195f3b 116,fVevent(0)
9d63a2f6 117,fGeom(0)
118,fOutputList(0)
119,fCuts(0)
120,fIdentifiedAsOutInz(kFALSE)
121,fPassTheEventCut(kFALSE)
122,fCFM(0)
123,fPID(0)
e7195f3b 124,ftpcpid(0)
9d63a2f6 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)
040cc6fa 134,fpTCut(0)
9d63a2f6 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)
20791315 179,fOpeningAngleLS(0)
180,fOpeningAngleULS(0)
181,fNonHFE(new AliSelectNonHFE)
182,fDCA(0)
183,fOpeningAngleCut(0)
184,fOP_angle(0)
8e8d9d2a 185,fAssoTPCCluster(0)
186,fAssoITSRefit(0)
ff36a82c 187,fMultCut(0)
188,fMultCorAfterCentrBeforeCuts(0)
189,fMultCorAfterVZTRKComp(0)
190,fCentralityBeforePileup(0)
191,fCentralityAfterVZTRK(0)
acd6f996 192,fCentralityAfterCorrCut(0)
193,fMultCorAfterCorrCut(0)
9f508565 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)
7fcd6b0f 203,fSparsephipsiULS(0)
204,fSparsephipsiLS(0)
205,fSparseMassULS(0)
206,fSparseMassLS(0)
52a564ee 207,fHistEPDistrWeight(0)
208,EPweights(0)
209,EPVzAftW(0)
e7195f3b 210,multCorrection(0)
211,fptminAsso(0)
914042c2 212{
96387d94 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());
914042c2 228}
229
230//________________________________________________________________________
96387d94 231AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
232: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
9d63a2f6 233,fDebug(0)
234,fAOD(0)
e7195f3b 235,fVevent(0)
9d63a2f6 236,fGeom(0)
237,fOutputList(0)
238,fCuts(0)
239,fIdentifiedAsOutInz(kFALSE)
240,fPassTheEventCut(kFALSE)
241,fCFM(0)
242,fPID(0)
e7195f3b 243,ftpcpid(0)
9d63a2f6 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)
040cc6fa 253,fpTCut(0)
9d63a2f6 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)
20791315 298,fOpeningAngleLS(0)
299,fOpeningAngleULS(0)
300,fNonHFE(new AliSelectNonHFE)
301,fDCA(0)
302,fOpeningAngleCut(0)
303,fOP_angle(0)
8e8d9d2a 304,fAssoTPCCluster(0)
305,fAssoITSRefit(0)
ff36a82c 306,fMultCut(0)
307,fMultCorAfterCentrBeforeCuts(0)
308,fMultCorAfterVZTRKComp(0)
309,fCentralityBeforePileup(0)
310,fCentralityAfterVZTRK(0)
acd6f996 311,fCentralityAfterCorrCut(0)
312,fMultCorAfterCorrCut(0)
9f508565 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)
7fcd6b0f 322,fSparsephipsiULS(0)
323,fSparsephipsiLS(0)
324,fSparseMassULS(0)
325,fSparseMassLS(0)
52a564ee 326,fHistEPDistrWeight(0)
327,EPweights(0)
328,EPVzAftW(0)
e7195f3b 329,multCorrection(0)
330,fptminAsso(0)
914042c2 331{
96387d94 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());
9d63a2f6 343 // DefineOutput(3, TTree::Class());
96387d94 344 if(fSideBandsFlow){
345 DefineOutput(3, AliFlowEventSimple::Class());
346 }
9d63a2f6 347 //DefineOutput(3, TTree::Class());
914042c2 348}
349//_________________________________________
350
a70c9e97 351AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
914042c2 352{
96387d94 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
914042c2 365}
366//_________________________________________
367
a70c9e97 368void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
914042c2 369{
96387d94 370 //Main loop
371 //Called for each event
372
373 // create pointer to event
374
375 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
e7195f3b 376 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
377
96387d94 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 }
914042c2 396
9f508565 397 // cout << "kTrigger == " << fTrigger <<endl;
914042c2 398 if(fTrigger==0){
868965a9 399 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral)) return;
914042c2 400 }
401 if(fTrigger==1){
868965a9 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
914042c2 409 }
410 if(fTrigger==2){
868965a9 411 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA)) return;
914042c2 412 }
413 if(fTrigger==3){
868965a9 414 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
914042c2 415 }
acd6f996 416 if(fTrigger==4){
417 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kCentral | AliVEvent::kSemiCentral))) return;
418 }
868965a9 419 if(fTrigger==5){
420 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
421 }
acd6f996 422
96387d94 423
424 //---------------CENTRALITY AND EVENT SELECTION-----------------------
425 Int_t fNOtrks = fAOD->GetNumberOfTracks();
9d63a2f6 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;
96387d94 436
437 // Event cut
9d63a2f6 438 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
439 if(fNOtrks<2) return;
440
ff36a82c 441 Bool_t pass = kFALSE; //to select centrality and pile up protection
9d63a2f6 442 CheckCentrality(fAOD,pass);
443 if(!pass)return;
ff36a82c 444 fvertex->Fill(vtxz);
a70c9e97 445
96387d94 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
9d63a2f6 452 if(fSideBandsFlow){
96387d94 453 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont); //Calculate event plane Qvector and EP resolution for inclusive
9d63a2f6 454 }
455
96387d94 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());
914042c2 461 }
914042c2 462
96387d94 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
9f508565 483
96387d94 484 //=============================================V0EP from Alex======================================================================
9d63a2f6 485 Double_t qxEPa = 0, qyEPa = 0;
486 Double_t qxEPc = 0, qyEPc = 0;
9f508565 487 Double_t qxEP = 0, qyEP = 0;
9d63a2f6 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);
9f508565 491 Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
9d63a2f6 492
493
494 Double_t Qx2 = 0, Qy2 = 0;
9f508565 495 Double_t Qx2p = 0, Qy2p = 0;
496 Double_t Qx2n = 0, Qy2n = 0;
9d63a2f6 497
498 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
499
f15c1f69 500 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iT));
501 if(!aodTrack) AliFatal("Not a standard AOD");
9d63a2f6 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
9f508565 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
9d63a2f6 525 Qx2 += TMath::Cos(2*aodTrack->Phi());
526 Qy2 += TMath::Sin(2*aodTrack->Phi());
9f508565 527
528
529
530
9d63a2f6 531 }
532
533 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
9f508565 534 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
535 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
9d63a2f6 536
537 EPVzA->Fill(evPlAngV0A);
538 EPVzC->Fill(evPlAngV0C);
539 EPTPC->Fill(evPlAngTPC);
540
9f508565 541 EPTPCn->Fill(evPlAngTPCn);
542 EPTPCp->Fill(evPlAngTPCp);
543 EPVz->Fill(evPlAngV0);
544
545
52a564ee 546 Double_t weightEP =1;
547 if(EPweights){
548 weightEP = GiveMeWeight(evPlAngV0);
549 EPVzAftW->Fill(evPlAngV0,weightEP);
550
551 }
552
9f508565 553
9d63a2f6 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
9f508565 557
558
52a564ee 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 }
96387d94 569 //====================================================================================================================
52a564ee 570
9d63a2f6 571
572
96387d94 573 AliAODTrack *track = NULL;
574
575 // Track loop
576 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
577 {
f15c1f69 578 track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
579 if(!track) AliFatal("Not a standard AOD");
96387d94 580 if (!track)
914042c2 581 {
96387d94 582 printf("ERROR: Could not receive track %d\n", iTracks);
583 continue;
914042c2 584 }
96387d94 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;
e7195f3b 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
96387d94 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};
9f508565 676 fSparseElectronpurity->Fill(valuepurity);
96387d94 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());
03b64410 687 sTrackCont->SetForRPSelection(kTRUE);
96387d94 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);
03b64410 701 // fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
96387d94 702 }
703 } //end of for loop on RPs
704 fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
03b64410 705 fFlowEventCont->SetNumberOfPOIs(fFlowEventCont->GetNumberOfPOIs()+1);
706
96387d94 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);
9f508565 721
722 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
723 Double_t v2PhiV0tot[2] = {
724 v2PhiVz,
725 pt};
52a564ee 726
727 if(EPweights) fV2Phivzerotot->Fill(v2PhiV0tot,weightEP);
728 if(!EPweights) fV2Phivzerotot->Fill(v2PhiV0tot);
729
730
731
96387d94 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());
03b64410 741 sTrack->SetForRPSelection(kTRUE);
96387d94 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);
03b64410 755 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
96387d94 756 }
757 } //end of for loop on RPs
758 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
03b64410 759 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
96387d94 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;
e7195f3b 803 SelectPhotonicElectron(iTracks,track,fEovP, evPlAngV0, fFlagPhotonicElec,weightEP,mult);
96387d94 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);
9d63a2f6 815 if(fSideBandsFlow){
96387d94 816 PostData(3, fFlowEventCont);
9d63a2f6 817 }
96387d94 818
819 //----------hfe end---------
914042c2 820}
821//_________________________________________
e7195f3b 822void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Double_t fEovP,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec, Double_t weightEPflat, Double_t multev)
914042c2 823{
96387d94 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++){
f15c1f69 829 AliAODTrack *trackAsso = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(jTracks));
830 if(!trackAsso) AliFatal("Not a standard AOD");
96387d94 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.;
9f508565 850 Double_t ptcutonmasshighpt = track->Pt();
acd6f996 851
96387d94 852 ptAsso = trackAsso->Pt();
853 Short_t chargeAsso = trackAsso->Charge();
854 Short_t charge = track->Charge();
acd6f996 855
e7195f3b 856 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
857 if(ptAsso <fptminAsso) continue;
96387d94 858 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
e7195f3b 859
860 Double_t CorrectTPCNSigma;
861 if(multCorrection){
862 CorrectTPCNSigma = ftpcpid->GetCorrectedTPCnSigma(trackAsso->Eta(), multev, nsigma);
863 nsigma = CorrectTPCNSigma;
864 }
96387d94 865
866 //80
867 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
868 if(nsigma < -3 || nsigma > 3) continue;
96387d94 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
7fcd6b0f 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
9f508565 912 if(ptcutonmasshighpt >= 8.){
acd6f996 913 if(fFlagLS) fInvmassLS1highpt->Fill(mass);
914 if(fFlagULS) fInvmassULS1highpt->Fill(mass);
9f508565 915 }
916
917
96387d94 918 if(mass<fInvmassCut){
919 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
920 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
921 }
922
7fcd6b0f 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(),
52a564ee 934 fEovP,
7fcd6b0f 935 DeltaPhi_eEP
936 };
52a564ee 937 if(EPweights) fSparsephipsiULS->Fill(ulsSparse,weightEPflat);
938 if(!EPweights) fSparsephipsiULS->Fill(ulsSparse);
7fcd6b0f 939 }
940 if(fFlagLS){
941 Double_t lsSparse[3] = {
942 track->Pt(),
52a564ee 943 fEovP,
7fcd6b0f 944 DeltaPhi_eEP
945 };
52a564ee 946 if(EPweights) fSparsephipsiLS->Fill(lsSparse,weightEPflat);
947 if(!EPweights)fSparsephipsiLS->Fill(lsSparse);
948 }
7fcd6b0f 949 }
950
951
952
96387d94 953 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
954 flagPhotonicElec = kTRUE;
955 }
956 }//track loop
957 fFlagPhotonicElec = flagPhotonicElec;
914042c2 958}
8e8d9d2a 959//__________________________________________________________________________________
a70c9e97 960void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
914042c2 961{
96387d94 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);
ff36a82c 1041
96387d94 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
914042c2 1048 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
1049 fOutputList->Add(fLSElecPt);
96387d94 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);
acd6f996 1056
9f508565 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);
96387d94 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
9f508565 1069 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1070 fOutputList->Add(fCentralityNoPassForFlattening);
1071
ff36a82c 1072 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1073 fOutputList->Add(fCentralityBeforePileup);
9f508565 1074
ff36a82c 1075 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1076 fOutputList->Add(fCentralityAfterVZTRK);
acd6f996 1077
1078 fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1079 fOutputList->Add(fCentralityAfterCorrCut);
9f508565 1080
96387d94 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
ff36a82c 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
acd6f996 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
96387d94 1114 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1115 fOutputList->Add(fMultvsCentr);
9d63a2f6 1116
20791315 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
96387d94 1123
1124 //----------------------------------------------------------------------------
52a564ee 1125 EPVzA = new TH1D("EPVzA", "EPVzA", 60, -TMath::Pi()/2, TMath::Pi()/2);
9d63a2f6 1126 fOutputList->Add(EPVzA);
52a564ee 1127 EPVzC = new TH1D("EPVzC", "EPVzC", 60, -TMath::Pi()/2, TMath::Pi()/2);
9d63a2f6 1128 fOutputList->Add(EPVzC);
52a564ee 1129 EPTPC = new TH1D("EPTPC", "EPTPC", 60, -TMath::Pi()/2, TMath::Pi()/2);
9d63a2f6 1130 fOutputList->Add(EPTPC);
9f508565 1131
1132
52a564ee 1133 EPVz = new TH1D("EPVz", "EPVz", 60, -TMath::Pi()/2, TMath::Pi()/2);
9f508565 1134 fOutputList->Add(EPVz);
52a564ee 1135 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 60, -TMath::Pi()/2, TMath::Pi()/2);
9f508565 1136 fOutputList->Add(EPTPCp);
52a564ee 1137 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 60, -TMath::Pi()/2, TMath::Pi()/2);
9f508565 1138 fOutputList->Add(EPTPCn);
1139
96387d94 1140 //----------------------------------------------------------------------------
9d63a2f6 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}>))");
52a564ee 1145 fSubEventDPhiv2->Sumw2();
9d63a2f6 1146 fOutputList->Add(fSubEventDPhiv2);
9f508565 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}>))");
52a564ee 1154 fSubEventDPhiv2new->Sumw2();
9f508565 1155 fOutputList->Add(fSubEventDPhiv2new);
52a564ee 1156
1157
9f508565 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};
9d63a2f6 1161 // v2A, v2C, pt
1162 Int_t bins[3] = { 50, 50, nPtBins};
1163 Double_t xmin[3] = { -1., -1., 0};
9f508565 1164 Double_t xmax[3] = { 1., 1., 13};
9d63a2f6 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)");
52a564ee 1172 fV2Phi->Sumw2();
9d63a2f6 1173 fOutputList->Add(fV2Phi);
9f508565 1174
9f508565 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)");
52a564ee 1188 fV2Phivzerotot->Sumw2();
1189
9f508565 1190 fOutputList->Add(fV2Phivzerotot);
1191
1192
1193
96387d94 1194 //----------------------------------------------------------------------------
acd6f996 1195
9d63a2f6 1196 if(fPhiminusPsi){
96387d94 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);
9d63a2f6 1210 }
96387d94 1211 //----------------------------------------------------------------------------
1212 //----------------------------------------------------------------------------
9d63a2f6 1213 if(fpurity){
96387d94 1214 Int_t binsvpurity[3]={ 600,200, 200}; //pt, E/p,TPCnSigma
1215 Double_t xminvpurity[3]={0, 0, -10};
ff36a82c 1216 Double_t xmaxvpurity[3]={30, 2, 10};
96387d94 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);
9d63a2f6 1222 }
96387d94 1223 //----------------------------------------------------------------------------
1224
7fcd6b0f 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");
52a564ee 1233 fSparsephipsiULS->Sumw2();
7fcd6b0f 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");
52a564ee 1240 fSparsephipsiLS->Sumw2();
7fcd6b0f 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
52a564ee 1256 EPVzAftW = new TH1D("EPVzAftW", "EPVzAftW",60, -TMath::Pi()/2, TMath::Pi()/2);
1257 fOutputList->Add(EPVzAftW);
7fcd6b0f 1258
52a564ee 1259 fOutputList->Add(fHistEPDistrWeight);
1260
7fcd6b0f 1261
96387d94 1262 PostData(1,fOutputList);
1263 // create and post flowevent
1264 fFlowEvent = new AliFlowEvent(10000);
1265 PostData(2, fFlowEvent);
9d63a2f6 1266
1267 if(fSideBandsFlow){
96387d94 1268 fFlowEventCont = new AliFlowEvent(10000);
1269 PostData(3, fFlowEventCont);
9d63a2f6 1270 }
96387d94 1271}
914042c2 1272
1273//________________________________________________________________________
a70c9e97 1274void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
914042c2 1275{
96387d94 1276 // Info("Terminate");
1277 AliAnalysisTaskSE::Terminate();
914042c2 1278}
1279//_____________________________________________________________________________
a70c9e97 1280template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
914042c2 1281{
96387d94 1282 // QA multiplicity plots
1283 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1284 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
914042c2 1285}
1286//_____________________________________________________________________________
a70c9e97 1287template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
914042c2 1288{
96387d94 1289 //Set null cuts
914042c2 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//_____________________________________________________________________________
96387d94 1298void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
914042c2 1299{
96387d94 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);
914042c2 1306}
1307//_____________________________________________________________________________
a70c9e97 1308Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
914042c2 1309{
96387d94 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;
914042c2 1314}
1315//_________________________________________
a70c9e97 1316void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
914042c2 1317{
9f508565 1318 //============================Multiplicity TPV vs Global===============================================================================
9d63a2f6 1319 const Int_t nGoodTracks = event->GetNumberOfTracks();
9d63a2f6 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
f15c1f69 1323 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
1324 if(!trackAOD) AliFatal("Not a standard AOD");
9d63a2f6 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
f15c1f69 1331 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
1332 if(!trackAOD) AliFatal("Not a standard AOD");
9d63a2f6 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++;
9f508565 1341 } //track loop
1342 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1343 //============================================================================================================================
ff36a82c 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);
9f508565 1347 // cout << "--------------Centrality evaluated-------------------------"<<endl;
ff36a82c 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){
9f508565 1359 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1360 fCentralityBeforePileup->Fill(fCentrality);
1361 }//...after centrality selectrion
1362 //============================================================================================================================
1363 //to remove the bias introduced by multeplicity outliers---------------------
ff36a82c 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){
9f508565 1371 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1372 fCentralityAfterVZTRK->Fill(fCentrality);
ff36a82c 1373 }//...after centrality selectrion
9f508565 1374 //============================================================================================================================
1375 if(fMultCut){
47633986 1376 if(fTrigger==1 || fTrigger==4 || fTrigger==5){
9f508565 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 }
acd6f996 1391 if (centralitypass){
1392 fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1393 fCentralityAfterCorrCut->Fill(fCentrality);
1394 }//...after CORR CUT
9f508565 1395 //=================================All cuts are passed==================++++==================================================
1396 //=================================Now Centrality flattening for central trigger==================++++==================================================
41919e27 1397 if(fTrigger==0 || fTrigger==4){
9f508565 1398 if(!IsEventSelectedForCentrFlattening(fCentrality)){
edaa818d 1399 centralitypass = kFALSE;
9f508565 1400 fCentralityNoPassForFlattening->Fill(fCentrality);
1401 }
1402 }
1403 //==============================fill histo after all cuts==============================++++==================================================
9d63a2f6 1404 if(centralitypass){
96387d94 1405 fCentralityPass->Fill(fCentrality);
1406 fMultCorAfterCuts->Fill(multGlob, multTPC);
1407 fMultvsCentr->Fill(fCentrality, multTPC);
a70c9e97 1408 }
914042c2 1409}
1410//_____________________________________________________________________________
a70c9e97 1411void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
914042c2 1412{
96387d94 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;
914042c2 1417}
1418//_____________________________________________________________________________
a70c9e97 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)
914042c2 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//_____________________________________________________________________________
9f508565 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//---------------------------------------------------------------------------
52a564ee 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
9f508565 1553