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