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