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