Change of cuts for two track loop
[u/mrichter/AliRoot.git] / PWGUD / UPC / AliAnalysisTaskUpcPsi2s.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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 // c++ headers
17 #include <iostream>
18 #include <string.h>
19
20 // root headers
21 #include "TH1I.h"
22 #include "TTree.h"
23 #include "TClonesArray.h"
24 #include "TParticle.h"
25 #include "TObjString.h"
26 #include "TFile.h"
27 #include "TDatabasePDG.h"
28 #include "TLorentzVector.h"
29
30 // aliroot headers
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 #include "AliESDEvent.h"
34 #include "AliAODEvent.h"
35 #include "AliMCEvent.h"
36 #include "AliAODVZERO.h"
37 #include "AliAODZDC.h"
38 #include "AliESDVZERO.h"
39 #include "AliESDZDC.h"
40 #include "AliPIDResponse.h"
41 #include "AliAODTrack.h"
42 #include "AliAODPid.h"
43 #include "AliAODVertex.h"
44 #include "AliESDVertex.h"
45 #include "AliMultiplicity.h"
46 #include "AliESDtrack.h"
47 #include "AliESDMuonTrack.h"
48 #include "AliAODMCParticle.h"
49 #include "AliMCParticle.h"
50
51 // my headers
52 #include "AliAnalysisTaskUpcPsi2s.h"
53
54 ClassImp(AliAnalysisTaskUpcPsi2s);
55
56 using std::cout;
57 using std::endl;
58
59 //trees for UPC analysis,
60 // michal.broz@cern.ch
61
62 //_____________________________________________________________________________
63 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s() 
64   : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
65     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
66     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
67     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
68     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
69     fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
70     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
71     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
72
73 {
74
75 //Dummy constructor
76
77 }//AliAnalysisTaskUpcPsi2s
78
79
80 //_____________________________________________________________________________
81 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name) 
82   : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
83     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
84     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
85     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
86     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
87     fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
88     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
89     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
90
91 {
92
93   // Constructor
94   if( strstr(name,"ESD") ) fType = 0;
95   if( strstr(name,"AOD") ) fType = 1;
96   
97   Init();
98
99   DefineOutput(1, TTree::Class());
100   DefineOutput(2, TTree::Class());
101   DefineOutput(3, TList::Class());
102   DefineOutput(4, TList::Class());
103
104 }//AliAnalysisTaskUpcPsi2s
105
106 //_____________________________________________________________________________
107 void AliAnalysisTaskUpcPsi2s::Init()
108 {
109   
110   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
111
112 }//Init
113
114 //_____________________________________________________________________________
115 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s() 
116 {
117   // Destructor
118   if(fJPsiTree){
119      delete fJPsiTree;
120      fJPsiTree = 0x0;
121   }
122   if(fPsi2sTree){
123      delete fPsi2sTree;
124      fPsi2sTree = 0x0;
125   }
126   if(fListTrig){
127      delete fListTrig;
128      fListTrig = 0x0;
129   }
130   if(fListHist){
131      delete fListHist;
132      fListHist = 0x0;
133   }
134
135 }//~AliAnalysisTaskUpcPsi2s
136
137
138 //_____________________________________________________________________________
139 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
140 {
141   //input file
142   fDataFilnam = new TObjString();
143   fDataFilnam->SetString("");
144
145     //tracks
146   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
147   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
148   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
149   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
150
151   //output tree with JPsi candidate events
152   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
153   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
154   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
155   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
156   
157   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
158   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
159   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
160   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
161   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
162   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
163   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
164   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
165   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
166   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
167   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
168   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
169   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");                        
170   if( fType == 0 ) {
171     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
172   }
173   if( fType == 1 ) {
174     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
175   }
176  
177  //output tree with Psi2s candidate events
178   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
179   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
180   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
181   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
182   
183   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
184   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
185   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
186   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
187   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
188   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
189   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
190   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
191   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
192   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
193   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
194   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
195   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");                       
196   if( fType == 0 ) {
197     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
198   }
199   if( fType == 1 ) {
200     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
201   }
202   
203   fListTrig = new TList();
204   fListTrig ->SetOwner();
205   
206   fHistUpcTriggersPerRun = new TH1D("fHistUpcTriggersPerRun", "fHistUpcTriggersPerRun", 3000, 167000.5, 170000.5);
207   fListTrig->Add(fHistUpcTriggersPerRun);
208   
209   fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 3000, 167000.5, 170000.5);
210   fListTrig->Add(fHistZedTriggersPerRun);
211
212   fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 3000, 167000.5, 170000.5);
213   fListTrig->Add(fHistCvlnTriggersPerRun);
214   
215   fListHist = new TList();
216   fListHist ->SetOwner();
217   
218   TString CutNameJPsi[12] = {"Analyzed","Triggered","Vertex cut","V0 decision","Two good tracks",
219                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
220   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",12,0.5,12.5);
221   for (Int_t i = 0; i<12; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
222   fListHist->Add(fHistNeventsJPsi);
223   
224   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
225   fListHist->Add(fHistTPCsignalJPsi);
226   
227   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
228   fListHist->Add(fHistDiLeptonPtJPsi);
229
230   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
231   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
232   fListHist->Add(fHistDiElectronMass);
233   
234   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
235   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
236   fListHist->Add(fHistDiMuonMass);
237
238   TString CutNamePsi2s[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Four good tracks",
239                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
240
241   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",13,0.5,13.5);
242   for (Int_t i = 0; i<13; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
243   fListHist->Add(fHistNeventsPsi2s);
244
245   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
246   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
247   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
248   fListHist->Add(fHistPsi2sMassVsPt);
249   
250   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
251   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
252   fListHist->Add(fHistPsi2sMassCoherent);
253   
254   PostData(1, fJPsiTree);
255   PostData(2, fPsi2sTree);
256   PostData(3, fListTrig);
257   PostData(4, fListHist);
258
259 }//UserCreateOutputObjects
260
261 //_____________________________________________________________________________
262 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *) 
263 {
264
265   //cout<<"#################### Next event ##################"<<endl;
266
267   if( fType == 0 ){
268         RunESDtrig(); 
269         if(fRunHist) RunESDhist();
270         if(fRunTree) RunESDtree();
271         }
272
273   if( fType == 1 ){
274         RunAODtrig(); 
275         if(fRunHist) RunAODhist();
276         if(fRunTree) RunAODtree();
277         }
278
279 }//UserExec
280 //_____________________________________________________________________________
281 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
282 {
283
284   //input event
285   AliAODEvent *aod = (AliAODEvent*) InputEvent();
286   if(!aod) return;
287
288   fRunNum = aod ->GetRunNumber();
289   //Trigger
290   TString trigger = aod->GetFiredTriggerClasses();
291   
292   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
293   
294   if(trigger.Contains("CVLN-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN-B triggers
295   
296   //if(aod->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
297   
298   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
299   if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
300
301 PostData(3, fListTrig);
302
303 }
304 //_____________________________________________________________________________
305 void AliAnalysisTaskUpcPsi2s::RunAODhist()
306 {
307
308   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
309   
310   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
311   Double_t muonMass = partMuon->Mass();
312   
313   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
314   Double_t electronMass = partElectron->Mass();
315   
316   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
317   Double_t pionMass = partPion->Mass();
318
319   //input event
320   AliAODEvent *aod = (AliAODEvent*) InputEvent();
321   if(!aod) return;
322
323   fHistNeventsJPsi->Fill(1);
324   fHistNeventsPsi2s->Fill(1);
325
326   //Trigger
327   TString trigger = aod->GetFiredTriggerClasses();
328   
329   if( !trigger.Contains("CCUP4-B") ) return;
330   
331   fHistNeventsJPsi->Fill(2);
332   fHistNeventsPsi2s->Fill(2);
333
334   //primary vertex
335   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
336   fVtxContrib = fAODVertex->GetNContributors();
337   if(fVtxContrib < 2) return;
338   
339   fHistNeventsJPsi->Fill(3);
340   fHistNeventsPsi2s->Fill(3);
341
342   //VZERO, ZDC
343   AliAODVZERO *fV0data = aod ->GetVZEROData();
344   AliAODZDC *fZDCdata = aod->GetZDCData();
345   
346   fV0Adecision = fV0data->GetV0ADecision();
347   fV0Cdecision = fV0data->GetV0CDecision();
348   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
349   
350   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
351   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
352
353   if( fZDCAenergy > 12000 || fZDCCenergy > 12000) return;
354   
355   fHistNeventsJPsi->Fill(4);
356   fHistNeventsPsi2s->Fill(4);
357
358    Int_t nGoodTracks=0;
359   //Two tracks loop
360   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
361   
362   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
363   Short_t qLepton[4], qPion[4];
364   UInt_t nLepton=0, nPion=0, nHighPt=0;
365   Double_t jRecTPCsignal[5];
366   Int_t mass[3]={-1,-1,-1};
367   
368   //Two track loop
369   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
370     AliAODTrack *trk = aod->GetTrack(itr);
371     if( !trk ) continue;
372
373       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
374       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
375       if(trk->GetTPCNcls() < 70)continue;
376       if(trk->Chi2perNDF() > 4)continue;
377       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
378       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
379       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
380       if(TMath::Abs(dca[1]) > 2) continue;
381       if(TMath::Abs(dca[0]) > 0.2) continue;
382      
383       TrackIndex[nGoodTracks] = itr;
384       nGoodTracks++;
385                                   
386       if(nGoodTracks > 2) break;  
387   }//Track loop
388
389   if(nGoodTracks == 2){
390           fHistNeventsJPsi->Fill(5);
391           for(Int_t i=0; i<2; i++){
392                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);                
393                 if(trk->Pt() > 1) nHighPt++;     
394                 jRecTPCsignal[nLepton] = trk->GetTPCsignal();     
395                 qLepton[nLepton] = trk->Charge();
396                 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
397                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
398                                 mass[nLepton] = 0;
399                                 }
400                 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
401                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
402                                 mass[nLepton] = 1;
403                                 }
404                 nLepton++;              
405                 }               
406         if(nLepton == 2){
407                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
408                 if(qLepton[0]*qLepton[1] < 0){
409                         fHistNeventsJPsi->Fill(7);
410                         if(nHighPt > 0){
411                                 fHistNeventsJPsi->Fill(8);
412                                 fHistTPCsignalJPsi->Fill(jRecTPCsignal[0],jRecTPCsignal[1]);
413                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
414                                 if(mass[0] == mass[1] && mass[0] != -1) {
415                                         fHistNeventsJPsi->Fill(10);
416                                         vCandidate = vLepton[0]+vLepton[1];
417                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
418                                         if(mass[0] == 0) {
419                                                 fHistDiMuonMass->Fill(vCandidate.M());
420                                                 fHistNeventsJPsi->Fill(11);
421                                                 }
422                                         if(mass[0] == 1) {
423                                                 fHistDiElectronMass->Fill(vCandidate.M());
424                                                 fHistNeventsJPsi->Fill(12);
425                                                 }
426                                         }
427                                 }
428                         }
429                 }
430   }
431   
432   nGoodTracks = 0;
433   //Four track loop
434   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
435     AliAODTrack *trk = aod->GetTrack(itr);
436     if( !trk ) continue;
437
438       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
439       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
440       if(trk->GetTPCNcls() < 50)continue;
441       if(trk->Chi2perNDF() > 4)continue;
442       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
443       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
444       if(TMath::Abs(dca[1]) > 2) continue;
445      
446       TrackIndex[nGoodTracks] = itr;
447       nGoodTracks++;
448                                   
449       if(nGoodTracks > 4) break;  
450   }//Track loop
451   
452   nLepton=0; nPion=0; nHighPt=0;
453   mass[0]= -1; mass[1]= -1, mass[2]= -1;
454   
455   if(nGoodTracks == 4){
456           fHistNeventsPsi2s->Fill(5);
457           for(Int_t i=0; i<4; i++){
458                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
459                 
460                 if(trk->Pt() > 1){   
461                         jRecTPCsignal[nLepton] = trk->GetTPCsignal();      
462                         qLepton[nLepton] = trk->Charge();
463                         if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
464                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
465                                         mass[nLepton] = 0;
466                                         }
467                         if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
468                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
469                                         mass[nLepton] = 1;
470                                         }
471                         nLepton++;
472                         }
473                 else{
474                         qPion[nPion] = trk->Charge();
475                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
476                         nPion++;
477                         }
478                                       
479                 if(nLepton > 2 || nPion > 2) break;
480                 }
481         if((nLepton == 2) && (nPion == 2)){
482                 fHistNeventsPsi2s->Fill(6);
483                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
484                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
485                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
486                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
487                         fHistNeventsPsi2s->Fill(10);
488                         if(mass[0] == mass[1]) {
489                                 fHistNeventsPsi2s->Fill(11); 
490                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
491                                 vDilepton = vLepton[0]+vLepton[1];
492                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
493                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
494                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
495                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
496                                 }
497                         }
498                 }
499   }
500   
501   PostData(4, fListHist);
502
503 }
504
505 //_____________________________________________________________________________
506 void AliAnalysisTaskUpcPsi2s::RunAODtree()
507 {
508   //input event
509   AliAODEvent *aod = (AliAODEvent*) InputEvent();
510   if(!aod) return;
511
512   //input data
513   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
514   fDataFilnam->Clear();
515   fDataFilnam->SetString(filnam);
516   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
517   fRunNum = aod ->GetRunNumber();
518
519   //Trigger
520   TString trigger = aod->GetFiredTriggerClasses();
521   
522   if( !trigger.Contains("CCUP4-B") ) return;
523
524   //trigger inputs
525   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
526   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
527
528   //Event identification
529   fPerNum = aod ->GetPeriodNumber();
530   fOrbNum = aod ->GetOrbitNumber();
531   fBCrossNum = aod ->GetBunchCrossNumber();
532
533   //primary vertex
534   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
535   fVtxContrib = fAODVertex->GetNContributors();
536
537   //Tracklets
538   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
539
540   //VZERO, ZDC
541   AliAODVZERO *fV0data = aod ->GetVZEROData();
542   AliAODZDC *fZDCdata = aod->GetZDCData();
543   
544   fV0Adecision = fV0data->GetV0ADecision();
545   fV0Cdecision = fV0data->GetV0CDecision();
546   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
547   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
548   
549   Int_t nGoodTracks=0;
550   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
551
552   //Two track loop
553   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
554     AliAODTrack *trk = aod->GetTrack(itr);
555     if( !trk ) continue;
556
557       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
558       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
559       if(trk->GetTPCNcls() < 70)continue;
560       if(trk->Chi2perNDF() > 4)continue;
561       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
562       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
563       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
564       if(TMath::Abs(dca[1]) > 2) continue;
565       if(TMath::Abs(dca[0]) > 0.2) continue;
566      
567       TrackIndex[nGoodTracks] = itr;
568       nGoodTracks++;
569                                   
570       if(nGoodTracks > 2) break;  
571   }//Track loop
572
573   if(nGoodTracks == 2){
574           for(Int_t i=0; i<2; i++){
575                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
576                 
577                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
578                 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
579                                 
580                 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
581                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk); 
582                 
583                 }
584   fJPsiTree ->Fill();
585   PostData(1, fJPsiTree);
586   }
587
588   nGoodTracks = 0;
589   //Four track loop
590   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
591     AliAODTrack *trk = aod->GetTrack(itr);
592     if( !trk ) continue;
593
594       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
595       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
596       if(trk->GetTPCNcls() < 50)continue;
597       if(trk->Chi2perNDF() > 4)continue;
598       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
599       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
600       if(TMath::Abs(dca[1]) > 2) continue;
601      
602       TrackIndex[nGoodTracks] = itr;
603       nGoodTracks++;
604                                   
605       if(nGoodTracks > 4) break;  
606   }//Track loop
607   
608   
609   if(nGoodTracks == 4){
610           for(Int_t i=0; i<4; i++){
611                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
612                 
613                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
614                 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
615                 
616                 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
617                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk); 
618                 
619                 }
620   fPsi2sTree ->Fill();
621   PostData(2, fPsi2sTree);
622   }
623
624 }//RunAOD
625
626 //_____________________________________________________________________________
627 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
628 {
629
630   //input event
631   AliESDEvent *esd = (AliESDEvent*) InputEvent();
632   if(!esd) return;
633
634   fRunNum = esd ->GetRunNumber();
635   //Trigger
636   TString trigger = esd->GetFiredTriggerClasses();
637   
638   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
639   
640   if(trigger.Contains("CVLN-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN-B triggers
641   
642   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
643   
644 PostData(3, fListTrig);
645
646 }
647 //_____________________________________________________________________________
648 void AliAnalysisTaskUpcPsi2s::RunESDhist()
649 {
650
651   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
652   
653   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
654   Double_t muonMass = partMuon->Mass();
655   
656   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
657   Double_t electronMass = partElectron->Mass();
658   
659   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
660   Double_t pionMass = partPion->Mass();
661
662   //input event
663   AliESDEvent *esd = (AliESDEvent*) InputEvent();
664   if(!esd) return;
665
666   fHistNeventsJPsi->Fill(1);
667   fHistNeventsPsi2s->Fill(1);
668
669   //Trigger
670   TString trigger = esd->GetFiredTriggerClasses();
671   
672   if( !trigger.Contains("CCUP4-B") ) return;
673   
674   fHistNeventsJPsi->Fill(2);
675   fHistNeventsPsi2s->Fill(2);
676
677   //primary vertex
678   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
679   fVtxContrib = fESDVertex->GetNContributors();
680   if(fVtxContrib < 2) return;
681   
682   fHistNeventsJPsi->Fill(3);
683   fHistNeventsPsi2s->Fill(3);
684
685   //VZERO, ZDC
686   AliESDVZERO *fV0data = esd->GetVZEROData();
687   AliESDZDC *fZDCdata = esd->GetESDZDC();
688   
689   fV0Adecision = fV0data->GetV0ADecision();
690   fV0Cdecision = fV0data->GetV0CDecision();
691   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
692   
693   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
694   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
695   if( fZDCAenergy > 12000 || fZDCCenergy > 12000) return;
696   
697   fHistNeventsJPsi->Fill(4);
698   fHistNeventsPsi2s->Fill(4);
699
700    Int_t nGoodTracks=0;
701   //Two tracks loop
702   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
703   
704   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
705   Short_t qLepton[4], qPion[4];
706   UInt_t nLepton=0, nPion=0, nHighPt=0;
707   Double_t jRecTPCsignal[5];
708   Int_t mass[3]={-1,-1,-1};
709
710   //Track loop
711   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
712     AliESDtrack *trk = esd->GetTrack(itr);
713     if( !trk ) continue;
714
715       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
716       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
717       if(trk->GetTPCNcls() < 50)continue;
718       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
719       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
720       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
721       trk->GetImpactParameters(dca[0],dca[1]);
722       if(TMath::Abs(dca[1]) > 2) continue;
723       
724       TrackIndex[nGoodTracks] = itr;
725       nGoodTracks++;
726       if(nGoodTracks > 4) break;   
727   }//Track loop
728   
729   if(nGoodTracks == 2){
730           fHistNeventsJPsi->Fill(5);
731           for(Int_t i=0; i<2; i++){
732                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);                
733                 if(trk->Pt() > 1) nHighPt++;     
734                 jRecTPCsignal[nLepton] = trk->GetTPCsignal();     
735                 qLepton[nLepton] = trk->Charge();
736                 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
737                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
738                                 mass[nLepton] = 0;
739                                 }
740                 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
741                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
742                                 mass[nLepton] = 1;
743                                 }
744                 nLepton++;              
745                 }               
746         if(nLepton == 2){
747                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
748                 if(qLepton[0]*qLepton[1] < 0){
749                         fHistNeventsJPsi->Fill(7);
750                         if(nHighPt > 0){
751                                 fHistNeventsJPsi->Fill(8);
752                                 fHistTPCsignalJPsi->Fill(jRecTPCsignal[0],jRecTPCsignal[1]);
753                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
754                                 if(mass[0] == mass[1] && mass[0] != -1) {
755                                         fHistNeventsJPsi->Fill(10);
756                                         vCandidate = vLepton[0]+vLepton[1];
757                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
758                                         if(mass[0] == 0) {
759                                                 fHistDiMuonMass->Fill(vCandidate.M());
760                                                 fHistNeventsJPsi->Fill(11);
761                                                 }
762                                         if(mass[0] == 1) {
763                                                 fHistDiElectronMass->Fill(vCandidate.M());
764                                                 fHistNeventsJPsi->Fill(12);
765                                                 }
766                                         }
767                                 }
768                         }
769                 }
770   }
771   nLepton=0; nPion=0; nHighPt=0;
772   mass[0]= -1; mass[1]= -1, mass[2]= -1;
773   
774   if(nGoodTracks == 4){
775           fHistNeventsPsi2s->Fill(5);
776           for(Int_t i=0; i<4; i++){
777                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
778                 
779                 if(trk->Pt() > 1){   
780                         jRecTPCsignal[nLepton] = trk->GetTPCsignal();      
781                         qLepton[nLepton] = trk->Charge();
782                         if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
783                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
784                                         mass[nLepton] = 0;
785                                         }
786                         if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
787                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
788                                         mass[nLepton] = 1;
789                                         }
790                         nLepton++;
791                         }
792                 else{
793                         qPion[nPion] = trk->Charge();
794                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
795                         nPion++;
796                         }
797                                       
798                 if(nLepton > 2 || nPion > 2) break;
799                 }
800         if((nLepton == 2) && (nPion == 2)){
801                 fHistNeventsPsi2s->Fill(6);
802                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
803                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
804                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
805                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
806                         fHistNeventsPsi2s->Fill(10);
807                         if(mass[0] == mass[1]) {
808                                 fHistNeventsPsi2s->Fill(11); 
809                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
810                                 vDilepton = vLepton[0]+vLepton[1];
811                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
812                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
813                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
814                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
815                                 }
816                         }
817                 }
818   }
819   
820   PostData(4, fListHist);
821
822 }
823
824 //_____________________________________________________________________________
825 void AliAnalysisTaskUpcPsi2s::RunESDtree()
826 {
827
828   //input event
829   AliESDEvent *esd = (AliESDEvent*) InputEvent();
830   if(!esd) return;
831
832   //input data
833   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
834   fDataFilnam->Clear();
835   fDataFilnam->SetString(filnam);
836   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
837   fRunNum = esd->GetRunNumber();
838
839   //Trigger
840   TString trigger = esd->GetFiredTriggerClasses();
841   
842   if( !trigger.Contains("CCUP4-B") ) return;
843
844   //trigger inputs
845   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
846   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
847
848   //Event identification
849   fPerNum = esd->GetPeriodNumber();
850   fOrbNum = esd->GetOrbitNumber();
851   fBCrossNum = esd->GetBunchCrossNumber();
852
853   //primary vertex
854   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
855   fVtxContrib = fESDVertex->GetNContributors();
856
857   //Tracklets
858   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
859
860   //VZERO, ZDC
861   AliESDVZERO *fV0data = esd->GetVZEROData();
862   AliESDZDC *fZDCdata = esd->GetESDZDC();
863   
864   fV0Adecision = fV0data->GetV0ADecision();
865   fV0Cdecision = fV0data->GetV0CDecision();
866   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
867   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
868   
869   Int_t nGoodTracks=0;
870   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
871   
872   //Track loop
873   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
874     AliESDtrack *trk = esd->GetTrack(itr);
875     if( !trk ) continue;
876
877       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
878       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
879       if(trk->GetTPCNcls() < 50)continue;
880       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
881       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
882       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
883       trk->GetImpactParameters(dca[0],dca[1]);
884       if(TMath::Abs(dca[1]) > 2) continue;
885       
886       TrackIndex[nGoodTracks] = itr;
887       nGoodTracks++;
888       if(nGoodTracks > 4) break;   
889   }//Track loop
890
891   if(nGoodTracks == 2){
892           for(Int_t i=0; i<2; i++){
893                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
894                 
895                 AliExternalTrackParam cParam;
896                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
897                                 
898                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk); 
899                 
900                 }
901   fJPsiTree ->Fill();
902   PostData(1, fJPsiTree);
903   }
904   
905   if(nGoodTracks == 4){
906           for(Int_t i=0; i<4; i++){
907                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
908                 
909                 AliExternalTrackParam cParam;
910                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
911
912                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk); 
913                 
914                 }
915   fPsi2sTree ->Fill();
916   PostData(2, fPsi2sTree);
917   }
918
919 }//RunESD
920
921 //_____________________________________________________________________________
922 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *) 
923 {
924
925   cout<<"Analysis complete."<<endl;
926 }//Terminate
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956