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