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