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