]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added pile up protection
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalQCSP.cxx
CommitLineData
914042c2 1/**************************************************************************
96387d94 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
914042c2 14 **************************************************************************/
15
16
17////////////////////////////////////////////////////////////////////////
18// //
9d63a2f6 19// Task for Heavy Flavour Electron Flow with TPC plus EMCal //
914042c2 20// Non-Photonic Electron identified with Invariant mass //
21// analysis methos in function SelectPhotonicElectron //
96387d94 22// //
914042c2 23// //
24// Author: Andrea Dubla (Utrecht University) //
96387d94 25// //
914042c2 26// //
27////////////////////////////////////////////////////////////////////////
28
29#include "TChain.h"
30#include "TTree.h"
31#include "TH2F.h"
32#include "TMath.h"
33#include "TCanvas.h"
34#include "THnSparse.h"
35#include "TLorentzVector.h"
36#include "TString.h"
37#include "TFile.h"
38#include "AliAnalysisTask.h"
39#include "AliAnalysisManager.h"
40#include "AliESDEvent.h"
41#include "AliESDHandler.h"
42#include "AliAODEvent.h"
43#include "AliAODHandler.h"
a70c9e97 44#include "AliAnalysisTaskFlowTPCEMCalQCSP.h"
914042c2 45#include "TGeoGlobalMagField.h"
46#include "AliLog.h"
47#include "AliAnalysisTaskSE.h"
48#include "TRefArray.h"
49#include "TVector.h"
50#include "AliESDInputHandler.h"
51#include "AliESDpid.h"
52#include "AliAODInputHandler.h"
53#include "AliAODPid.h"
54#include "AliESDtrackCuts.h"
55#include "AliPhysicsSelection.h"
56#include "AliCentralitySelectionTask.h"
57#include "AliESDCaloCluster.h"
58#include "AliAODCaloCluster.h"
59#include "AliESDCaloTrigger.h"
60#include "AliEMCALRecoUtils.h"
61#include "AliEMCALGeometry.h"
62#include "AliGeomManager.h"
63#include "stdio.h"
64#include "TGeoManager.h"
65#include "iostream"
66#include "fstream"
67#include "AliEMCALTrack.h"
68//#include "AliEMCALTracker.h"
69#include "AliMagF.h"
70#include "AliKFParticle.h"
71#include "AliKFVertex.h"
72#include "AliPID.h"
73#include "AliPIDResponse.h"
74#include "AliHFEcontainer.h"
75#include "AliHFEcuts.h"
76#include "AliHFEpid.h"
77#include "AliHFEpidBase.h"
78#include "AliHFEpidQAmanager.h"
79#include "AliHFEtools.h"
80#include "AliCFContainer.h"
81#include "AliCFManager.h"
82#include "AliKFParticle.h"
83#include "AliKFVertex.h"
84#include "AliCentrality.h"
85#include "AliVEvent.h"
86#include "AliStack.h"
87#include "AliMCEvent.h"
88#include "TProfile.h"
89#include "AliFlowCandidateTrack.h"
90#include "AliFlowTrackCuts.h"
91#include "AliFlowEventSimple.h"
92#include "AliFlowCommonConstants.h"
93#include "AliFlowEvent.h"
94#include "TVector3.h"
95#include "TRandom2.h"
96#include "AliESDVZERO.h"
97#include "AliAODVZERO.h"
98#include "AliPID.h"
99#include "AliPIDResponse.h"
100#include "AliFlowTrack.h"
101#include "AliAnalysisTaskVnV0.h"
20791315 102#include "AliSelectNonHFE.h"
914042c2 103
9d63a2f6 104
914042c2 105class AliFlowTrackCuts;
106
9d63a2f6 107using namespace std;
914042c2 108
a70c9e97 109ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
914042c2 110//________________________________________________________________________
96387d94 111AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
112: AliAnalysisTaskSE(name)
9d63a2f6 113,fDebug(0)
114,fAOD(0)
115,fGeom(0)
116,fOutputList(0)
117,fCuts(0)
118,fIdentifiedAsOutInz(kFALSE)
119,fPassTheEventCut(kFALSE)
120,fCFM(0)
121,fPID(0)
122,fPIDqa(0)
123,fCutsRP(0) // track cuts for reference particles
124,fNullCuts(0) // dummy cuts for flow event tracks
125,fFlowEvent(0) //! flow events (one for each inv mass band)
126,fkCentralityMethod(0)
127,fCentrality(0)
128,fCentralityMin(0)
129,fCentralityMax(0)
130,fInvmassCut(0)
040cc6fa 131,fpTCut(0)
9d63a2f6 132,fTrigger(0)
133,fPhi(0)
134,fEta(0)
135,fVZEROA(0)
136,fVZEROC(0)
137,fTPCM(0)
138,fNoEvents(0)
139,fTrkEovPBef(0)
140//,fdEdxBef(0)
141,fInclusiveElecPt(0)
142,fTPCnsigma(0)
143,fTPCnsigmaAft(0)
144,fCentralityPass(0)
145,fCentralityNoPass(0)
146,fInvmassLS1(0)
147,fInvmassULS1(0)
148,fPhotoElecPt(0)
149,fSemiInclElecPt(0)
150,fULSElecPt(0)
151,fLSElecPt(0)
152,fminTPC(-1)
153,fmaxTPC(3)
154,fminEovP(0.8)
155,fmaxEovP(1.2)
156,fminM20(0.03)
157,fmaxM20(0.3)
158,fminM02(0.03)
159,fmaxM02(0.5)
160,fDispersion(1)
161,fMultCorAfterCuts(0)
162,fMultvsCentr(0)
163,fSubEventDPhiv2(0)
164,EPVzA(0)
165,EPVzC(0)
166,EPTPC(0)
167,fV2Phi(0)
168,fSparseElectronHadron(0)
169,fvertex(0)
170,fMultCorBeforeCuts(0)
171,fSideBandsFlow(kFALSE)
172,fPhiminusPsi(kFALSE)
173,fFlowEventCont(0) //! flow events (one for each inv mass band)
174,fpurity(kFALSE)
175,fSparseElectronpurity(0)
20791315 176,fOpeningAngleLS(0)
177,fOpeningAngleULS(0)
178,fNonHFE(new AliSelectNonHFE)
179,fDCA(0)
180,fOpeningAngleCut(0)
181,fOP_angle(0)
8e8d9d2a 182,fAssoTPCCluster(0)
183,fAssoITSRefit(0)
ff36a82c 184,fMultCut(0)
185,fMultCorAfterCentrBeforeCuts(0)
186,fMultCorAfterVZTRKComp(0)
187,fCentralityBeforePileup(0)
188,fCentralityAfterVZTRK(0)
914042c2 189{
96387d94 190 //Named constructor
191
192 fPID = new AliHFEpid("hfePid");
193 // Define input and output slots here
194 // Input slot #0 works with a TChain
195 DefineInput(0, TChain::Class());
196 // Output slot #0 id reserved by the base class for AOD
197 // Output slot #1 writes into a TH1 container
198 // DefineOutput(1, TH1I::Class());
199 DefineOutput(1, TList::Class());
200 DefineOutput(2, AliFlowEventSimple::Class());
201 if(fSideBandsFlow){
202 DefineOutput(3, AliFlowEventSimple::Class());
203 }
204 // DefineOutput(3, TTree::Class());
914042c2 205}
206
207//________________________________________________________________________
96387d94 208AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
209: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
9d63a2f6 210,fDebug(0)
211,fAOD(0)
212,fGeom(0)
213,fOutputList(0)
214,fCuts(0)
215,fIdentifiedAsOutInz(kFALSE)
216,fPassTheEventCut(kFALSE)
217,fCFM(0)
218,fPID(0)
219,fPIDqa(0)
220,fCutsRP(0) // track cuts for reference particles
221,fNullCuts(0) // dummy cuts for flow event tracks
222,fFlowEvent(0) //! flow events (one for each inv mass band)
223,fkCentralityMethod(0)
224,fCentrality(0)
225,fCentralityMin(0)
226,fCentralityMax(0)
227,fInvmassCut(0)
040cc6fa 228,fpTCut(0)
9d63a2f6 229,fTrigger(0)
230,fPhi(0)
231,fEta(0)
232,fVZEROA(0)
233,fVZEROC(0)
234,fTPCM(0)
235,fNoEvents(0)
236,fTrkEovPBef(0)
237//,fdEdxBef(0)
238,fInclusiveElecPt(0)
239,fTPCnsigma(0)
240,fTPCnsigmaAft(0)
241,fCentralityPass(0)
242,fCentralityNoPass(0)
243,fInvmassLS1(0)
244,fInvmassULS1(0)
245,fPhotoElecPt(0)
246,fSemiInclElecPt(0)
247,fULSElecPt(0)
248,fLSElecPt(0)
249,fminTPC(-1)
250,fmaxTPC(3)
251,fminEovP(0.8)
252,fmaxEovP(1.2)
253,fminM20(0.03)
254,fmaxM20(0.3)
255,fminM02(0.03)
256,fmaxM02(0.5)
257,fDispersion(1)
258,fMultCorAfterCuts(0)
259,fMultvsCentr(0)
260,fSubEventDPhiv2(0)
261,EPVzA(0)
262,EPVzC(0)
263,EPTPC(0)
264,fV2Phi(0)
265,fSparseElectronHadron(0)
266,fvertex(0)
267,fMultCorBeforeCuts(0)
268,fSideBandsFlow(kFALSE)
269,fPhiminusPsi(kFALSE)
270,fFlowEventCont(0) //! flow events (one for each inv mass band)
271,fpurity(kFALSE)
272,fSparseElectronpurity(0)
20791315 273,fOpeningAngleLS(0)
274,fOpeningAngleULS(0)
275,fNonHFE(new AliSelectNonHFE)
276,fDCA(0)
277,fOpeningAngleCut(0)
278,fOP_angle(0)
8e8d9d2a 279,fAssoTPCCluster(0)
280,fAssoITSRefit(0)
ff36a82c 281,fMultCut(0)
282,fMultCorAfterCentrBeforeCuts(0)
283,fMultCorAfterVZTRKComp(0)
284,fCentralityBeforePileup(0)
285,fCentralityAfterVZTRK(0)
914042c2 286{
96387d94 287 //Default constructor
288 fPID = new AliHFEpid("hfePid");
289 // Constructor
290 // Define input and output slots here
291 // Input slot #0 works with a TChain
292 DefineInput(0, TChain::Class());
293 // Output slot #0 id reserved by the base class for AOD
294 // Output slot #1 writes into a TH1 container
295 // DefineOutput(1, TH1I::Class());
296 DefineOutput(1, TList::Class());
297 DefineOutput(2, AliFlowEventSimple::Class());
9d63a2f6 298 // DefineOutput(3, TTree::Class());
96387d94 299 if(fSideBandsFlow){
300 DefineOutput(3, AliFlowEventSimple::Class());
301 }
9d63a2f6 302 //DefineOutput(3, TTree::Class());
914042c2 303}
304//_________________________________________
305
a70c9e97 306AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
914042c2 307{
96387d94 308 //Destructor
309
310 delete fOutputList;
311 delete fGeom;
312 delete fPID;
313 delete fCFM;
314 delete fPIDqa;
315 if (fDCA) delete fNonHFE;
316 if (fOutputList) delete fOutputList;
317 if (fFlowEvent) delete fFlowEvent;
318 if (fFlowEventCont) delete fFlowEventCont;
319
914042c2 320}
321//_________________________________________
322
a70c9e97 323void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
914042c2 324{
96387d94 325 //Main loop
326 //Called for each event
327
328 // create pointer to event
329
330 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
331 if (!fAOD)
332 {
333 printf("ERROR: fAOD not available\n");
334 return;
335 }
336
337 if(!fCuts)
338 {
339 AliError("HFE cuts not available");
340 return;
341 }
342
343 if(!fPID->IsInitialized())
344 {
345 // Initialize PID with the given run number
346 AliWarning("PID not initialised, get from Run no");
347 fPID->InitializePID(fAOD->GetRunNumber());
348 }
914042c2 349
edaa818d 350 // cout << "kTrigger == " << fTrigger <<endl;
914042c2 351
352 if(fTrigger==0){
96387d94 353 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
914042c2 354 }
355 if(fTrigger==1){
96387d94 356 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
914042c2 357 }
358 if(fTrigger==2){
96387d94 359 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
914042c2 360 }
361 if(fTrigger==3){
96387d94 362 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
914042c2 363 }
9d63a2f6 364
96387d94 365
366 //---------------CENTRALITY AND EVENT SELECTION-----------------------
367 Int_t fNOtrks = fAOD->GetNumberOfTracks();
9d63a2f6 368 Float_t vtxz = -999;
369 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
370 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
371 TString vtxTtl = trkVtx->GetTitle();
372 if (!vtxTtl.Contains("VertexerTracks"))return;
373 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
374 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
375 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
376 vtxz = trkVtx->GetZ();
377 if(TMath::Abs(vtxz)>10)return;
96387d94 378
379 // Event cut
9d63a2f6 380 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
381 if(fNOtrks<2) return;
382
ff36a82c 383 Bool_t pass = kFALSE; //to select centrality and pile up protection
9d63a2f6 384 CheckCentrality(fAOD,pass);
385 if(!pass)return;
ff36a82c 386 fvertex->Fill(vtxz);
a70c9e97 387
96387d94 388 fNoEvents->Fill(0);
389 PlotVZeroMultiplcities(fAOD);
390
391 SetNullCuts(fAOD);
392 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
393
9d63a2f6 394 if(fSideBandsFlow){
96387d94 395 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont); //Calculate event plane Qvector and EP resolution for inclusive
9d63a2f6 396 }
397
96387d94 398 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
399 if(!pidResponse)
400 {
401 AliDebug(1, "Using default PID Response");
402 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
914042c2 403 }
914042c2 404
96387d94 405 fPID->SetPIDResponse(pidResponse);
406 fCFM->SetRecEventInfo(fAOD);
407
408 // Look for kink mother
409 Int_t numberofvertices = fAOD->GetNumberOfVertices();
410 Double_t listofmotherkink[numberofvertices];
411 Int_t numberofmotherkink = 0;
412 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
413 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
414 if(!aodvertex) continue;
415 if(aodvertex->GetType()==AliAODVertex::kKink) {
416 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
417 if(!mother) continue;
418 Int_t idmother = mother->GetID();
419 listofmotherkink[numberofmotherkink] = idmother;
420 //printf("ID %d\n",idmother);
421 numberofmotherkink++;
422 }
423 }
424
425 //=============================================V0EP from Alex======================================================================
9d63a2f6 426 Double_t qxEPa = 0, qyEPa = 0;
427 Double_t qxEPc = 0, qyEPc = 0;
428
429 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
430 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
431
432
433 Double_t Qx2 = 0, Qy2 = 0;
434
435 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
436
437 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
438
439 if (!aodTrack)
440 continue;
441
442 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
443 continue;
444
445 if (!aodTrack->TestFilterBit(128))
446 continue;
447
448 Qx2 += TMath::Cos(2*aodTrack->Phi());
449 Qy2 += TMath::Sin(2*aodTrack->Phi());
450 }
451
452 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
453
454 EPVzA->Fill(evPlAngV0A);
455 EPVzC->Fill(evPlAngV0C);
456 EPTPC->Fill(evPlAngTPC);
457
458 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
459 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
460 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
96387d94 461 //====================================================================================================================
9d63a2f6 462
463
96387d94 464 AliAODTrack *track = NULL;
465
466 // Track loop
467 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
468 {
469 track = fAOD->GetTrack(iTracks);
470 if (!track)
914042c2 471 {
96387d94 472 printf("ERROR: Could not receive track %d\n", iTracks);
473 continue;
914042c2 474 }
96387d94 475
476 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
477 //----------hfe begin---------
478 if(track->Eta()<-0.7 || track->Eta()>0.7) continue; //eta cuts on candidates
479
480 // RecKine: ITSTPC cuts
481 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
482
483 // Reject kink mother
484 Bool_t kinkmotherpass = kTRUE;
485 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
486 if(track->GetID() == listofmotherkink[kinkmother]) {
487 kinkmotherpass = kFALSE;
488 continue;
489 }
490 }
491 if(!kinkmotherpass) continue;
492
493 // RecPrim
494 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
495 // HFEcuts: ITS layers cuts
496 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
497 // HFE cuts: TPC PID cleanup
498 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
499
500 Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
501 // Track extrapolation
502 Int_t fClsId = track->GetEMCALcluster();
503 if(fClsId < 0) continue;
504 AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
505 if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
506
507 pt = track->Pt(); //pt track after cuts
508 if(pt<fpTCut) continue;
509 fClsE = cluster->E();
510 p = track->P();
511 // dEdx = track->GetTPCsignal();
512 fEovP = fClsE/p;
513 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
514 Double_t m20 =cluster->GetM20();
515 Double_t m02 =cluster->GetM02();
516 Double_t disp=cluster->GetDispersion();
517 if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
518 fTPCnsigma->Fill(p,fTPCnSigma);
519 // fdEdxBef->Fill(p,dEdx);
520 Double_t eta = track->Eta();
521 Double_t phi = track->Phi();
522 //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
523 //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
524 Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
525 if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
526 Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
527 if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
528
529 if(fPhiminusPsi){
530 Double_t valueElh[8] = {
531 pt,
532 fEovP,
533 fTPCnSigma,
534 m20,
535 m02,
536 disp,
537 dPhi_aeh,
538 dPhi_ceh};
539 fSparseElectronHadron->Fill(valueElh);
540 }
541 //----------------------------------------------------------------------------------------------------------
542 //---------------------------From here usual electron selection---------------------------------------------
543 //----------------------------------------------------------------------------------------------------------
544 if(m20 < fminM20 || m20 > fmaxM20) continue;
545 if(m02 < fminM02 || m02 > fmaxM02) continue;
546 if(disp > fDispersion ) continue;
547 //---------------------------------for purity---------------------------------------------------------------
548 if(fpurity){
549 Double_t valuepurity[3] = {
550 pt,
551 fEovP,
552 fTPCnSigma};
ff36a82c 553 fSparseElectronpurity->Fill(valuepurity);
96387d94 554 }
555 //----------------------------------------------------------------------------------------------------------
556 //----------------------------------------------------------------------------------------------------------
557 if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
558 //===============================Flow Event for Contamination=============================================
559 if(fSideBandsFlow){
560 if(fEovP>0 && fEovP<0.6){
561 AliFlowTrack *sTrackCont = new AliFlowTrack();
562 sTrackCont->Set(track);
563 sTrackCont->SetID(track->GetID());
564 sTrackCont->SetForRPSelection(kFALSE);
565 sTrackCont->SetForPOISelection(kTRUE);
566 sTrackCont->SetMass(2637);
567 for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
568 {
569 // cout << " no of rps " << iRPs << endl;
570 AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
571 if (!iRPCont) continue;
572 if (!iRPCont->InRPSelection()) continue;
573 if( sTrackCont->GetID() == iRPCont->GetID())
574 {
575 if(fDebug) printf(" was in RP set");
576 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
577 iRPCont->SetForRPSelection(kFALSE);
578 fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
579 }
580 } //end of for loop on RPs
581 fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
582 }
583 }
584 //==========================================================================================================
585 //===============================From here eovP cut is used fro QC, SP and EPV0=============================
586 if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
587 //==========================================================================================================
588 //============================Event Plane Method with V0====================================================
589 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
590 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
591 Double_t v2Phi[3] = {
592 v2PhiV0A,
593 v2PhiV0C,
594 pt};
595 fV2Phi->Fill(v2Phi);
596 //=========================================================================================================
597 fTPCnsigmaAft->Fill(p,fTPCnSigma);
598 fInclusiveElecPt->Fill(pt);
599 fPhi->Fill(phi);
600 fEta->Fill(eta);
601 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
602 AliFlowTrack *sTrack = new AliFlowTrack();
603 sTrack->Set(track);
604 sTrack->SetID(track->GetID());
605 sTrack->SetForRPSelection(kFALSE);
606 sTrack->SetForPOISelection(kTRUE);
607 sTrack->SetMass(263732);
608 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
609 {
610 // cout << " no of rps " << iRPs << endl;
611 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
612 if (!iRP) continue;
613 if (!iRP->InRPSelection()) continue;
614 if( sTrack->GetID() == iRP->GetID())
615 {
616 if(fDebug) printf(" was in RP set");
617 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
618 iRP->SetForRPSelection(kFALSE);
619 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
620 }
621 } //end of for loop on RPs
622 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
623
624
625
626 if(fDCA){
627 //----------------------Selection of Photonic Electrons DCA-----------------------------
628 fNonHFE = new AliSelectNonHFE();
629 fNonHFE->SetAODanalysis(kTRUE);
630 fNonHFE->SetInvariantMassCut(fInvmassCut);
631 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
632 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
633 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
634 fNonHFE->SetAlgorithm("DCA"); //KF
635 fNonHFE->SetPIDresponse(pidResponse);
636 fNonHFE->SetTrackCuts(-3,3);
637
638 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
639 fNonHFE->SetHistAngle(fOpeningAngleULS);
640 //fNonHFE->SetHistDCABack(fDCABack);
641 //fNonHFE->SetHistDCA(fDCA);
642 fNonHFE->SetHistMassBack(fInvmassLS1);
643 fNonHFE->SetHistMass(fInvmassULS1);
644
645 fNonHFE->FindNonHFE(iTracks,track,fAOD);
646
647 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
648 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
649 // Bool_t fUlsIsPartner = kFALSE;
650 // Bool_t fLsIsPartner = kFALSE;
651 if(fNonHFE->IsULS()){
652 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
653 fULSElecPt->Fill(track->Pt());
654 }
655 }
656
657 if(fNonHFE->IsLS()){
658 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
659 fLSElecPt->Fill(track->Pt());
660 }
661 }
662 }
663
664 if(!fDCA){
665 //----------------------Selection of Photonic Electrons KFParticle-----------------------------
666 Bool_t fFlagPhotonicElec = kFALSE;
667 SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
668 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
669 // Semi inclusive electron
670 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
671 }
672
673
674
675 }//end loop on track
676
677 PostData(1, fOutputList);
678 PostData(2, fFlowEvent);
9d63a2f6 679 if(fSideBandsFlow){
96387d94 680 PostData(3, fFlowEventCont);
9d63a2f6 681 }
96387d94 682
683 //----------hfe end---------
914042c2 684}
685//_________________________________________
a70c9e97 686void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
914042c2 687{
96387d94 688 //Identify non-heavy flavour electrons using Invariant mass method KF
689
690 Bool_t flagPhotonicElec = kFALSE;
691
692 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
693 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
694 if (!trackAsso) {
695 printf("ERROR: Could not receive track %d\n", jTracks);
696 continue;
697 }
698 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
699 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
700 // if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit) || (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
701
702 if(fAssoITSRefit){
703 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
704 }
705
706 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
707
708 if(jTracks == itrack) continue;
709 Double_t ptAsso=-999., nsigma=-999.0;
710 Double_t mass=-999., width = -999;
711 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
712 Double_t openingAngle = -999.;
713
714 ptAsso = trackAsso->Pt();
715 Short_t chargeAsso = trackAsso->Charge();
716 Short_t charge = track->Charge();
717 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
718
719 //80
720 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
721 if(nsigma < -3 || nsigma > 3) continue;
722 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
723 if(ptAsso <0.3) continue;
724
725 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
726 if(charge>0) fPDGe1 = -11;
727 if(chargeAsso>0) fPDGe2 = -11;
728
729 if(charge == chargeAsso) fFlagLS = kTRUE;
730 if(charge != chargeAsso) fFlagULS = kTRUE;
731
732 AliKFParticle::SetField(fAOD->GetMagneticField());
733 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
734 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
735 AliKFParticle recg(ge1, ge2);
736
737 if(recg.GetNDF()<1) continue;
738 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
739 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
740 recg.GetMass(mass,width);
741
742 openingAngle = ge1.GetAngle(ge2);
743 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
744 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
745 if(fOP_angle)if(openingAngle > fOpeningAngleCut) continue;
746
747
748 if(fFlagLS) fInvmassLS1->Fill(mass);
749 if(fFlagULS) fInvmassULS1->Fill(mass);
750
751 if(mass<fInvmassCut){
752 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
753 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
754 }
755
756 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
757 flagPhotonicElec = kTRUE;
758 }
759 }//track loop
760 fFlagPhotonicElec = flagPhotonicElec;
914042c2 761}
8e8d9d2a 762//__________________________________________________________________________________
a70c9e97 763void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
914042c2 764{
96387d94 765 //Create histograms
766
767 //----------hfe initialising begin---------
768 fNullCuts = new AliFlowTrackCuts("null_cuts");
769
770 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
771 cc->SetNbinsMult(10000);
772 cc->SetMultMin(0);
773 cc->SetMultMax(10000);
774
775 cc->SetNbinsPt(100);
776 cc->SetPtMin(0);
777 cc->SetPtMax(50);
778
779 cc->SetNbinsPhi(180);
780 cc->SetPhiMin(0.0);
781 cc->SetPhiMax(TMath::TwoPi());
782
783 cc->SetNbinsEta(30);
784 cc->SetEtaMin(-7.0);
785 cc->SetEtaMax(+7.0);
786
787 cc->SetNbinsQ(500);
788 cc->SetQMin(0.0);
789 cc->SetQMax(3.0);
790
791 //--------Initialize PID
792 fPID->SetHasMCData(kFALSE);
793 if(!fPID->GetNumberOfPIDdetectors())
794 {
795 fPID->AddDetector("TPC", 0);
796 fPID->AddDetector("EMCAL", 1);
797 }
798
799 fPID->SortDetectors();
800 fPIDqa = new AliHFEpidQAmanager();
801 fPIDqa->Initialize(fPID);
802
803 //--------Initialize correction Framework and Cuts
804 fCFM = new AliCFManager;
805 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
806 fCFM->SetNStepParticle(kNcutSteps);
807 for(Int_t istep = 0; istep < kNcutSteps; istep++)
808 fCFM->SetParticleCutsList(istep, NULL);
809
810 if(!fCuts){
811 AliWarning("Cuts not available. Default cuts will be used");
812 fCuts = new AliHFEcuts;
813 fCuts->CreateStandardCuts();
814 }
815
816 fCuts->SetAOD();
817 fCuts->Initialize(fCFM);
818 //----------hfe initialising end--------
819 //---------Output Tlist
820 fOutputList = new TList();
821 fOutputList->SetOwner();
822 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
823
824 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
825 fOutputList->Add(fNoEvents);
826
827 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
828 fOutputList->Add(fTPCnsigma);
829
830 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
831 fOutputList->Add(fTPCnsigmaAft);
832
833 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
834 fOutputList->Add(fTrkEovPBef);
835
836 // fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
837 // fOutputList->Add(fdEdxBef);
838
839 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
840 fOutputList->Add(fInclusiveElecPt);
841
842 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
843 fOutputList->Add(fPhotoElecPt);
ff36a82c 844
96387d94 845 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
846 fOutputList->Add(fSemiInclElecPt);
847
848 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
849 fOutputList->Add(fULSElecPt);
850
914042c2 851 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
852 fOutputList->Add(fLSElecPt);
96387d94 853
854 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
855 fOutputList->Add(fInvmassLS1);
856
857 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
858 fOutputList->Add(fInvmassULS1);
859
860 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
861 fOutputList->Add(fCentralityPass);
862
863 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
864 fOutputList->Add(fCentralityNoPass);
865
ff36a82c 866 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
867 fOutputList->Add(fCentralityBeforePileup);
868
869 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
870 fOutputList->Add(fCentralityAfterVZTRK);
871
96387d94 872 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
873 fOutputList->Add(fPhi);
874
875 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
876 fOutputList->Add(fEta);
877
878 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
879 fOutputList->Add(fVZEROA);
880
881 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
882 fOutputList->Add(fVZEROC);
883
884 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
885 fOutputList->Add(fTPCM);
886
887 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
888 fOutputList->Add(fvertex);
889
890 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
891 fOutputList->Add(fMultCorBeforeCuts);
892
893 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
894 fOutputList->Add(fMultCorAfterCuts);
895
ff36a82c 896 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
897 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
898
899 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
900 fOutputList->Add(fMultCorAfterVZTRKComp);
901
96387d94 902 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
903 fOutputList->Add(fMultvsCentr);
9d63a2f6 904
20791315 905 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
906 fOutputList->Add(fOpeningAngleLS);
907
908 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
909 fOutputList->Add(fOpeningAngleULS);
910
96387d94 911
912 //----------------------------------------------------------------------------
9d63a2f6 913 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
914 fOutputList->Add(EPVzA);
915 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
916 fOutputList->Add(EPVzC);
917 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
918 fOutputList->Add(EPTPC);
96387d94 919 //----------------------------------------------------------------------------
9d63a2f6 920 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
921 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
922 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
923 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
924 fOutputList->Add(fSubEventDPhiv2);
96387d94 925 //================================Event Plane with VZERO=====================
9d63a2f6 926 const Int_t nPtBins = 10;
927 Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
928 // v2A, v2C, pt
929 Int_t bins[3] = { 50, 50, nPtBins};
930 Double_t xmin[3] = { -1., -1., 0};
931 Double_t xmax[3] = { 1., 1., 8};
932 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
933 // Set bin limits for axes which are not standard binned
934 fV2Phi->SetBinEdges(2, binsPt);
935 // set axes titles
936 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
937 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
938 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
939 fOutputList->Add(fV2Phi);
96387d94 940 //----------------------------------------------------------------------------
9d63a2f6 941 if(fPhiminusPsi){
96387d94 942 Int_t binsvElectH[8]={ 600, 200, 200 ,100, 100, 100, 10, 10}; //pt, E/p,TPCnSigma,M20,M02,Disp Phi-psiV0A ,Phi-PsiV0C,eta (commented)
943 Double_t xminvElectH[8]={0, 0, -10 , 0, 0, 0, 0, 0};
944 Double_t xmaxvElectH[8]={20, 2, 10 , 2, 2, 2, TMath::Pi(), TMath::Pi()};
945 fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
946 fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
947 fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
948 fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
949 fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
950 fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
951 fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
952 fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
953 fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
954 fOutputList->Add(fSparseElectronHadron);
9d63a2f6 955 }
96387d94 956 //----------------------------------------------------------------------------
957 //----------------------------------------------------------------------------
9d63a2f6 958 if(fpurity){
96387d94 959 Int_t binsvpurity[3]={ 600,200, 200}; //pt, E/p,TPCnSigma
960 Double_t xminvpurity[3]={0, 0, -10};
ff36a82c 961 Double_t xmaxvpurity[3]={30, 2, 10};
96387d94 962 fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
963 fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
964 fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
965 fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
966 fOutputList->Add(fSparseElectronpurity);
9d63a2f6 967 }
96387d94 968 //----------------------------------------------------------------------------
969
970 PostData(1,fOutputList);
971 // create and post flowevent
972 fFlowEvent = new AliFlowEvent(10000);
973 PostData(2, fFlowEvent);
9d63a2f6 974
975 if(fSideBandsFlow){
96387d94 976 fFlowEventCont = new AliFlowEvent(10000);
977 PostData(3, fFlowEventCont);
9d63a2f6 978 }
96387d94 979}
914042c2 980
981//________________________________________________________________________
a70c9e97 982void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
914042c2 983{
96387d94 984 // Info("Terminate");
985 AliAnalysisTaskSE::Terminate();
914042c2 986}
987//_____________________________________________________________________________
a70c9e97 988template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
914042c2 989{
96387d94 990 // QA multiplicity plots
991 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
992 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
914042c2 993}
994//_____________________________________________________________________________
a70c9e97 995template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
914042c2 996{
96387d94 997 //Set null cuts
914042c2 998 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
999 fCutsRP->SetEvent(event, MCEvent());
1000 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1001 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1002 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1003 fNullCuts->SetEvent(event, MCEvent());
1004}
1005//_____________________________________________________________________________
96387d94 1006void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
914042c2 1007{
96387d94 1008 //Prepare flow events
1009 FlowEv->ClearFast();
1010 FlowEv->Fill(fCutsRP, fNullCuts);
1011 FlowEv->SetReferenceMultiplicity(iMulti);
1012 FlowEv->DefineDeadZone(0, 0, 0, 0);
1013 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
914042c2 1014}
1015//_____________________________________________________________________________
a70c9e97 1016Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
914042c2 1017{
96387d94 1018 // Check single track cuts for a given cut step
1019 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1020 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1021 return kTRUE;
914042c2 1022}
1023//_________________________________________
a70c9e97 1024void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
914042c2 1025{
ff36a82c 1026//============================Multiplicity TPV vs Global===============================================================================
9d63a2f6 1027 const Int_t nGoodTracks = event->GetNumberOfTracks();
9d63a2f6 1028 Float_t multTPC(0.); // tpc mult estimate
1029 Float_t multGlob(0.); // global multiplicity
1030 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1031 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1032 if (!trackAOD) continue;
1033 if (!(trackAOD->TestFilterBit(1))) continue;
1034 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;
1035 multTPC++;
1036 }
1037 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1038 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1039 if (!trackAOD) continue;
1040 if (!(trackAOD->TestFilterBit(16))) continue;
1041 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;
1042 Double_t b[2] = {-99., -99.};
1043 Double_t bCov[3] = {-99., -99., -99.};
1044 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1045 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1046 multGlob++;
ff36a82c 1047 } //track loop
1048 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1049//============================================================================================================================
1050 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1051 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1052 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1053 // cout << "--------------Centrality evaluated-------------------------"<<endl;
1054 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1055 {
1056 fCentralityNoPass->Fill(fCentrality);
1057 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1058 centralitypass = kFALSE;
1059 }else
1060 {
1061 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1062 centralitypass = kTRUE;
1063 }
1064 if (centralitypass){
1065 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1066 fCentralityBeforePileup->Fill(fCentrality);
1067 }//...after centrality selectrion
1068//============================================================================================================================
1069//to remove the bias introduced by multeplicity outliers---------------------
1070 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1071 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1072 if (TMath::Abs(centv0 - centTrk) > 5.0){
1073 centralitypass = kFALSE;
1074 fCentralityNoPass->Fill(fCentrality);
1075 }
1076 if (centralitypass){
1077 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1078 fCentralityAfterVZTRK->Fill(fCentrality);
1079 }//...after centrality selectrion
1080//============================================================================================================================
1081 if(fMultCut){
1082 if(fTrigger==1){
96387d94 1083 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
edaa818d 1084 // cout <<" Trigger ==" <<fTrigger<< endl;
a70c9e97 1085 centralitypass = kFALSE;
1086 fCentralityNoPass->Fill(fCentrality);
edaa818d 1087 }//2011 Semicentral
1088 }
1089 if(fTrigger==0){
edaa818d 1090 if(! (multTPC > (77.6 + 1.399*multGlob) && multTPC < (371.3 + 1.50*multGlob))){
1091 // cout <<" Trigger ==" <<fTrigger<< endl;
edaa818d 1092 centralitypass = kFALSE;
1093 fCentralityNoPass->Fill(fCentrality);
1094 }//2011
1095 }//2011 Central
ff36a82c 1096 }
1097//=================================All cuts are passed==================++++==================================================
9d63a2f6 1098 if(centralitypass){
96387d94 1099 fCentralityPass->Fill(fCentrality);
1100 fMultCorAfterCuts->Fill(multGlob, multTPC);
1101 fMultvsCentr->Fill(fCentrality, multTPC);
a70c9e97 1102 }
914042c2 1103}
1104//_____________________________________________________________________________
a70c9e97 1105void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
914042c2 1106{
96387d94 1107 // Set a centrality range ]min, max] and define the method to use for centrality selection
1108 fCentralityMin = CentralityMin;
1109 fCentralityMax = CentralityMax;
1110 fkCentralityMethod = CentralityMethod;
914042c2 1111}
1112//_____________________________________________________________________________
a70c9e97 1113void AliAnalysisTaskFlowTPCEMCalQCSP::SetIDCuts(Double_t minTPC, Double_t maxTPC, Double_t minEovP, Double_t maxEovP, Double_t minM20, Double_t maxM20, Double_t minM02, Double_t maxM02, Double_t Dispersion)
914042c2 1114{
1115 //Set ID cuts
1116 fminTPC = minTPC;
1117 fmaxTPC = maxTPC;
1118 fminEovP = minEovP;
1119 fmaxEovP = maxEovP;
1120 fminM20 = minM20;
1121 fmaxM20 = maxM20;
1122 fminM02 = minM02;
1123 fmaxM02 = maxM02;
1124 fDispersion = Dispersion;
1125}
1126//_____________________________________________________________________________