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),hCounter(0),fJPsiTree(0),fPsi2sTree(0),
  65     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
  66     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  67     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  68     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
  69     fListQA(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
  70     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
  71     fHistK0sMass(0)
  72
  73 {
  74
  75 //Dummy constructor
  76
  77 }//AliAnalysisTaskUpcPsi2s
  78
  79
  80 //_____________________________________________________________________________
  81 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
  82   : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),hCounter(0),fJPsiTree(0),fPsi2sTree(0),
  83     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
  84     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  85     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  86     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
  87     fListQA(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
  88     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
  89     fHistK0sMass(0)
  90
  91 {
  92
  93   // Constructor
  94   if( strstr(name,"ESD") ) fType = 0;
  95   if( strstr(name,"AOD") ) fType = 1;
  96
  97   Init();
  98
  99   DefineOutput(1, TTree::Class());
 100   DefineOutput(2, TTree::Class());
 101   DefineOutput(3, TH1I::Class());
 102   DefineOutput(4, TList::Class());
 103
 104 }//AliAnalysisTaskUpcPsi2s
 105
 106 //_____________________________________________________________________________
 107 void AliAnalysisTaskUpcPsi2s::Init()
 108 {
 109
 110   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
 111
 112 }//Init
 113
 114 //_____________________________________________________________________________
 115 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
 116 {
 117   // Destructor
 118   if(fJPsiTree){
 119      delete fJPsiTree;
 120      fJPsiTree = 0x0;
 121   }
 122   if(fPsi2sTree){
 123      delete fPsi2sTree;
 124      fPsi2sTree = 0x0;
 125   }
 126   if(hCounter){
 127      delete hCounter;
 128      hCounter = 0x0;
 129   }
 130   if(fListQA){
 131      delete fListQA;
 132      fListQA = 0x0;
 133   }
 134
 135 }//~AliAnalysisTaskUpcPsi2s
 136
 137
 138 //_____________________________________________________________________________
 139 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
 140 {
 141    hCounter = new TH1I("hCounter", "hCounter", 34000, 1., 34001.);
 142
 143   //input file
 144   fDataFilnam = new TObjString();
 145   fDataFilnam->SetString("");
 146
 147     //tracks
 148   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
 149   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
 150   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
 151   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
 152
 153   //output tree with JPsi candidate events
 154   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
 155   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 156   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 157   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 158
 159   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 160   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 161   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 162   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 163   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 164   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 165   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 166   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 167   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 168   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 169   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
 170   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 171   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 172   if( fType == 0 ) {
 173     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
 174   }
 175   if( fType == 1 ) {
 176     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
 177   }
 178
 179  //output tree with Psi2s candidate events
 180   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
 181   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 182   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 183   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 184
 185   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 186   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 187   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 188   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 189   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 190   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 191   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 192   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 193   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 194   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 195   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
 196   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 197   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 198   if( fType == 0 ) {
 199     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
 200   }
 201   if( fType == 1 ) {
 202     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
 203   }
 204
 205   fListQA = new TList();
 206   fListQA ->SetOwner();
 207
 208   TString CutNameJPsi[12] = {"Analyzed","Triggered","Vertex cut","V0 decision","Two good tracks",
 209                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
 210   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",12,0.5,12.5);
 211   for (Int_t i = 0; i<12; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
 212   fListQA->Add(fHistNeventsJPsi);
 213
 214   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
 215   fListQA->Add(fHistTPCsignalJPsi);
 216
 217   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
 218   fListQA->Add(fHistDiLeptonPtJPsi);
 219
 220   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
 221   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
 222   fListQA->Add(fHistDiElectronMass);
 223
 224   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
 225   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
 226   fListQA->Add(fHistDiMuonMass);
 227
 228   TString CutNamePsi2s[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Four good tracks",
 229                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
 230
 231   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",13,0.5,13.5);
 232   for (Int_t i = 0; i<13; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
 233   fListQA->Add(fHistNeventsPsi2s);
 234
 235   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
 236   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 237   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 238   fListQA->Add(fHistPsi2sMassVsPt);
 239
 240   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
 241   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 242   fListQA->Add(fHistPsi2sMassCoherent);
 243
 244   fHistK0sMass = new TH1D("fHistK0sMass","fHistK0sMass",200,0.4,0.6);
 245   fListQA->Add(fHistK0sMass);
 246
 247   PostData(1, fJPsiTree);
 248   PostData(2, fPsi2sTree);
 249   PostData(3, hCounter);
 250   PostData(4, fListQA);
 251
 252 }//UserCreateOutputObjects
 253
 254 //_____________________________________________________________________________
 255 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
 256 {
 257
 258   //cout<<"#################### Next event ##################"<<endl;
 259
 260   if( fType == 0 ) RunESD();
 261   if( fType == 1 ){
 262         if(fRunHist) RunAODhist();
 263         if(fRunTree) RunAODtree();
 264         }
 265
 266 }//UserExec
 267 //_____________________________________________________________________________
 268 void AliAnalysisTaskUpcPsi2s::RunAODhist()
 269 {
 270
 271   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 272
 273   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 274   Double_t muonMass = partMuon->Mass();
 275
 276   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 277   Double_t electronMass = partElectron->Mass();
 278
 279   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 280   Double_t pionMass = partPion->Mass();
 281
 282   //input event
 283   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 284   if(!aod) return;
 285
 286   fHistNeventsJPsi->Fill(1);
 287   fHistNeventsPsi2s->Fill(1);
 288
 289   //Trigger
 290   TString trigger = aod->GetFiredTriggerClasses();
 291
 292   if( !trigger.Contains("CCUP4-B") ) return;
 293
 294   fHistNeventsJPsi->Fill(2);
 295   fHistNeventsPsi2s->Fill(2);
 296
 297   //primary vertex
 298   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 299   fVtxContrib = fAODVertex->GetNContributors();
 300   if(fVtxContrib < 2) return;
 301
 302   fHistNeventsJPsi->Fill(3);
 303   fHistNeventsPsi2s->Fill(3);
 304
 305
 306   //VZERO, ZDC
 307   AliAODVZERO *fV0data = aod ->GetVZEROData();
 308   //AliAODZDC *fZDCdata = aod->GetZDCData();
 309
 310   fV0Adecision = fV0data->GetV0ADecision();
 311   fV0Cdecision = fV0data->GetV0CDecision();
 312   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 313
 314   fHistNeventsJPsi->Fill(4);
 315   fHistNeventsPsi2s->Fill(4);
 316
 317    Int_t nGoodTracks=0;
 318   //Two tracks loop
 319   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 320
 321   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 322   Short_t qLepton[4], qPion[4];
 323   UInt_t nLepton=0, nPion=0, nHighPt=0;
 324   Double_t jRecTPCsignal[5];
 325   Int_t mass[3]={-1,-1,-1};
 326
 327   //Track loop
 328   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 329     AliAODTrack *trk = aod->GetTrack(itr);
 330     if( !trk ) continue;
 331
 332       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 333       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 334       if(trk->GetTPCNcls() < 50)continue;
 335       if(trk->Chi2perNDF() > 4)continue;
 336       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 337       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 338       if(TMath::Abs(dca[1]) > 2) continue;
 339
 340       TrackIndex[nGoodTracks] = itr;
 341       nGoodTracks++;
 342
 343       if(nGoodTracks > 4) break;
 344   }//Track loop
 345
 346   if(nGoodTracks == 2){
 347           fHistNeventsJPsi->Fill(5);
 348           for(Int_t i=0; i<2; i++){
 349                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 350                 if(trk->Pt() > 1) nHighPt++;
 351                 jRecTPCsignal[nLepton] = trk->GetTPCsignal();
 352                 qLepton[nLepton] = trk->Charge();
 353                 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
 354                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 355                                 mass[nLepton] = 0;
 356                                 }
 357                 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
 358                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 359                                 mass[nLepton] = 1;
 360                                 }
 361                 nLepton++;
 362                 }
 363         if(nLepton == 2){
 364                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
 365                 if(qLepton[0]*qLepton[1] < 0){
 366                         fHistNeventsJPsi->Fill(7);
 367                         if(nHighPt > 0){
 368                                 fHistNeventsJPsi->Fill(8);
 369                                 fHistTPCsignalJPsi->Fill(jRecTPCsignal[0],jRecTPCsignal[1]);
 370                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
 371                                 if(mass[0] == mass[1] && mass[0] != -1) {
 372                                         fHistNeventsJPsi->Fill(10);
 373                                         vCandidate = vLepton[0]+vLepton[1];
 374                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
 375                                         if(mass[0] == 0) {
 376                                                 fHistDiMuonMass->Fill(vCandidate.M());
 377                                                 fHistNeventsJPsi->Fill(11);
 378                                                 }
 379                                         if(mass[0] == 1) {
 380                                                 fHistDiElectronMass->Fill(vCandidate.M());
 381                                                 fHistNeventsJPsi->Fill(12);
 382                                                 }
 383                                         }
 384                                 }
 385                         }
 386                 }
 387   }
 388   nLepton=0; nPion=0; nHighPt=0;
 389   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 390
 391   if(nGoodTracks == 4){
 392           fHistNeventsPsi2s->Fill(5);
 393           for(Int_t i=0; i<4; i++){
 394                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 395
 396                 if(trk->Pt() > 1){
 397                         jRecTPCsignal[nLepton] = trk->GetTPCsignal();
 398                         qLepton[nLepton] = trk->Charge();
 399                         if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
 400                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 401                                         mass[nLepton] = 0;
 402                                         }
 403                         if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
 404                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 405                                         mass[nLepton] = 1;
 406                                         }
 407                         nLepton++;
 408                         }
 409                 else{
 410                         qPion[nPion] = trk->Charge();
 411                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
 412                         nPion++;
 413                         }
 414
 415                 if(nLepton > 2 || nPion > 2) break;
 416                 }
 417         if((nLepton == 2) && (nPion == 2)){
 418                 fHistNeventsPsi2s->Fill(6);
 419                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
 420                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
 421                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
 422                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
 423                         fHistNeventsPsi2s->Fill(10);
 424                         if(mass[0] == mass[1]) {
 425                                 fHistNeventsPsi2s->Fill(11);
 426                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
 427                                 vDilepton = vLepton[0]+vLepton[1];
 428                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
 429                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 430                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
 431                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
 432                                 }
 433                         }
 434                 }
 435   }
 436
 437
 438    //---------------------------------K0s + K0s loop - very experimental--------------------
 439    nGoodTracks = 0;
 440   //V0s loop
 441   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 442     AliAODv0 *v0 = aod->GetV0(iV0);
 443     if( !v0 ) continue;
 444     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 445     if (lOnFlyStatus) continue;
 446
 447     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 448     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 449     if (!pTrack || !nTrack) continue;
 450
 451     if ( pTrack->Charge() == nTrack->Charge())continue;
 452
 453       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 454       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 455       if(pTrack->GetTPCNcls() < 50)continue;
 456       if(nTrack->GetTPCNcls() < 50)continue;
 457       if(pTrack->Chi2perNDF() > 4)continue;
 458       if(nTrack->Chi2perNDF() > 4)continue;
 459
 460       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 461       if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 462       if(TMath::Abs(dca[1]) > 2) continue;
 463       if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 464       if(TMath::Abs(dca[1]) > 2) continue;
 465
 466       TrackIndex[nGoodTracks] = iV0;
 467       nGoodTracks++;
 468       if(nGoodTracks > 2) break;
 469   }//V0s loop
 470   if(nGoodTracks == 2){
 471         for(Int_t i=0; i<2; i++){
 472                 AliAODv0 *v0 = aod->GetV0(TrackIndex[i]);
 473                 fHistK0sMass->Fill(v0->MassK0Short());
 474                 }
 475   }
 476
 477   PostData(4, fListQA);
 478
 479 }
 480
 481 //_____________________________________________________________________________
 482 void AliAnalysisTaskUpcPsi2s::RunAODtree()
 483 {
 484   //input event
 485   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 486   if(!aod) return;
 487
 488   //input data
 489   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 490   fDataFilnam->Clear();
 491   fDataFilnam->SetString(filnam);
 492   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 493   fRunNum = aod ->GetRunNumber();
 494
 495   hCounter->Fill( 1 );
 496
 497   //Trigger
 498   TString trigger = aod->GetFiredTriggerClasses();
 499
 500   fTrigger[0]   = trigger.Contains("CINT7-B");
 501   fTrigger[1]   = trigger.Contains("CCUP4-B"); // CE
 502   fTrigger[2]   = trigger.Contains("CCUP4-E"); // CE
 503
 504   Bool_t isTRG = kFALSE;
 505   for(Int_t i=1; i<ntrg; i++) {
 506     if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
 507   }
 508   if( !isTRG ) {PostData(3, hCounter); return;}
 509
 510   hCounter->Fill( 2 );
 511
 512   //trigger inputs
 513   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 514   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 515
 516   //Event identification
 517   fPerNum = aod ->GetPeriodNumber();
 518   fOrbNum = aod ->GetOrbitNumber();
 519   fBCrossNum = aod ->GetBunchCrossNumber();
 520
 521   //primary vertex
 522   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 523   fVtxContrib = fAODVertex->GetNContributors();
 524
 525   //Tracklets
 526   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 527
 528   //VZERO, ZDC
 529   AliAODVZERO *fV0data = aod ->GetVZEROData();
 530   AliAODZDC *fZDCdata = aod->GetZDCData();
 531
 532   fV0Adecision = fV0data->GetV0ADecision();
 533   fV0Cdecision = fV0data->GetV0CDecision();
 534   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 535   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 536
 537   Int_t nGoodTracks=0;
 538   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 539
 540   //Track loop
 541   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 542     AliAODTrack *trk = aod->GetTrack(itr);
 543     if( !trk ) continue;
 544
 545       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 546       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 547       if(trk->GetTPCNcls() < 50)continue;
 548       if(trk->Chi2perNDF() > 4)continue;
 549       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 550       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 551       if(TMath::Abs(dca[1]) > 2) continue;
 552
 553       TrackIndex[nGoodTracks] = itr;
 554       nGoodTracks++;
 555
 556       if(nGoodTracks > 4) break;
 557   }//Track loop
 558
 559   if(nGoodTracks == 2){
 560           for(Int_t i=0; i<2; i++){
 561                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 562
 563                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 564                 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
 565
 566                 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
 567                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
 568
 569                 }
 570   fJPsiTree ->Fill();
 571   PostData(1, fJPsiTree);
 572   }
 573
 574   if(nGoodTracks == 4){
 575           for(Int_t i=0; i<4; i++){
 576                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 577
 578                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 579                 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
 580
 581                 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
 582                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
 583
 584                 }
 585   fPsi2sTree ->Fill();
 586   PostData(2, fPsi2sTree);
 587   }
 588
 589   PostData(3, hCounter);
 590
 591 }//RunAOD
 592
 593 //_____________________________________________________________________________
 594 void AliAnalysisTaskUpcPsi2s::RunESD()
 595 {
 596
 597   //input event
 598   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 599   if(!esd) return;
 600
 601   //input data
 602   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 603   fDataFilnam->Clear();
 604   fDataFilnam->SetString(filnam);
 605   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 606   fRunNum = esd->GetRunNumber();
 607
 608   hCounter->Fill( 1 );
 609
 610   //Trigger
 611   TString trigger = esd->GetFiredTriggerClasses();
 612
 613   fTrigger[0]   = trigger.Contains("CINT7-B");
 614   fTrigger[1]   = trigger.Contains("CCUP4-B"); // CE
 615   fTrigger[2]   = trigger.Contains("CCUP4-E"); // CE
 616
 617   Bool_t isTRG = kFALSE;
 618   for(Int_t i=1; i<ntrg; i++) {
 619     if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
 620   }
 621   if( !isTRG ) {PostData(3, hCounter); return;}
 622
 623   hCounter->Fill( 2 );
 624
 625   //trigger inputs
 626   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
 627   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
 628
 629   //Event identification
 630   fPerNum = esd->GetPeriodNumber();
 631   fOrbNum = esd->GetOrbitNumber();
 632   fBCrossNum = esd->GetBunchCrossNumber();
 633
 634   //primary vertex
 635   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
 636   fVtxContrib = fESDVertex->GetNContributors();
 637
 638   //Tracklets
 639   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
 640
 641   //VZERO, ZDC
 642   AliESDVZERO *fV0data = esd->GetVZEROData();
 643   AliESDZDC *fZDCdata = esd->GetESDZDC();
 644
 645   fV0Adecision = fV0data->GetV0ADecision();
 646   fV0Cdecision = fV0data->GetV0CDecision();
 647   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
 648   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
 649
 650   Int_t nGoodTracks=0;
 651   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 652
 653   //Track loop
 654   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
 655     AliESDtrack *trk = esd->GetTrack(itr);
 656     if( !trk ) continue;
 657
 658       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 659       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 660       if(trk->GetTPCNcls() < 50)continue;
 661       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
 662       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
 663       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
 664       trk->GetImpactParameters(dca[0],dca[1]);
 665       if(TMath::Abs(dca[1]) > 2) continue;
 666
 667       TrackIndex[nGoodTracks] = itr;
 668       nGoodTracks++;
 669       if(nGoodTracks > 4) break;
 670   }//Track loop
 671
 672   if(nGoodTracks == 2){
 673           for(Int_t i=0; i<2; i++){
 674                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
 675
 676                 AliExternalTrackParam cParam;
 677                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
 678
 679                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
 680
 681                 }
 682   fJPsiTree ->Fill();
 683   PostData(1, fJPsiTree);
 684   }
 685
 686   if(nGoodTracks == 4){
 687           for(Int_t i=0; i<4; i++){
 688                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
 689
 690                 AliExternalTrackParam cParam;
 691                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
 692
 693                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
 694
 695                 }
 696   fPsi2sTree ->Fill();
 697   PostData(2, fPsi2sTree);
 698   }
 699
 700   PostData(3, hCounter);
 701
 702 }//RunESD
 703
 704 //_____________________________________________________________________________
 705 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
 706 {
 707
 708   cout<<"Analysis complete."<<endl;
 709 }//Terminate
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739