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