updated for nonHFE reconstruction study
[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 ,fpTCut(0)
127 ,fTrigger(0)
128 ,fPhi(0)
129 ,fEta(0)
130 ,fVZEROA(0)
131 ,fVZEROC(0)
132 ,fTPCM(0)
133 ,fNoEvents(0)
134 ,fInclusiveElecPt(0)
135 //,fTPCnsigma(0)
136 ,fTPCnsigmaAft(0)
137 //,fTOFns(0)
138 ,fTOFnsAft(0)
139 ,fTOFBetaAft(0)
140 ,fCentralityPass(0)
141 ,fCentralityNoPass(0)
142 ,fInvmassLS1(0)
143 ,fInvmassULS1(0)
144 ,fPhotoElecPt(0)
145 ,fSemiInclElecPt(0)
146 ,fULSElecPt(0)
147 ,fLSElecPt(0)
148 ,fMultCorAfterCuts(0)
149 ,fMultvsCentr(0)
150 ,fSubEventDPhiv2(0)
151 ,EPVzA(0)
152 ,EPVzC(0)
153 ,EPTPC(0)
154 ,fV2Phi(0)
155 ,fvertex(0)
156 ,fMultCorBeforeCuts(0)
157 ,fPIDResponse(0)
158 ,fQAPid(0)
159 ,fminTPCnsigma(0)
160 ,fmaxTPCnsigma(3)
161 ,fQAPIDSparse(kFALSE)
162 ,fminTOFnSigma(-2)
163 ,fmaxTOFnSigma(2)
164 ,fQAPidSparse(0)
165 ,fTPCS(0)
166 ,fVz(0)
167 ,fPhiCut(0)
168 ,fOpeningAngleCut(0.1)
169 ,fOP_angle(0)
170 ,fOpeningAngleLS(0)
171 ,fOpeningAngleULS(0)
172 {
173   //Named constructor
174
175 //  fPID = new AliHFEpid("hfePid");
176   // Define input and output slots here
177   // Input slot #0 works with a TChain
178   DefineInput(0, TChain::Class());
179   // Output slot #0 id reserved by the base class for AOD
180   // Output slot #1 writes into a TH1 container
181   // DefineOutput(1, TH1I::Class());
182   DefineOutput(1, TList::Class());
183   DefineOutput(2, AliFlowEventSimple::Class());
184 }
185
186 //________________________________________________________________________
187 AliAnalysisTaskFlowTPCTOFQCSP::AliAnalysisTaskFlowTPCTOFQCSP() 
188   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElectFlow")
189 ,fDebug(0)
190 ,fAOD(0)
191 ,fOutputList(0)
192 ,fCuts(0)
193 ,fIdentifiedAsOutInz(kFALSE)
194 ,fPassTheEventCut(kFALSE)
195 ,fCFM(0)
196 //,fPID(0)
197 //,fPIDqa(0)
198 ,fCutsRP(0)     // track cuts for reference particles
199 ,fNullCuts(0) // dummy cuts for flow event tracks
200 ,fFlowEvent(0) //! flow events (one for each inv mass band)
201 ,fkCentralityMethod(0)
202 ,fCentrality(0)
203 ,fCentralityMin(0)
204 ,fCentralityMax(0)
205 ,fInvmassCut(0)
206 ,fpTCut(0)
207 ,fTrigger(0)
208 ,fPhi(0)
209 ,fEta(0)
210 ,fVZEROA(0)
211 ,fVZEROC(0)
212 ,fTPCM(0)
213 ,fNoEvents(0)
214 ,fInclusiveElecPt(0)
215 //,fTPCnsigma(0)
216 ,fTPCnsigmaAft(0)
217 //,fTOFns(0)
218 ,fTOFnsAft(0)
219 ,fTOFBetaAft(0)
220 ,fCentralityPass(0)
221 ,fCentralityNoPass(0)
222 ,fInvmassLS1(0)
223 ,fInvmassULS1(0)
224 ,fPhotoElecPt(0)
225 ,fSemiInclElecPt(0)
226 ,fULSElecPt(0)
227 ,fLSElecPt(0)
228 ,fMultCorAfterCuts(0)
229 ,fMultvsCentr(0)
230 ,fSubEventDPhiv2(0)
231 ,EPVzA(0)
232 ,EPVzC(0)
233 ,EPTPC(0)
234 ,fV2Phi(0)
235 ,fvertex(0)
236 ,fMultCorBeforeCuts(0)
237 ,fPIDResponse(0)
238 ,fQAPid(0)
239 ,fminTPCnsigma(0)
240 ,fmaxTPCnsigma(3)
241 ,fQAPIDSparse(kFALSE)
242 ,fminTOFnSigma(-2)
243 ,fmaxTOFnSigma(2)
244 ,fQAPidSparse(0)
245 ,fTPCS(0)
246 ,fVz(0)
247 ,fPhiCut(0)
248 ,fOpeningAngleCut(0.1)
249 ,fOP_angle(0)
250 ,fOpeningAngleLS(0)
251 ,fOpeningAngleULS(0)
252 {
253   //Default constructor
254 //  fPID = new AliHFEpid("hfePid");
255   // Constructor
256   // Define input and output slots here
257   // Input slot #0 works with a TChain
258   DefineInput(0, TChain::Class());
259   // Output slot #0 id reserved by the base class for AOD
260   // Output slot #1 writes into a TH1 container
261   // DefineOutput(1, TH1I::Class());
262   DefineOutput(1, TList::Class());
263   DefineOutput(2, AliFlowEventSimple::Class());
264     //DefineOutput(3, TTree::Class());
265 }
266 //_________________________________________
267
268 AliAnalysisTaskFlowTPCTOFQCSP::~AliAnalysisTaskFlowTPCTOFQCSP()
269 {
270   //Destructor 
271
272   delete fOutputList;
273 //  delete fPID;
274   delete fPIDResponse;
275   delete fCFM;
276 //  delete fPIDqa;
277   if (fOutputList) delete fOutputList;
278   if (fFlowEvent) delete fFlowEvent;
279
280 }
281 //_________________________________________
282
283 void AliAnalysisTaskFlowTPCTOFQCSP::UserExec(Option_t*)
284 {
285   //Main loop
286   //Called for each event
287
288   // create pointer to event
289
290   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
291   if (!fAOD)
292   {
293     printf("ERROR: fAOD not available\n");
294     return;
295   }
296
297   if(!fCuts)
298   {
299     AliError("HFE cuts not available");
300     return;
301   }
302
303 //  if(!fPID->IsInitialized())
304 //  {
305     // Initialize PID with the given run number
306 //    AliWarning("PID not initialised, get from Run no");
307 //    fPID->InitializePID(fAOD->GetRunNumber());
308 //  }
309     
310   // cout << "kTrigger   ==   " << fTrigger <<endl;
311     
312     if(fTrigger==0){
313 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
314     }
315     if(fTrigger==1){
316 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
317     }
318     if(fTrigger==2){
319 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
320     }
321     if(fTrigger==3){
322 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
323     }
324   
325     
326 //---------------CENTRALITY AND EVENT SELECTION-----------------------
327     
328     
329     
330    Int_t fNOtrks =  fAOD->GetNumberOfTracks();
331     Float_t vtxz = -999;
332     const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
333     if (!trkVtx || trkVtx->GetNContributors()<=0)return;
334     TString vtxTtl = trkVtx->GetTitle();
335     if (!vtxTtl.Contains("VertexerTracks"))return;
336     const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
337     if (!spdVtx || spdVtx->GetNContributors()<=0)return;
338     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
339     vtxz = trkVtx->GetZ();
340     if(TMath::Abs(vtxz)>fVz)return;
341
342 // Event cut
343     if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
344     if(fNOtrks<2) return;
345     
346
347     Bool_t pass = kFALSE; //to select centrality
348     CheckCentrality(fAOD,pass);
349     if(!pass)return;
350     
351     fvertex->Fill(vtxz);
352
353     
354     fNoEvents->Fill(0);
355     PlotVZeroMultiplcities(fAOD);
356
357     SetNullCuts(fAOD);
358     PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
359     
360  
361     fCFM->SetRecEventInfo(fAOD);
362
363   // Look for kink mother
364   Int_t numberofvertices = fAOD->GetNumberOfVertices();
365   Double_t listofmotherkink[numberofvertices];
366   Int_t numberofmotherkink = 0;
367   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
368     AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
369     if(!aodvertex) continue;
370     if(aodvertex->GetType()==AliAODVertex::kKink) {
371       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
372       if(!mother) continue;
373       Int_t idmother = mother->GetID();
374       listofmotherkink[numberofmotherkink] = idmother;
375       //printf("ID %d\n",idmother);
376       numberofmotherkink++;
377     }
378   }
379     
380 //=============================================V0EP from Alex======================================================================
381     Double_t qxEPa = 0, qyEPa = 0;
382     Double_t qxEPc = 0, qyEPc = 0;
383     
384     Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
385     Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
386     
387     
388     Double_t Qx2 = 0, Qy2 = 0;
389     
390     for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
391         
392         AliAODTrack* aodTrack = fAOD->GetTrack(iT);
393         
394         if (!aodTrack)
395             continue;
396         
397         if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
398             continue;
399         
400         if (!aodTrack->TestFilterBit(128))
401             continue;
402         
403         Qx2 += TMath::Cos(2*aodTrack->Phi());
404         Qy2 += TMath::Sin(2*aodTrack->Phi());
405     }
406     
407     Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
408     
409     EPVzA->Fill(evPlAngV0A);
410     EPVzC->Fill(evPlAngV0C);
411     EPTPC->Fill(evPlAngTPC);
412     
413     fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
414     fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
415     fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
416 //====================================================================================================================
417     
418  AliAODTrack *track = NULL;
419
420 // Track loop 
421  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) 
422  {
423    track = fAOD->GetTrack(iTracks);
424    if (!track) 
425    {
426      printf("ERROR: Could not receive track %d\n", iTracks);
427      continue;
428    }
429      
430    if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
431      
432 //--------------------------------------hfe begin-----------------------------------------------------------
433 //==========================================================================================================
434 //======================================track cuts==========================================================
435    if(track->Eta()<-0.8 || track->Eta()>0.8)    continue;    //eta cuts on candidates
436
437    // RecKine: ITSTPC cuts  
438    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
439
440    // Reject kink mother
441    Bool_t kinkmotherpass = kTRUE;
442    for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
443      if(track->GetID() == listofmotherkink[kinkmother]) {
444        kinkmotherpass = kFALSE;
445        continue;
446      }
447    }
448    if(!kinkmotherpass) continue;
449
450    // RecPrim
451  //  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;  //deleted for DCA absence
452    // HFEcuts: ITS layers cuts
453    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
454    // HFE cuts: TPC PID cleanup
455    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
456 //==========================================================================================================
457      Double_t eta = track->Eta();
458      Double_t phi = track->Phi();
459      Double_t pt = track->Pt();         //pt track after cuts
460      if(pt<fpTCut) continue;
461 //==========================================================================================================
462 //=========================================PID==============================================================
463      if(track->GetTPCsignalN() < fTPCS) continue;
464      Float_t fTPCnSigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
465      Float_t fTOFnSigma = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
466      
467   //   Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
468
469 //     fTPCnsigma->Fill(track->P(),fTPCnSigma);
470 //     fTOFns->Fill(track->P(),fTOFnSigma);
471      
472      Double_t valPidSparse[4] = {
473          pt,
474          fTPCnSigma,
475          fTOFnSigma,
476          phi
477         };
478      fQAPidSparse->Fill(valPidSparse);
479
480      /* || track->GetTPCsignal() < 72 || track->GetTPCsignal() > 95*/
481
482      if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
483      if(fTPCnSigma < fminTPCnsigma || fTPCnSigma > fmaxTPCnsigma) continue; //cuts on nsigma tpc
484 //==========================================================================================================
485 //=========================================QA PID SPARSE====================================================
486      Float_t timeTOF = track->GetTOFsignal();
487      Double_t intTime[5] = {-99., -99., -99., -99., -99.};
488      track->GetIntegratedTimes(intTime);
489      Float_t timeElec = intTime[0];
490      Float_t intLength = 2.99792458e-2* timeElec;
491      Double_t beta = 0.1;
492      if ((intLength > 0) && (timeTOF > 0))
493          beta = intLength/2.99792458e-2/timeTOF;
494    
495     if(fQAPIDSparse){
496      Double_t valPid[4] = {
497          track->P(),
498          track->GetTPCsignal(),
499          beta,
500          track->Charge()
501          };
502      fQAPid->Fill(valPid);
503      }
504    
505      
506       if(fPhiCut){
507       if(phi<1.4 || phi >3.14)continue; //to have same EMCal phi acceptance
508       }
509      
510 //==========================================================================================================
511 //============================Event Plane Method with V0====================================================
512      Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
513      Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
514      Double_t v2Phi[3] = {
515          v2PhiV0A,
516          v2PhiV0C,
517          pt};
518      fV2Phi->Fill(v2Phi);
519 //=========================================================================================================
520    fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
521    fTOFnsAft->Fill(track->P(),fTOFnSigma);
522    fTOFBetaAft->Fill(track->P(),beta);
523    fInclusiveElecPt->Fill(pt);
524    fPhi->Fill(phi);
525    fEta->Fill(eta);
526 //=========================================================================================================
527 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
528    AliFlowTrack *sTrack = new AliFlowTrack();
529      sTrack->Set(track);
530      sTrack->SetID(track->GetID());
531      sTrack->SetForRPSelection(kFALSE);
532      sTrack->SetForPOISelection(kTRUE);
533      sTrack->SetMass(263732);
534      for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
535      {
536       //   cout << " no of rps " << iRPs << endl;
537          AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
538          if (!iRP) continue;
539          if (!iRP->InRPSelection()) continue;
540         if( sTrack->GetID() == iRP->GetID())
541         {
542           if(fDebug) printf(" was in RP set");
543           //  cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
544             iRP->SetForRPSelection(kFALSE);
545             fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
546         }
547       } //end of for loop on RPs
548    fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
549 //=========================================================================================================
550 //----------------------Selection and Flow of Photonic Electrons-----------------------------
551    Bool_t fFlagPhotonicElec = kFALSE;
552    SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
553    if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
554    // Semi inclusive electron 
555    if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
556 //-------------------------------------------------------------------------------------------
557
558  }//end loop on track
559  PostData(1, fOutputList);
560  PostData(2, fFlowEvent);
561  //----------hfe end---------
562 }
563 //_________________________________________
564 void AliAnalysisTaskFlowTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
565 {
566     
567   //Identify non-heavy flavour electrons using Invariant mass method
568   Bool_t flagPhotonicElec = kFALSE;
569
570   for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
571     AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
572     if (!trackAsso) {
573       printf("ERROR: Could not receive track %d\n", jTracks);
574       continue;
575     }
576     //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
577       if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
578       if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
579
580       
581     if(jTracks == itrack) continue;
582     Double_t ptAsso=-999., nsigma=-999.0;
583     Double_t mass=-999., width = -999;
584       Double_t openingAngle = -999.;
585     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
586
587   
588     ptAsso = trackAsso->Pt();
589     Short_t chargeAsso = trackAsso->Charge();
590     Short_t charge = track->Charge();
591    // nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
592       nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
593
594     if(trackAsso->GetTPCNcls() < 80) continue;
595     if(nsigma < -3 || nsigma > 3) continue;
596     if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
597     if(ptAsso <0.3) continue;
598
599     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
600     if(charge>0) fPDGe1 = -11;
601     if(chargeAsso>0) fPDGe2 = -11;
602
603     if(charge == chargeAsso) fFlagLS = kTRUE;
604     if(charge != chargeAsso) fFlagULS = kTRUE;
605
606     AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
607     AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
608     AliKFParticle recg(ge1, ge2);
609
610     if(recg.GetNDF()<1) continue;
611     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
612     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
613     recg.GetMass(mass,width);
614
615       if(fOP_angle){
616     openingAngle = ge1.GetAngle(ge2);
617     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
618     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
619     if(openingAngle > fOpeningAngleCut) continue;
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