o updates for PbPb analysis of B->J/psi (Fiorella)
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliAnalysisTaskDielectronReadAODBranch.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 //  Class to retrieve the branch of the dielectron candidates stored in   //
18 //  filtered AODs file (AliAOD.Dielectron.root). It is possible to        //
19 //  apply tighter cuts to candidates stored in the branch and do the      //
20 //  matching with MC truth.                                               //
21 ////////////////////////////////////////////////////////////////////////////
22
23 #include "TChain.h"
24 #include "TNtuple.h"
25 #include "TList.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TDatabasePDG.h"
29 #include "TF1.h"
30 #include "AliAnalysisManager.h"
31 #include "AliAODHandler.h"
32 #include "AliAODEvent.h"
33 #include "AliAODVertex.h"
34 #include "AliAODTrack.h"
35 #include "AliAODMCParticle.h"
36 #include "AliAnalysisTaskDielectronReadAODBranch.h"
37 #include "AliDielectronPair.h"
38 #include "AliDielectronMC.h"
39
40 ClassImp(AliAnalysisTaskDielectronReadAODBranch)
41         //________________________________________________________________________
42         AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch():
43                 AliAnalysisTaskSE(),
44                 fOutput(0),
45                 fNtupleJPSI(0),
46                 fNentries(0),
47                 fInvMass(0),
48                 fInvMassNoCuts(0),
49                 fpsproperSignal(0),
50                 fpsproperSidebands(0),            
51                 fpsproperAll(0),                  
52                 fpsproperUnder(0),
53                 fpsproperUpper(0),
54                 fLxyVsPtleg1(0),
55                 fLxyVsPtleg2(0),
56                 fLxyVsPt(0),
57                 fMeeVsPt(0),
58                 fMeeVsLxy(0),
59                 fprimvtxZ(0),  
60                 fsecvtxZ(0),  
61                 fprimvtxX(0),  
62                 fsecvtxX(0),  
63                 fprimvtxY(0),  
64                 fsecvtxY(0),  
65                 fPt(0),
66                 fPtLeg1(0),
67                 fPtLeg2(0),
68                 fdEdxP(0),
69                 fHasMC(0),
70                 fobj(0),
71                 fobjMC(0),
72                 fPtCut(1.),
73                 fSpdFirstRequired(kFALSE),
74                 fClsTPC(90),
75                 fPairType(1),
76                 fPtJpsi(1.3),
77                 fInvMassSignalLimits(0),
78                 fInvMassSideBandsLimits(0),
79                 fSecondary(0)
80 {
81         // Default constructor
82 }
83
84 //________________________________________________________________________
85 AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(const char *name):
86         AliAnalysisTaskSE(name),
87         fOutput(0),
88         fNtupleJPSI(0),
89         fNentries(0),
90         fInvMass(0),
91         fInvMassNoCuts(0),
92         fpsproperSignal(0),
93         fpsproperSidebands(0),            
94         fpsproperAll(0),                  
95         fpsproperUnder(0),
96         fpsproperUpper(0),
97         fLxyVsPtleg1(0),
98         fLxyVsPtleg2(0),
99         fLxyVsPt(0),
100         fMeeVsPt(0),
101         fMeeVsLxy(0),
102         fprimvtxZ(0), 
103         fsecvtxZ(0),  
104         fprimvtxX(0), 
105         fsecvtxX(0),  
106         fprimvtxY(0), 
107         fsecvtxY(0),  
108         fPt(0),
109         fPtLeg1(0),
110         fPtLeg2(0),
111         fdEdxP(0),
112         fHasMC(0),
113         fobj(0),
114         fobjMC(0),
115         fPtCut(1.),
116         fSpdFirstRequired(kFALSE),
117         fClsTPC(90),
118         fPairType(1),
119         fPtJpsi(1.3),
120         fInvMassSignalLimits(0),
121         fInvMassSideBandsLimits(0), 
122         fSecondary(0)
123 {
124         // Default constructor
125         DefineInput(0,TChain::Class());
126         DefineOutput(1,TList::Class());  //My private output
127
128         fInvMassSignalLimits = new Double_t[2]; 
129         fInvMassSignalLimits[0] = 0.;  fInvMassSignalLimits[1] = 0.; 
130         fInvMassSideBandsLimits = new Double_t[2]; 
131         fInvMassSideBandsLimits[0] = 0.;  fInvMassSideBandsLimits[1] = 0.;
132 }
133
134 //________________________________________________________________________
135 AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
136 {
137         if(fOutput)   delete fOutput;
138         if(fobj)      delete fobj;
139         if(fobjMC)    delete fobjMC;
140 }
141
142 //________________________________________________________________________
143 void AliAnalysisTaskDielectronReadAODBranch::Init()
144 {
145         // Initialization
146         if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
147         return;
148 }
149
150 //________________________________________________________________________
151 void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
152 {
153         //
154         // Create the output container
155         //
156         if(fDebug > 1) printf("AnalysisTaskReadAOD::UserCreateOutputObjects() \n");
157         fOutput = new TList();
158         fOutput->SetOwner();
159         //// Histogram booking
160         
161         // invariant mass
162         fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
163         fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
164         fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);
165
166         //pseudoproper 
167         fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
168         fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f );  X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
169         fpsproperAll = new TH1F("psproper_decay_length_all",Form("psproper_decay_length_distrib_all(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
170         fpsproperUnder = new TH1F("psproper_decay_length_under",Form("psproper_decay_length_distrib_under(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.); 
171         fpsproperUpper = new TH1F("psproper_decay_length_upper",Form("psproper_decay_length_distrib_upper(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);  
172         //
173         fLxyVsPtleg1 = new TH2F("Lxy_vs_Pt_leg1",Form("Lxy_vs_Pt_leg1_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
174         fLxyVsPtleg2 = new TH2F("Lxy_vs_Pt_leg2",Form("Lxy_vs_Pt_leg2_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
175         fLxyVsPt = new TH2F("Lxy_vs_Pt_jpsi",Form("Lxy_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
176         fMeeVsPt = new TH2F("Mee_vs_Pt_jpsi",Form("Mee_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; M_{ee}",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,200, 1.99,4.1);
177         fMeeVsLxy = new TH2F("Mee_vs_Lxy_jpsi",Form("Mee_vs_Lxy_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV),pT(J/#psi)>%f GeV/c; X; M_{ee}",fClsTPC,fPtCut,fPtJpsi), 400, -3000, 3000.,200, 1.99,4.1);
178
179         // QA plots
180         fprimvtxZ = new TH1F("prim_vtx_Z",Form("prim_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
181         fsecvtxZ = new TH1F("sec_vtx_Z",Form("sec_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f);  Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
182         fprimvtxX = new TH1F("prim_vtx_X",Form("prim_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
183         fsecvtxX = new TH1F("sec_vtx_X",Form("sec_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  X[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
184         fprimvtxY = new TH1F("prim_vtx_Y",Form("prim_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
185         fsecvtxY = new TH1F("sec_vtx_Y",Form("sec_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f);  Y[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
186         //
187         fPt = new TH1F("Pt(J/psi)",Form("Pt_Jpsi_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); p_{T}(J/#psi)[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
188         fPtLeg1 = new TH1F("Pt_leg1",Form("Pt_leg1_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
189         fPtLeg2 = new TH1F("Pt_leg2",Form("Pt_leg2_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f);  p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
190         //dE/dx plots
191         fdEdxP = new TH2F("dE/dx_vs_Ptpc",Form("dE/dx_vs_Ptpc(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f)",fClsTPC,fPtCut,fPtJpsi),400,0.2,20.,200,0.,200.);
192         
193         //J/psi n-tuple (X and M values)
194         fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
195  
196         // add histograms to list
197         fOutput->Add(fNentries);
198         fOutput->Add(fInvMass);
199         fOutput->Add(fInvMassNoCuts); 
200         fOutput->Add(fpsproperSignal);
201         fOutput->Add(fpsproperSidebands);
202         fOutput->Add(fpsproperAll);
203         fOutput->Add(fpsproperUnder);
204         fOutput->Add(fpsproperUpper);
205         fOutput->Add(fLxyVsPtleg1);
206         fOutput->Add(fLxyVsPtleg2);
207         fOutput->Add(fLxyVsPt);
208         fOutput->Add(fMeeVsPt);
209         fOutput->Add(fMeeVsLxy);
210         fOutput->Add(fprimvtxZ);
211         fOutput->Add(fsecvtxZ);
212         fOutput->Add(fprimvtxX);
213         fOutput->Add(fsecvtxX);
214         fOutput->Add(fprimvtxY);
215         fOutput->Add(fsecvtxY);
216         fOutput->Add(fPt);
217         fOutput->Add(fPtLeg1);
218         fOutput->Add(fPtLeg2);
219         fOutput->Add(fNtupleJPSI);
220         fOutput->Add(fdEdxP);
221         return;
222 }
223
224 //________________________________________________________________________
225 void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
226 {
227         // Execute analysis for current event:
228         AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
229         if (!aod) return;
230
231         Double_t vtxPrim[3] = {0.,0.,0.};
232
233         AliAODVertex* primvtx = aod->GetPrimaryVertex();
234         vtxPrim[0] = primvtx->GetX();
235         vtxPrim[1] = primvtx->GetY();
236         vtxPrim[2] = primvtx->GetZ();
237
238         AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
239
240         TTree *aodTree = aodHandler->GetTree();
241         if(!fobj) aodTree->SetBranchAddress("dielectrons",&fobj);
242
243         if(fHasMC){ aodTree->SetBranchAddress("mcparticles",&fobjMC);
244         if(!fobjMC) printf("AliAnalysisTaskDielectronReadAODBranch::UserExec: MC particles branch not found!\n"); }
245
246         fNentries->Fill(0);
247         aodTree->GetEvent(Entry());
248
249         // loop over candidates
250         if(fobj) {
251                 TObjArray *objArr = (TObjArray*)fobj->UncheckedAt(fPairType);
252                 for(int j=0;j<objArr->GetEntriesFast();j++)
253                 {
254                         AliDielectronPair *pairObj = (AliDielectronPair*)objArr->UncheckedAt(j);
255                         fInvMassNoCuts->Fill(pairObj->M()); 
256                         Double_t vtxSec[3] = {0.,0.,0.};
257                         fprimvtxX->Fill(vtxPrim[0]);
258                         fprimvtxY->Fill(vtxPrim[1]);
259                         fprimvtxZ->Fill(vtxPrim[2]);
260                         pairObj->XvYvZv(vtxSec);
261                         fsecvtxX->Fill(vtxSec[0]);
262                         fsecvtxY->Fill(vtxSec[1]);
263                         fsecvtxZ->Fill(vtxSec[2]);
264
265                         Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(pairObj->Px()) + (vtxSec[1]-vtxPrim[1])*(pairObj->Py()))/pairObj->Pt();
266                         Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/pairObj->Pt();
267
268                         AliAODTrack *trk = (AliAODTrack*)pairObj->GetFirstDaughter();
269                         AliAODTrack *trk1 = (AliAODTrack*)pairObj->GetSecondDaughter();
270                         if(!trk || !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}
271
272                         //check in case of MC analysis if candidate is a true J/psi->ee
273                           if(fHasMC && ((AliDielectronMC::Instance()->IsJpsiPrimary(pairObj))!=fSecondary)) continue;
274                            AliAODPid *pid  = trk->GetDetPid();
275                                 AliAODPid *pid1 = trk1->GetDetPid();
276                                 if(!pid || !pid1){
277                                         if(fDebug>1) printf("No AliAODPid found\n");
278                                         continue;
279                                 }
280                                  // ptLeg cut
281                                 if((trk->Pt()<fPtCut) || (trk1->Pt()<fPtCut)) continue;
282                            
283                                 // pt jpsi cut
284                                 if((pairObj->Pt())<fPtJpsi) continue; 
285  
286                               //spd first required
287                               if(fSpdFirstRequired)
288                               {
289                                if((!trk->HasPointOnITSLayer(0)) || (!trk1->HasPointOnITSLayer(0))) continue;
290                               }
291                            //
292                           if((trk->GetTPCNcls()) <= fClsTPC || (trk1->GetTPCNcls()) <= fClsTPC) continue;
293                                 if(trk->Charge()==-1){
294                                         fPtLeg1->Fill(trk->Pt());
295                                 } else {fPtLeg2->Fill(trk->Pt());}
296                                 if(trk1->Charge()==-1){
297                                         fPtLeg1->Fill(trk1->Pt());
298                                 } else {fPtLeg2->Fill(trk1->Pt());}
299                                 //Fill dE/dx related plots (before PID cuts)
300                                 fdEdxP->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
301                                 fdEdxP->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
302                                         fInvMass->Fill(pairObj->M());
303                                         fPt->Fill(pairObj->Pt());
304                                         fMeeVsLxy->Fill(10000*psProperDecayLength,pairObj->M());
305                                         fMeeVsPt->Fill(pairObj->Pt(),pairObj->M());
306                                         fpsproperAll->Fill(10000*psProperDecayLength);
307                                            
308                                            //psproper in the signal region 
309                                            if((pairObj->M())<fInvMassSignalLimits[1] && (pairObj->M())>fInvMassSignalLimits[0]){
310                                                 fpsproperSignal->Fill(10000*psProperDecayLength);
311                                                 fLxyVsPt->Fill(pairObj->Pt(),10000*psProperDecayLength);  //jpsi 
312                                                 fLxyVsPtleg1->Fill(trk->Pt(),10000*psProperDecayLength);  //leg1 (warning: mixture of pos and neg tracks)
313                                                 fLxyVsPtleg2->Fill(trk1->Pt(),10000*psProperDecayLength); //leg2 (warning: mixture of pos and neg tracks)
314                                                 fNtupleJPSI->Fill(10000*psProperDecayLength,pairObj->M()); // fill the N-tuple (imput of the minimization algorithm)
315                                                  } 
316
317                                         // psproper in the sidebands 
318                                         if((pairObj->M())> fInvMassSideBandsLimits[1] || (pairObj->M())< fInvMassSideBandsLimits[0])
319                                         fpsproperSidebands->Fill(10000*psProperDecayLength);
320                                           
321                                         // psproper in the lower sideband
322                                         if(pairObj->M()< fInvMassSideBandsLimits[0]) fpsproperUnder->Fill(10000*psProperDecayLength);
323                                           
324                                         // psproper in the upper sideband
325                                         if(pairObj->M()> fInvMassSideBandsLimits[1]) fpsproperUpper->Fill(10000*psProperDecayLength);
326                         } // end loop over candidates
327                 }// end loop over pair types
328         // Post the data
329         PostData(1,fOutput);
330         return;
331
332
333 //________________________________________________________________________
334 void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
335 {
336         if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
337         return;
338 }
339