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