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