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