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