]> 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
modified
[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
497 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
498
499 if (!aodTrack)
500 continue;
501
502 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
503 continue;
504
505 if (!aodTrack->TestFilterBit(128))
506 continue;
507
508
509 if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
510
511 Qx2p += TMath::Cos(2*aodTrack->Phi());
512 Qy2p += TMath::Sin(2*aodTrack->Phi());
513 }
514 if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
515
516 Qx2n += TMath::Cos(2*aodTrack->Phi());
517 Qy2n += TMath::Sin(2*aodTrack->Phi());
518 }
519
520
521 Qx2 += TMath::Cos(2*aodTrack->Phi());
522 Qy2 += TMath::Sin(2*aodTrack->Phi());
523
524
525
526
527 }
528
529 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
530 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
531 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
532
533 EPVzA->Fill(evPlAngV0A);
534 EPVzC->Fill(evPlAngV0C);
535 EPTPC->Fill(evPlAngTPC);
536
537 EPTPCn->Fill(evPlAngTPCn);
538 EPTPCp->Fill(evPlAngTPCp);
539 EPVz->Fill(evPlAngV0);
540
541 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
542 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
543 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
544
545
546 fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
547 fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
548 fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
549
550
551 //====================================================================================================================
552
553 AliAODTrack *track = NULL;
554
555 // Track loop
556 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
557 {
558 track = fAOD->GetTrack(iTracks);
559 if (!track)
560 {
561 printf("ERROR: Could not receive track %d\n", iTracks);
562 continue;
563 }
564 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
565
566 //--------------------------------------hfe begin-----------------------------------------------------------
567 //==========================================================================================================
568 //======================================track cuts==========================================================
569 if(track->Eta()<-0.8 || track->Eta()>0.8) continue; //eta cuts on candidates
570
571 // RecKine: ITSTPC cuts
572 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
573
574 // Reject kink mother
575 Bool_t kinkmotherpass = kTRUE;
576 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
577 if(track->GetID() == listofmotherkink[kinkmother]) {
578 kinkmotherpass = kFALSE;
579 continue;
580 }
581 }
582 if(!kinkmotherpass) continue;
583
584 // RecPrim
585 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
586 // HFEcuts: ITS layers cuts
587 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
588 // HFE cuts: TPC PID cleanup
589 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
590 //==========================================================================================================
591 Double_t eta = track->Eta();
592 Double_t phi = track->Phi();
47633986 593
594 if(fPhiCut){
595 if(phi<1.4 || phi >3.14)continue; //to have same EMCal phi acceptance
596 }
597
598
599
4c12b1ec
R		 600 Double_t pt = track->Pt(); //pt track after cuts
601 if(pt<fpTCutmin || pt>fpTCutmax) continue;
602 //==========================================================================================================
603 //=========================================PID==============================================================
604 if(track->GetTPCsignalN() < fTPCS) continue;
605 Float_t fITSnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron) : 1000;
606 Float_t fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
607 Float_t fTOFnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron) : 1000;
608 // Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
609 fITSnsigma->Fill(track->P(),fITSnSigma);
610 fTPCnsigma->Fill(track->P(),fTPCnSigma);
611 fTOFns->Fill(track->P(),fTOFnSigma);
612 fITSvsTOF->Fill(fTOFnSigma,fITSnSigma);
613 fTPCvsITS->Fill(fTPCnSigma,fITSnSigma);
614 fTPCvsTOF->Fill(fTPCnSigma,fTOFnSigma);
615
616 if( pt >= 0.3){
617 if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
618 }//cuts on nsigma tof full pt range
619
620 fITSnsigmaAftTOF->Fill(track->P(),fITSnSigma);
621 fTPCnsigmaAftTOF->Fill(track->P(),fTPCnSigma);
622 fTPCvsITSafterTOF->Fill(fTPCnSigma,fITSnSigma);
623
624 Double_t valPidSparse[3] = {
625 pt,
626 fITSnSigma,
627 fTPCnSigma,
628 };
629 fQAPidSparse->Fill(valPidSparse);
630
631
632 if( pt < 1.5){
633 if(fITSnSigma < fminITSnsigmaLowpT || fITSnSigma > fmaxITSnsigmaLowpT)continue;
634 }//cuts on nsigma its low pt
635 if( pt >= 1.5){
636 if(fITSnSigma < fminITSnsigmaHighpT || fITSnSigma > fmaxITSnsigmaHighpT)continue;
637 }//cuts on nsigma its high pt
638 fTPCnsigmaAftITSTOF->Fill(track->P(),fTPCnSigma);
639 if(pt >= 0.25 && pt < 1.5){
640 if(fTPCnSigma < fminTPCnsigmaLowpT || fTPCnSigma > fmaxTPCnsigmaLowpT) continue;
641 }//cuts on nsigma tpc lowpt
642 if( pt >= 1.5){
643 if(fTPCnSigma < fminTPCnsigmaHighpT || fTPCnSigma > fmaxTPCnsigmaHighpT) continue;
644 }//cuts on nsigma tpc high pt
645 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
646 fTPCnsigmaVSptAft->Fill(pt,fTPCnSigma);
647
648 //==========================================================================================================
649 //=========================================QA PID SPARSE====================================================
650 Float_t timeTOF = track->GetTOFsignal();
651 Double_t intTime[5] = {-99., -99., -99., -99., -99.};
652 track->GetIntegratedTimes(intTime);
653 Float_t timeElec = intTime[0];
654 Float_t intLength = 2.99792458e-2* timeElec;
655 Double_t beta = 0.1;
656 if ((intLength > 0) && (timeTOF > 0))
657 beta = intLength/2.99792458e-2/timeTOF;
658
659 // if(fQAPIDSparse){
660 // Double_t valPid[4] = {
661 // track->P(),
662 // track->GetTPCsignal(),
663 // beta,
664 // track->Charge()
665 // };
666 // fQAPid->Fill(valPid);
667 // }
668
669
670 fITSnsigmaAft->Fill(track->P(),fITSnSigma);
671 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
672 fTOFnsAft->Fill(track->P(),fTOFnSigma);
673 fTOFBetaAft->Fill(track->P(),beta);
674 fInclusiveElecPt->Fill(pt);
675 fPhi->Fill(phi);
676 fEta->Fill(eta);
677 //=========================================================================================================
678 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
679 AliFlowTrack *sTrack = new AliFlowTrack();
680 sTrack->Set(track);
681 sTrack->SetID(track->GetID());
682 sTrack->SetForRPSelection(kTRUE);
683 sTrack->SetForPOISelection(kTRUE);
684 sTrack->SetMass(263732);
685 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
686 {
687 // cout << " no of rps " << iRPs << endl;
688 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
689 if (!iRP) continue;
690 if (!iRP->InRPSelection()) continue;
691 if( sTrack->GetID() == iRP->GetID())
692 {
693 if(fDebug) printf(" was in RP set");
694 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
695 iRP->SetForRPSelection(kFALSE);
696 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
697 }
698 } //end of for loop on RPs
699 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
700 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
701 //============================Event Plane Method with V0====================================================
702 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
703 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
704 Double_t v2Phi[3] = {
705 v2PhiV0A,
706 v2PhiV0C,
707 pt};
708 fV2Phi->Fill(v2Phi);
709
710 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
711 Double_t v2PhiV0tot[2] = {
712 v2PhiVz,
713 pt};
714 fV2Phivzerotot->Fill(v2PhiV0tot);
715
716 //==========================================================================================================
717 //=========================================================================================================
718
719
720 if(fDCA){
721 fNonHFE = new AliSelectNonHFE();
722 fNonHFE->SetAODanalysis(kTRUE);
723 fNonHFE->SetInvariantMassCut(fInvmassCut);
724 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
725 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
726 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
727 fNonHFE->SetAlgorithm("DCA"); //KF
728 fNonHFE->SetPIDresponse(pidResponse);
729 fNonHFE->SetTrackCuts(-3,3);
730
731 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
732 fNonHFE->SetHistAngle(fOpeningAngleULS);
733 //fNonHFE->SetHistDCABack(fDCABack);
734 //fNonHFE->SetHistDCA(fDCA);
735 fNonHFE->SetHistMassBack(fInvmassLS1);
736 fNonHFE->SetHistMass(fInvmassULS1);
737
738 fNonHFE->FindNonHFE(iTracks,track,fAOD);
739
740 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
741 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
742 // Bool_t fUlsIsPartner = kFALSE;
743 // Bool_t fLsIsPartner = kFALSE;
744 if(fNonHFE->IsULS()){
745 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
746 fULSElecPt->Fill(track->Pt());
747 }
748 }
749
750 if(fNonHFE->IsLS()){
751 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
752 fLSElecPt->Fill(track->Pt());
753 }
754 }
755 }
756 if(!fDCA){
757 //=========================================================================================================
758 //----------------------Selection and Flow of Photonic Electrons-----------------------------
759 Bool_t fFlagPhotonicElec = kFALSE;
7fcd6b0f 760 SelectPhotonicElectron(iTracks,track,fTPCnSigma,evPlAngV0,fFlagPhotonicElec);
4c12b1ec
R		 761 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
762 // Semi inclusive electron
763 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
764 }
765 //-------------------------------------------------------------------------------------------
766
767 }//end loop on track
768 PostData(1, fOutputList);
769 PostData(2, fFlowEvent);
770
771 //----------hfe end---------
772}
773//_________________________________________
7fcd6b0f 774void AliAnalysisTaskFlowITSTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Float_t fTPCnSigma,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec)
4c12b1ec
R		 775{
776
777 //Identify non-heavy flavour electrons using Invariant mass method
778 Bool_t flagPhotonicElec = kFALSE;
779
780 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
781 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
782 if (!trackAsso) {
783 printf("ERROR: Could not receive track %d\n", jTracks);
784 continue;
785 }
786 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
787 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
7fcd6b0f 788
789 if(fAssoITSRefit){
790 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
791 }
792
793 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
794
795// if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
4c12b1ec
R		 796
797
798 if(jTracks == itrack) continue;
799 Double_t ptAsso=-999., nsigma=-999.0;
800 Double_t mass=-999., width = -999;
801 Double_t openingAngle = -999.;
802 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
803
804
805 ptAsso = trackAsso->Pt();
806 Short_t chargeAsso = trackAsso->Charge();
807 Short_t charge = track->Charge();
808 // nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
809 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
810
811
7fcd6b0f 812 //if(trackAsso->GetTPCNcls() < 80) continue;
813 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
4c12b1ec
R		 814 if(nsigma < -3 || nsigma > 3) continue;
815 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
164d5623 816 if(ptAsso < fptminAsso) continue;
4c12b1ec
R		 817
818 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
819 if(charge>0) fPDGe1 = -11;
820 if(chargeAsso>0) fPDGe2 = -11;
821
822 if(charge == chargeAsso) fFlagLS = kTRUE;
823 if(charge != chargeAsso) fFlagULS = kTRUE;
824
825 AliKFParticle::SetField(fAOD->GetMagneticField());
826 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
827 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
828 AliKFParticle recg(ge1, ge2);
829
830 if(recg.GetNDF()<1) continue;
831 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
832 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
833 recg.GetMass(mass,width);
834
835 openingAngle = ge1.GetAngle(ge2);
836 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
837 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
838 if(fOP_angle){if(openingAngle > fOpeningAngleCut) continue;}
839
840
841 if(fFlagLS) fInvmassLS1->Fill(mass);
842 if(fFlagULS) fInvmassULS1->Fill(mass);
843
7fcd6b0f 844 if(fFlagULS){
845 Double_t MassSparseULS[3] = {
846 track->Pt(),
847 mass
848 };
849 fSparseMassULS->Fill(MassSparseULS);
850 }
851 if(fFlagLS){
852 Double_t MassSparseLS[3] = {
853 track->Pt(),
854 mass
855 };
856 fSparseMassLS->Fill(MassSparseLS);
857 }
858
859
4c12b1ec
R		 860 if(mass<fInvmassCut){
861 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
862 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
863 }
864
7fcd6b0f 865
866 Double_t phi = track->Phi();
867 Float_t DeltaPhi_eEP = TVector2::Phi_0_2pi(phi - evPlAngV0);
868 if(DeltaPhi_eEP > TMath::Pi()) {DeltaPhi_eEP = DeltaPhi_eEP - TMath::Pi();}
869
870
871 if(mass<fInvmassCut){
872 if(fFlagULS){
873 Double_t ulsSparse[3] = {
874 track->Pt(),
875 fTPCnSigma,
876 DeltaPhi_eEP
877 };
878 fSparsephipsiULS->Fill(ulsSparse);
879 }
880 if(fFlagLS){
881 Double_t lsSparse[3] = {
882 track->Pt(),
883 fTPCnSigma,
884 DeltaPhi_eEP
885 };
886 fSparsephipsiLS->Fill(lsSparse);
887 }
888 }
889
890
891
4c12b1ec
R		 892 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
893 flagPhotonicElec = kTRUE;
894 }
895 }//track loop
896
897 fFlagPhotonicElec = flagPhotonicElec;
898}
899//___________________________________________
900void AliAnalysisTaskFlowITSTPCTOFQCSP::UserCreateOutputObjects()
901{
902
903 //Create histograms
904 //----------hfe initialising begin---------
905 fNullCuts = new AliFlowTrackCuts("null_cuts");
906
907 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
908 cc->SetNbinsMult(10000);
909 cc->SetMultMin(0);
910 cc->SetMultMax(10000);
911
912 cc->SetNbinsPt(100);
913 cc->SetPtMin(0);
914 cc->SetPtMax(5);
915
916 cc->SetNbinsPhi(180);
917 cc->SetPhiMin(0.0);
918 cc->SetPhiMax(TMath::TwoPi());
919
920 cc->SetNbinsEta(30);
921 cc->SetEtaMin(-8.0);
922 cc->SetEtaMax(+8.0);
923
924 cc->SetNbinsQ(500);
925 cc->SetQMin(0.0);
926 cc->SetQMax(3.0);
927
928
929 // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
930 // AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
931 // if (!inputHandler){
932 // AliFatal("Input handler needed");
933 // }
934 // else{
935 // fPIDResponse=inputHandler->GetPIDResponse();
936 // }
937 //pid response object
938 // if (!fPIDResponse)AliError("PIDResponse object was not created");
939
940
941 //--------Initialize PID
942 fPID->SetHasMCData(kFALSE);
943 if(!fPID->GetNumberOfPIDdetectors())
944 {
945 fPID->AddDetector("ITS", 0);
946 fPID->AddDetector("TOF", 1);
947 fPID->AddDetector("TPC", 2);
948
949 }
950
951 fPID->SortDetectors();
952 fPIDqa = new AliHFEpidQAmanager();
953 fPIDqa->Initialize(fPID);
954
955
956
957 //--------Initialize correction Framework and Cuts
958 fCFM = new AliCFManager;
959 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
960 fCFM->SetNStepParticle(kNcutSteps);
961 for(Int_t istep = 0; istep < kNcutSteps; istep++)
962 fCFM->SetParticleCutsList(istep, NULL);
963
964 if(!fCuts){
965 AliWarning("Cuts not available. Default cuts will be used");
966 fCuts = new AliHFEcuts;
967 fCuts->CreateStandardCuts();
968 }
969
970 fCuts->SetAOD();
971 fCuts->Initialize(fCFM);
972 //----------hfe initialising end--------
973 //---------Output Tlist
974 fOutputList = new TList();
975 fOutputList->SetOwner();
976 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
977
978 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
979 fOutputList->Add(fNoEvents);
980
981 fITSnsigma = new TH2F("fITSnsigma", "ITS - n sigma before HFE pid",600,0,6,400,-20,20);
982 fOutputList->Add(fITSnsigma);
983
984 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",600,0,6,400,-20,20);
985 fOutputList->Add(fTPCnsigma);
986
987 fITSnsigmaAft = new TH2F("fITSnsigmaAft", "ITS - n sigma after HFE pid",1000,0,10,300,-10,20);
988 fOutputList->Add(fITSnsigmaAft);
989 fITSvsTOF = new TH2F("fITSvsTOF", "ITS tof",400,-20,20,400,-20,20);
990 fOutputList->Add(fITSvsTOF);
991 fTPCvsITS = new TH2F("TPCvsITS", "TPC ITS",400,-20,20,400,-20,20);
992 fOutputList->Add(fTPCvsITS);
993 fTPCvsTOF = new TH2F("TPCvsTOF", "TPC TOF",400,-20,20,400,-20,20);
994 fOutputList->Add(fTPCvsTOF);
995 fTPCvsITSafterTOF = new TH2F("TPCvsITSafterTOF", "TPC ITS",400,-20,20,400,-20,20);
996 fOutputList->Add(fTPCvsITSafterTOF);
997
998
999 fITSnsigmaAftTOF = new TH2F("fITSnsigmaAftTOF", "ITS - n sigma after HFE pid",600,0,6,400,-20,20);
1000 fOutputList->Add(fITSnsigmaAftTOF);
1001
1002 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1003 fOutputList->Add(fTPCnsigmaAft);
1004
1005 fTPCnsigmaVSptAft = new TH2F("fTPCnsigmaVSptAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1006 fOutputList->Add(fTPCnsigmaVSptAft);
1007
1008 fTPCnsigmaAftITSTOF = new TH2F("fTPCnsigmaAftITSTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1009 fOutputList->Add(fTPCnsigmaAftITSTOF);
1010
1011 fTPCnsigmaAftTOF = new TH2F("fTPCnsigmaAftTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1012 fOutputList->Add(fTPCnsigmaAftTOF);
1013
1014 fTOFns = new TH2F("fTOFns","track TOFnSigma",600,0,6,400,-20,20);
1015 fOutputList->Add(fTOFns);
1016
1017 fTOFnsAft = new TH2F("fTOFnsAft","track TOFnSigma",600,0,6,400,-20,20);
1018 fOutputList->Add(fTOFnsAft);
1019
1020 fTOFBetaAft = new TH2F("fTOFBetaAft","track TOFBeta",600,0,6,120,0,1.2);
1021 fOutputList->Add(fTOFBetaAft);
1022
1023 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",100,0,5);
1024 fOutputList->Add(fInclusiveElecPt);
1025
1026 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,5);
1027 fOutputList->Add(fPhotoElecPt);
1028
1029 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,5);
1030 fOutputList->Add(fSemiInclElecPt);
1031
1032 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",100,0,5);
1033 fOutputList->Add(fULSElecPt);
1034
1035 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",100,0,5);
1036 fOutputList->Add(fLSElecPt);
1037
1038 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1039 fOutputList->Add(fInvmassLS1);
1040
1041 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1042 fOutputList->Add(fInvmassULS1);
1043
1044 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
1045 fOutputList->Add(fCentralityPass);
1046
1047 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
1048 fOutputList->Add(fCentralityNoPass);
1049
1050 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1051 fOutputList->Add(fCentralityNoPassForFlattening);
1052
1053 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1054 fOutputList->Add(fCentralityBeforePileup);
1055
1056 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1057 fOutputList->Add(fCentralityAfterVZTRK);
1058
1059 fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1060 fOutputList->Add(fCentralityAfterCorrCut);
1061
1062 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
1063 fOutputList->Add(fPhi);
1064
1065 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
1066 fOutputList->Add(fEta);
1067
1068 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
1069 fOutputList->Add(fVZEROA);
1070
1071 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
1072 fOutputList->Add(fVZEROC);
1073
1074 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
1075 fOutputList->Add(fTPCM);
1076
1077 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
1078 fOutputList->Add(fvertex);
1079
1080 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1081 fOutputList->Add(fMultCorBeforeCuts);
1082
1083 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1084 fOutputList->Add(fMultCorAfterCuts);
1085
1086 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1087 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
1088
1089 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1090 fOutputList->Add(fMultCorAfterVZTRKComp);
1091
1092 fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1093 fOutputList->Add(fMultCorAfterCorrCut);
1094
1095 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1096 fOutputList->Add(fMultvsCentr);
1097
1098 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1099 fOutputList->Add(fOpeningAngleLS);
1100
1101 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1102 fOutputList->Add(fOpeningAngleULS);
1103
1104
1105
1106 //----------------------------------------------------------------------------
1107 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1108 fOutputList->Add(EPVzA);
1109 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1110 fOutputList->Add(EPVzC);
1111 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1112 fOutputList->Add(EPTPC);
1113 EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1114 fOutputList->Add(EPVz);
1115 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1116 fOutputList->Add(EPTPCp);
1117 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1118 fOutputList->Add(EPTPCn);
1119
1120
1121 //----------------------------------------------------------------------------
1122 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1123 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1124 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1125 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1126 fOutputList->Add(fSubEventDPhiv2);
1127
1128 fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1129 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1130 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1131 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1132 fOutputList->Add(fSubEventDPhiv2new);
1133
1134 //================================Event Plane with VZERO=====================
1135 const Int_t nPtBins = 12;
1136 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};
1137 // v2A, v2C, pt
1138 Int_t bins[3] = { 50, 50, nPtBins};
1139 Double_t xmin[3] = { -1., -1., 0};
1140 Double_t xmax[3] = { 1., 1., 5};
1141 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1142 // Set bin limits for axes which are not standard binned
1143 fV2Phi->SetBinEdges(2, binsPt);
1144 // set axes titles
1145 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1146 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1147 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1148 fOutputList->Add(fV2Phi);
1149 //================================Event Plane with VZERO=====================
1150 Int_t binsV[2] = { 50, 100};
1151 Double_t xminV[2] = { -1., 0};
1152 Double_t xmaxV[2] = { 1., 5};
1153 fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1154 // Set bin limits for axes which are not standard binned
1155 //fV2Phivzerotot->SetBinEdges(1, binsPt);
1156 // set axes titles
1157 fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1158 fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1159 fOutputList->Add(fV2Phivzerotot);
1160
1161
1162
1163 //----------------------------------------------------------------------------
1164 //----------------------------------------------------------------------------
1165 // if(fQAPIDSparse){
1166 // Int_t binsQA[4] = { 150, 100, 120, 3};
1167 // Double_t xminQA[4] = { 0., 50, 0, -1.5};
1168 // Double_t xmaxQA[4] = { 15., 150, 1.2, 1.5};
1169 // fQAPid = new THnSparseF("fQAPid", "p:dEdx:beta:ch", 4, binsQA, xminQA, xmaxQA);
1170 // fQAPid->GetAxis(0)->SetTitle("p (Gev/c");
1171 // fQAPid->GetAxis(1)->SetTitle("dE/dx");
1172// fQAPid->GetAxis(2)->SetTitle("#beta (TOF)");
1173// fQAPid->GetAxis(3)->SetTitle("charge");
1174// fOutputList->Add(fQAPid);
1175// }
1176 //===========================================================================
1177 Int_t binsQA2[3] = { 100, 40, 150/*, 60*/};
1178 Double_t xminQA2[3] = { 0., -2, -15/*, -3*/};
1179 Double_t xmaxQA2[3] = { 5., 2, 15/*, 3*/};
1180 fQAPidSparse = new THnSparseF("fQAPidSparse", "pt:itsnsigma:tpcnsigma", 3, binsQA2, xminQA2, xmaxQA2);
1181 fQAPidSparse->GetAxis(0)->SetTitle("pt (Gev/c)");
1182 fQAPidSparse->GetAxis(1)->SetTitle("itsnsigma");
1183 fQAPidSparse->GetAxis(2)->SetTitle("tpcnsigma");
1184 fOutputList->Add(fQAPidSparse);
1185 //===========================================================================
7fcd6b0f 1186
1187
1188
1189 Int_t binsphipsi[3] = { 100, 150, 6};
1190 Double_t xminphipsi[3] = { 0., -15, 0};
1191 Double_t xmaxphipsi[3] = { 5., 15, TMath::Pi()};
1192 fSparsephipsiULS = new THnSparseF("fSparsephipsiULS", "pt:tpcnsigma:DeltaPhiULS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1193 fSparsephipsiULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1194 fSparsephipsiULS->GetAxis(1)->SetTitle("tpcnsigma");
1195 fSparsephipsiULS->GetAxis(2)->SetTitle("DeltaPhiULS");
1196 fOutputList->Add(fSparsephipsiULS);
1197
1198 fSparsephipsiLS = new THnSparseF("fSparsephipsiLS", "pt:tpcnsigma:DeltaPhiLS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1199 fSparsephipsiLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1200 fSparsephipsiLS->GetAxis(1)->SetTitle("tpcnsigma");
1201 fSparsephipsiLS->GetAxis(2)->SetTitle("DeltaPhiLS");
1202 fOutputList->Add(fSparsephipsiLS);
1203
1204
1205 Int_t binsmass[2] = { 100, 200};
1206 Double_t xminmass[2] = { 0., 0};
1207 Double_t xmaxmass[2] = { 5., 1.};
1208 fSparseMassULS = new THnSparseF("fSparseMassULS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1209 fSparseMassULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1210 fSparseMassULS->GetAxis(1)->SetTitle("mass");
1211 fOutputList->Add(fSparseMassULS);
1212
1213 fSparseMassLS = new THnSparseF("fSparseMassLS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1214 fSparseMassLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1215 fSparseMassLS->GetAxis(1)->SetTitle("mass");
1216 fOutputList->Add(fSparseMassLS);
1217
1218
4c12b1ec
R		 1219 PostData(1,fOutputList);
1220 // create and post flowevent
1221 fFlowEvent = new AliFlowEvent(10000);
1222 PostData(2, fFlowEvent);
1223
1224}
1225//________________________________________________________________________
1226void AliAnalysisTaskFlowITSTPCTOFQCSP::Terminate(Option_t *)
1227{
1228 // Info("Terminate");
1229 AliAnalysisTaskSE::Terminate();
1230}
1231//_____________________________________________________________________________
1232template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::PlotVZeroMultiplcities(const T* event) const
1233{
1234 // QA multiplicity plots
1235 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1236 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1237}
1238//_____________________________________________________________________________
1239template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::SetNullCuts(T* event)
1240{
1241 //Set null cuts
1242 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1243 fCutsRP->SetEvent(event, MCEvent());
1244 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1245 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1246 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1247 fNullCuts->SetEvent(event, MCEvent());
1248}
1249//_____________________________________________________________________________
1250void AliAnalysisTaskFlowITSTPCTOFQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1251{
1252 //Prepare flow events
1253 FlowEv->ClearFast();
1254 FlowEv->Fill(fCutsRP, fNullCuts);
1255 FlowEv->SetReferenceMultiplicity(iMulti);
1256 FlowEv->DefineDeadZone(0, 0, 0, 0);
1257 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1258}
1259//_____________________________________________________________________________
1260Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1261{
1262 // Check single track cuts for a given cut step
1263 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1264 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1265 return kTRUE;
1266}
1267//_________________________________________
1268void AliAnalysisTaskFlowITSTPCTOFQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
1269{
1270 //============================Multiplicity TPV vs Global===============================================================================
1271 const Int_t nGoodTracks = event->GetNumberOfTracks();
1272 Float_t multTPC(0.); // tpc mult estimate
1273 Float_t multGlob(0.); // global multiplicity
1274 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1275 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1276 if (!trackAOD) continue;
1277 if (!(trackAOD->TestFilterBit(1))) continue;
1278 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;
1279 multTPC++;
1280 }
1281 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1282 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1283 if (!trackAOD) continue;
1284 if (!(trackAOD->TestFilterBit(16))) continue;
1285 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;
1286 Double_t b[2] = {-99., -99.};
1287 Double_t bCov[3] = {-99., -99., -99.};
1288 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1289 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1290 multGlob++;
1291 } //track loop
1292 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1293 //============================================================================================================================
1294 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1295 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1296 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1297 // cout << "--------------Centrality evaluated-------------------------"<<endl;
1298 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1299 {
1300 fCentralityNoPass->Fill(fCentrality);
1301 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1302 centralitypass = kFALSE;
1303 }else
1304 {
1305 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1306 centralitypass = kTRUE;
1307 }
1308 if (centralitypass){
1309 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1310 fCentralityBeforePileup->Fill(fCentrality);
1311 }//...after centrality selectrion
1312 //============================================================================================================================
1313 //to remove the bias introduced by multeplicity outliers---------------------
1314 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1315 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1316 if (TMath::Abs(centv0 - centTrk) > 5.0){
1317 centralitypass = kFALSE;
1318 fCentralityNoPass->Fill(fCentrality);
1319 }
1320 if (centralitypass){
1321 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1322 fCentralityAfterVZTRK->Fill(fCentrality);
1323 }//...after centrality selectrion
1324 //============================================================================================================================
1325 if(fMultCut){
47633986 1326 if(fTrigger==1 || fTrigger==4 || fTrigger==5){
4c12b1ec
R		 1327 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1328 // cout <<" Trigger ==" <<fTrigger<< endl;
1329 centralitypass = kFALSE;
1330 fCentralityNoPass->Fill(fCentrality);
1331 }//2011 Semicentral
1332 }
1333 if(fTrigger==0){
1334 if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1335 // cout <<" Trigger ==" <<fTrigger<< endl;
1336 centralitypass = kFALSE;
1337 fCentralityNoPass->Fill(fCentrality);
1338 }//2011
1339 }//2011 Central
1340 }
1341 if (centralitypass){
1342 fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1343 fCentralityAfterCorrCut->Fill(fCentrality);
1344 }//...after CORR CUT
1345 //=================================All cuts are passed==================++++==================================================
1346 //=================================Now Centrality flattening for central trigger==================++++==================================================
41919e27 1347 if(fTrigger==0 || fTrigger==4){
4c12b1ec
R		 1348 if(!IsEventSelectedForCentrFlattening(fCentrality)){
1349 centralitypass = kFALSE;
1350 fCentralityNoPassForFlattening->Fill(fCentrality);
1351 }
1352 }
1353 //==============================fill histo after all cuts==============================++++==================================================
1354 if(centralitypass){
1355 fCentralityPass->Fill(fCentrality);
1356 fMultCorAfterCuts->Fill(multGlob, multTPC);
1357 fMultvsCentr->Fill(fCentrality, multTPC);
1358 }
1359}
1360//_____________________________________________________________________________
1361void AliAnalysisTaskFlowITSTPCTOFQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1362{
1363 // Set a centrality range ]min, max] and define the method to use for centrality selection
1364 fCentralityMin = CentralityMin;
1365 fCentralityMax = CentralityMax;
1366 fkCentralityMethod = CentralityMethod;
1367}
1368//_____________________________________________________________________________
1369void 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)
1370{
1371 //Set PID cuts
1372 fminTOFnSigma = minTOFnSigma;
1373 fmaxTOFnSigma = maxTOFnSigma;
1374 fminITSnsigmaLowpT = minITSnsigmalowpt;
1375 fmaxITSnsigmaLowpT = maxITSnsigmalowpt;
1376 fminITSnsigmaHighpT = minITSnsigmahighpt;
1377 fmaxITSnsigmaHighpT = maxITSnsigmahighpt;
1378 fminTPCnsigmaLowpT = minTPCnsigmalowpt;
1379 fmaxTPCnsigmaLowpT = maxTPCnsigmalowpt;
1380 fminTPCnsigmaHighpT = minTPCnsigmahighpt;
1381 fmaxTPCnsigmaHighpT = maxTPCnsigmahighpt;
1382
1383}
1384//_____________________________________________________________________________
1385//_____________________________________________________________________________
1386void AliAnalysisTaskFlowITSTPCTOFQCSP::SetpTCuttrack(Double_t ptmin, Double_t ptmax)
1387{
1388 //Set pt cuts
1389 fpTCutmin = ptmin;
1390 fpTCutmax = ptmax;
164d5623 1391}
4c12b1ec
R		 1392//_____________________________________________________________________________
1393//_____________________________________________________________________________
1394void AliAnalysisTaskFlowITSTPCTOFQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1395 // set the histo for centrality flattening
1396 // the centrality is flatten in the range minCentr,maxCentr
1397 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1398 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1399 // 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)
1400 // 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
1401
1402 if(maxCentr<minCentr){
1403 AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1404 }
1405
1406 if(fHistCentrDistr)delete fHistCentrDistr;
1407 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1408 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1409 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1410 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1411 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1412 Double_t ref=0.,bincont=0.,binrefwidth=1.;
1413 Int_t binref=0;
1414 if(TMath::Abs(centrRef)<0.0001){
1415 binref=fHistCentrDistr->GetMinimumBin();
1416 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1417 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1418 }
1419 else if(centrRef>0.){
1420 binref=h->FindBin(centrRef);
1421 if(binref<1||binref>h->GetNbinsX()){
1422 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1423 }
1424 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1425 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1426 }
1427 else{
1428 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1429 binref=fHistCentrDistr->GetMaximumBin();
1430 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1431 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1432 }
1433
1434 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1435 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1436 bincont=h->GetBinContent(j);
1437 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1438 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1439 }
1440 else{
1441 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1442 }
1443 }
1444
1445 fHistCentrDistr->SetBinContent(0,switchTRand);
1446 return;
1447
1448}
1449
1450//-------------------------------------------------
1451Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1452 //
1453 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1454 // Can be faster if it was required that fHistCentrDistr covers
1455 // exactly the desired centrality range (e.g. part of the lines below should be done during the
1456 // setting of the histo) and TH1::SetMinimum called
1457 //
1458
1459 if(!fHistCentrDistr) return kTRUE;
1460 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1461 // if(maxbin>fHistCentrDistr->GetNbinsX()){
1462 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1463 // }
1464
1465 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1466 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1467 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1468
1469 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1470 if(gRandom->Uniform(1.)<bincont)return kTRUE;
1471 return kFALSE;
1472 }
1473
1474 if(centDigits*100.<bincont)return kTRUE;
1475 return kFALSE;
1476
1477}
1478//---------------------------------------------------------------------------
1479
1480
1481//_____________________________________________________________________________
1482