[u/mrichter/AliRoot.git] / PWGUD / UPC / AliAnalysisTaskUpcK0sK0s.cxx

/*************************************************************************
* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  * 
**************************************************************************/

// c++ headers
#include <iostream>
#include <string.h>
#include <bitset>

// root headers
#include "TH1I.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TParticle.h"
#include "TString.h"
#include "TFile.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"

// aliroot headers
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliMCEvent.h"
#include "AliAODVZERO.h"
#include "AliAODZDC.h"
#include "AliESDVZERO.h"
#include "AliESDZDC.h"
#include "AliPIDResponse.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "AliAODVertex.h"
#include "AliESDVertex.h"
#include "AliMultiplicity.h"
#include "AliESDtrack.h"
#include "AliESDMuonTrack.h"
#include "AliAODMCParticle.h"
#include "AliMCParticle.h"

// my headers
#include "AliAnalysisTaskUpcK0sK0s.h"

ClassImp(AliAnalysisTaskUpcK0sK0s);

using std::cout;
using std::endl;

//trees for UPC analysis,
// michal.broz@cern.ch

//_____________________________________________________________________________
AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s() 
  : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
    fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
    fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
    fDataFilnam(0),fRecoPass(0),fEvtNum(0),
    fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
    fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)

{

//Dummy constructor

}//AliAnalysisTaskUpcK0sK0s


//_____________________________________________________________________________
AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s(const char *name) 
  : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
    fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
    fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
    fDataFilnam(0),fRecoPass(0),fEvtNum(0),
    fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
    fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)

{

  // Constructor
  if( strstr(name,"ESD") ) fType = 0;
  if( strstr(name,"AOD") ) fType = 1;
  
  Init();

  DefineOutput(1, TTree::Class());
  DefineOutput(2, TList::Class());

}//AliAnalysisTaskUpcK0sK0s

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::Init()
{
  
  for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;

}//Init

//_____________________________________________________________________________
AliAnalysisTaskUpcK0sK0s::~AliAnalysisTaskUpcK0sK0s() 
{
  // Destructor
  if(fK0sTree){
     delete fK0sTree;
     fK0sTree = 0x0;
  }
  if(fListHist){
     delete fListHist;
     fListHist = 0x0;
  }

}//~AliAnalysisTaskUpcK0sK0s


//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::UserCreateOutputObjects()
{

    //vertices
  fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
  fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
  
    //tracks
  fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
  fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);

  //output tree with K0s candidate events
  fK0sTree = new TTree("fK0sTree", "fK0sTree");
  fK0sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
  fK0sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
  fK0sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
  
  fK0sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
  fK0sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
  fK0sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
  fK0sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
  fK0sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
  fK0sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
  fK0sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
  fK0sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
  fK0sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
  fK0sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
  fK0sTree ->Branch("fDataFilnam", &fDataFilnam);
  fK0sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
  fK0sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L"); 		       
  if( fType == 0 ) {
    fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
    fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
  }
  if( fType == 1 ) {
    fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
    fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
  }
 
  
  fListHist = new TList();
  fListHist ->SetOwner();
  
  fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
  fListHist->Add(fHistTriggersPerRun);
    
  fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
  fListHist->Add(fHistK0sMass);
  
  PostData(1, fK0sTree);
  PostData(2, fListHist);

}//UserCreateOutputObjects

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *) 
{

  //cout<<"#################### Next event ##################"<<endl;

  if( fType == 0 ) RunESD();
  if( fType == 1 ){ 
  	if(fRunHist) RunAODhist();
	if(fRunTree) RunAODtree();
	}

}//UserExec
//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::RunAODhist()
{
  
  //input event
  AliAODEvent *aod = (AliAODEvent*) InputEvent();
  if(!aod) return;
  
  //Trigger
  TString trigger = aod->GetFiredTriggerClasses();
  
  if( !trigger.Contains("CCUP") ) return;
  
  fRunNum = aod ->GetRunNumber();
  fHistTriggersPerRun->Fill(fRunNum);

  //primary vertex
  AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
  fVtxContrib = fAODVertex->GetNContributors();
  if(fVtxContrib < 2) return;


  //VZERO, ZDC
  AliAODVZERO *fV0data = aod ->GetVZEROData();
  
  fV0Adecision = fV0data->GetV0ADecision();
  fV0Cdecision = fV0data->GetV0CDecision();
  if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;

  AliAODZDC *fZDCdata = aod->GetZDCData();
  
  fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
  fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
  if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
  
  Int_t nGoodV0s=0;
  Int_t V0Index[3] = {-1,-1,-1};
  Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};

  Int_t nGoodTracks=0;
  Int_t TrackID[5] = {-1,-1,-1,-1,-1};


  //V0s loop
  for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
    AliAODv0 *v0 = aod->GetV0(iV0);
    if( !v0 ) continue;
    Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
    if (lOnFlyStatus) continue;
    
    AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
    AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
    if (!pTrack || !nTrack) continue;

    if ( pTrack->Charge() == nTrack->Charge())continue;

      if(!(pTrack->TestFilterBit(1<<0))) continue;
      if(!(nTrack->TestFilterBit(1<<0))) continue;
      if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(pTrack->GetTPCNcls() < 50)continue;
      if(nTrack->GetTPCNcls() < 50)continue;
      if(pTrack->Chi2perNDF() > 4)continue;
      if(nTrack->Chi2perNDF() > 4)continue;
      
      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
      if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
      if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
      
      V0Index[nGoodV0s] = iV0;
      V0TrackID[2*nGoodV0s] = pTrack->GetID();
      V0TrackID[2*nGoodV0s+1] = nTrack->GetID();

      nGoodV0s++; 
      if(nGoodV0s > 2) break;
  }//V0s loop
  
  //Track loop
  for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
    AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
    if(!trk) AliFatal("Not a standard AOD");
    if( !trk ) continue;
    if(!(trk->TestFilterBit(1<<0))) continue;

      if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(trk->GetTPCNcls() < 50)continue;
      if(trk->Chi2perNDF() > 4)continue;
      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
      if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
     
      TrackID[nGoodTracks] = trk->GetID();
      nGoodTracks++;
				  
      if(nGoodTracks > 4) break;  
  }//Track loop

  if(nGoodV0s == 2 && nGoodTracks == 4){
  	SortArray(TrackID);
  	SortArray(V0TrackID);
	for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
	
	AliAODv0 *v00 = aod->GetV0(V0Index[0]);	
	AliAODv0 *v01 = aod->GetV0(V0Index[1]);				
	fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
			
  }
  
  PostData(2, fListHist);

}

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::RunAODtree()
{
  //input event
  AliAODEvent *aod = (AliAODEvent*) InputEvent();
  if(!aod) return;

  //input data
  AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
if(!header) AliFatal("Not a standard AOD");

  fDataFilnam = header->GetESDFileName();
  fEvtNum = header->GetEventNumberESDFile();
  fRunNum = aod->GetRunNumber();
  

  //Trigger
  //Trigger
  TString trigger = aod->GetFiredTriggerClasses();
  
  fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
  fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
  fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
  
  Bool_t isTriggered = kFALSE;
  for(Int_t i=0; i<ntrg; i++) {
    if( fTrigger[i] ) isTriggered = kTRUE;
  }
  if(!isTriggered ) return;

  //trigger inputs
  fL0inputs = ((AliVAODHeader*)aod->GetHeader())->GetL0TriggerInputs();
  fL1inputs = ((AliVAODHeader*)aod->GetHeader())->GetL1TriggerInputs();

  //Event identification
  fPerNum = aod ->GetPeriodNumber();
  fOrbNum = aod ->GetOrbitNumber();
  fBCrossNum = aod ->GetBunchCrossNumber();

  //primary vertex
  AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
  fVtxContrib = fAODVertex->GetNContributors();

  //Tracklets
  fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();

  //VZERO, ZDC
  AliAODVZERO *fV0data = aod ->GetVZEROData();
  AliAODZDC *fZDCdata = aod->GetZDCData();
  
  fV0Adecision = fV0data->GetV0ADecision();
  fV0Cdecision = fV0data->GetV0CDecision();
  fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
  fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
  
  Int_t nGoodV0s=0;
  Int_t V0Index[3] = {-1,-1,-1};
  Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};

  Int_t nGoodTracks=0;
  Int_t TrackID[5] = {-1,-1,-1,-1,-1};

  //V0s loop
  for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
    AliAODv0 *v0 = aod->GetV0(iV0);
    if( !v0 ) continue;
    Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
    if (lOnFlyStatus) continue;
    
    AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
    AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
    if (!pTrack || !nTrack) continue;

    if ( pTrack->Charge() == nTrack->Charge())continue;

      if(!(pTrack->TestFilterBit(1<<0))) continue;
      if(!(nTrack->TestFilterBit(1<<0))) continue;
      if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(pTrack->GetTPCNcls() < 50)continue;
      if(nTrack->GetTPCNcls() < 50)continue;
      if(pTrack->Chi2perNDF() > 4)continue;
      if(nTrack->Chi2perNDF() > 4)continue;
      
      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
      if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
      if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
      
      V0Index[nGoodV0s] = iV0;
      V0TrackID[2*nGoodV0s] = pTrack->GetID();
      V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
      nGoodV0s++; 
      if(nGoodV0s > 2) break;
  }//V0s loop
  
  //Track loop
  for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
    AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
    if(!trk) AliFatal("Not a standard AOD");
    if( !trk ) continue;
    if(!(trk->TestFilterBit(1<<0))) continue;

      if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(trk->GetTPCNcls() < 50)continue;
      if(trk->Chi2perNDF() > 4)continue;
      Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
      if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
      if(TMath::Abs(dca[1]) > 2) continue;
     
      TrackID[nGoodTracks] = trk->GetID();
      nGoodTracks++;
				  
      if(nGoodTracks > 4) break;  
  }//Track loop
    
  if(nGoodV0s == 2 && nGoodTracks == 4){
  	//SortArray(TrackID);
  	//SortArray(V0TrackID);
	//for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
  	for(Int_t i=0; i<2; i++){
	  	AliAODv0 *v0 = aod->GetV0(V0Index[i]);
		AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
    		AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
				
		new((*fK0sAODv0s)[i]) AliAODv0(*v0);
		new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
		new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
		}
  fK0sTree ->Fill();
  
  }
  
  PostData(1, fK0sTree);   

}//RunAOD

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
    for (Int_t i = 3; i > 0; --i) {
        for (Int_t j = 0; j < i; ++j) {
            if (dArray[j] > dArray[j+1]) {
                Int_t dTemp = dArray[j];
                dArray[j] = dArray[j+1];
                dArray[j+1] = dTemp;
            }
        }
    }
}

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::RunESD()
{

}//RunESD

//_____________________________________________________________________________
void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *) 
{

  cout<<"Analysis complete."<<endl;
}//Terminate
Commit	Line	Data
8b17ae4a MB	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>
86e51cbb	19	#include <bitset>
8b17ae4a MB	20
	21	// root headers
	22	#include "TH1I.h"
	23	#include "TTree.h"
	24	#include "TClonesArray.h"
	25	#include "TParticle.h"
18e1cc2e	26	#include "TString.h"
8b17ae4a MB	27	#include "TFile.h"
	28	#include "TDatabasePDG.h"
	29	#include "TLorentzVector.h"
	30
	31	// aliroot headers
	32	#include "AliAnalysisManager.h"
	33	#include "AliInputEventHandler.h"
	34	#include "AliESDEvent.h"
	35	#include "AliAODEvent.h"
	36	#include "AliMCEvent.h"
	37	#include "AliAODVZERO.h"
	38	#include "AliAODZDC.h"
	39	#include "AliESDVZERO.h"
	40	#include "AliESDZDC.h"
	41	#include "AliPIDResponse.h"
	42	#include "AliAODTrack.h"
	43	#include "AliAODPid.h"
	44	#include "AliAODVertex.h"
	45	#include "AliESDVertex.h"
	46	#include "AliMultiplicity.h"
	47	#include "AliESDtrack.h"
	48	#include "AliESDMuonTrack.h"
	49	#include "AliAODMCParticle.h"
	50	#include "AliMCParticle.h"
	51
	52	// my headers
	53	#include "AliAnalysisTaskUpcK0sK0s.h"
	54
	55	ClassImp(AliAnalysisTaskUpcK0sK0s);
	56
	57	using std::cout;
	58	using std::endl;
	59
	60	//trees for UPC analysis,
	61	// michal.broz@cern.ch
	62
	63	//_____________________________________________________________________________
	64	AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s()
489951cd	65	: AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
8b17ae4a MB	66	fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
	67	fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
	68	fDataFilnam(0),fRecoPass(0),fEvtNum(0),
	69	fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
	70	fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
	71
	72	{
	73
	74	//Dummy constructor
	75
	76	}//AliAnalysisTaskUpcK0sK0s
	77
	78
	79	//_____________________________________________________________________________
	80	AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s(const char *name)
489951cd	81	: AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
8b17ae4a MB	82	fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
	83	fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
	84	fDataFilnam(0),fRecoPass(0),fEvtNum(0),
	85	fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
	86	fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
	87
	88	{
	89
	90	// Constructor
	91	if( strstr(name,"ESD") ) fType = 0;
	92	if( strstr(name,"AOD") ) fType = 1;
	93
	94	Init();
	95
	96	DefineOutput(1, TTree::Class());
489951cd	97	DefineOutput(2, TList::Class());
8b17ae4a MB	98
	99	}//AliAnalysisTaskUpcK0sK0s
	100
	101	//_____________________________________________________________________________
	102	void AliAnalysisTaskUpcK0sK0s::Init()
	103	{
	104
	105	for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
	106
	107	}//Init
	108
	109	//_____________________________________________________________________________
	110	AliAnalysisTaskUpcK0sK0s::~AliAnalysisTaskUpcK0sK0s()
	111	{
	112	// Destructor
	113	if(fK0sTree){
	114	delete fK0sTree;
	115	fK0sTree = 0x0;
	116	}
8b17ae4a MB	117	if(fListHist){
	118	delete fListHist;
	119	fListHist = 0x0;
	120	}
	121
	122	}//~AliAnalysisTaskUpcK0sK0s
	123
	124
	125	//_____________________________________________________________________________
	126	void AliAnalysisTaskUpcK0sK0s::UserCreateOutputObjects()
	127	{
8b17ae4a	128
8b17ae4a MB	129	//vertices
	130	fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
	131	fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
	132
	133	//tracks
	134	fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
	135	fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);
	136
	137	//output tree with K0s candidate events
	138	fK0sTree = new TTree("fK0sTree", "fK0sTree");
	139	fK0sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
	140	fK0sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
	141	fK0sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
	142
	143	fK0sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
	144	fK0sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
	145	fK0sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
	146	fK0sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
	147	fK0sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
	148	fK0sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
	149	fK0sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
	150	fK0sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
	151	fK0sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
	152	fK0sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
	153	fK0sTree ->Branch("fDataFilnam", &fDataFilnam);
	154	fK0sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
	155	fK0sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
	156	if( fType == 0 ) {
	157	fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
	158	fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
	159	}
	160	if( fType == 1 ) {
	161	fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
	162	fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
	163	}
	164
	165
	166	fListHist = new TList();
	167	fListHist ->SetOwner();
	168
	169	fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
	170	fListHist->Add(fHistTriggersPerRun);
	171
489951cd	172	fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
8b17ae4a MB	173	fListHist->Add(fHistK0sMass);
	174
	175	PostData(1, fK0sTree);
489951cd	176	PostData(2, fListHist);
8b17ae4a MB	177
	178	}//UserCreateOutputObjects
	179
	180	//_____________________________________________________________________________
	181	void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *)
	182	{
	183
	184	//cout<<"#################### Next event ##################"<<endl;
	185
	186	if( fType == 0 ) RunESD();
	187	if( fType == 1 ){
	188	if(fRunHist) RunAODhist();
	189	if(fRunTree) RunAODtree();
	190	}
	191
	192	}//UserExec
	193	//_____________________________________________________________________________
	194	void AliAnalysisTaskUpcK0sK0s::RunAODhist()
	195	{
18e1cc2e	196
8b17ae4a MB	197	//input event
	198	AliAODEvent aod = (AliAODEvent) InputEvent();
	199	if(!aod) return;
86e51cbb	200
8b17ae4a MB	201	//Trigger
	202	TString trigger = aod->GetFiredTriggerClasses();
	203
8a2486ba	204	if( !trigger.Contains("CCUP") ) return;
8b17ae4a MB	205
	206	fRunNum = aod ->GetRunNumber();
	207	fHistTriggersPerRun->Fill(fRunNum);
	208
8b17ae4a MB	209	//primary vertex
	210	AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
	211	fVtxContrib = fAODVertex->GetNContributors();
	212	if(fVtxContrib < 2) return;
	213
	214
	215	//VZERO, ZDC
	216	AliAODVZERO *fV0data = aod ->GetVZEROData();
	217
	218	fV0Adecision = fV0data->GetV0ADecision();
	219	fV0Cdecision = fV0data->GetV0CDecision();
	220	if(fV0Adecision != AliAODVZERO::kV0Empty \|\| fV0Cdecision != AliAODVZERO::kV0Empty) return;
	221
	222	AliAODZDC *fZDCdata = aod->GetZDCData();
	223
	224	fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
	225	fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
18e1cc2e	226	if( fZDCAenergy > 8200 \|\| fZDCCenergy > 8200) return;
8b17ae4a MB	227
	228	Int_t nGoodV0s=0;
	229	Int_t V0Index[3] = {-1,-1,-1};
	230	Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
	231
	232	Int_t nGoodTracks=0;
	233	Int_t TrackID[5] = {-1,-1,-1,-1,-1};
	234
	235
	236	//V0s loop
	237	for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
	238	AliAODv0 *v0 = aod->GetV0(iV0);
	239	if( !v0 ) continue;
	240	Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
	241	if (lOnFlyStatus) continue;
	242
	243	AliAODTrack pTrack=(AliAODTrack )v0->GetDaughter(0); //0->Positive Daughter
	244	AliAODTrack nTrack=(AliAODTrack )v0->GetDaughter(1); //1->Negative Daughter
	245	if (!pTrack \|\| !nTrack) continue;
	246
	247	if ( pTrack->Charge() == nTrack->Charge())continue;
	248
86e51cbb	249	if(!(pTrack->TestFilterBit(1<<0))) continue;
86e51cbb	250	if(!(nTrack->TestFilterBit(1<<0))) continue;
18e1cc2e	251	if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
18e1cc2e	252	if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a MB	253	if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
	254	if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
	255	if(pTrack->GetTPCNcls() < 50)continue;
	256	if(nTrack->GetTPCNcls() < 50)continue;
	257	if(pTrack->Chi2perNDF() > 4)continue;
	258	if(nTrack->Chi2perNDF() > 4)continue;
	259
	260	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
	261	if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	262	if(TMath::Abs(dca[1]) > 2) continue;
	263	if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	264	if(TMath::Abs(dca[1]) > 2) continue;
	265
	266	V0Index[nGoodV0s] = iV0;
	267	V0TrackID[2*nGoodV0s] = pTrack->GetID();
	268	V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
86e51cbb	269
8b17ae4a MB	270	nGoodV0s++;
	271	if(nGoodV0s > 2) break;
	272	}//V0s loop
	273
	274	//Track loop
	275	for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69	276	AliAODTrack trk = dynamic_cast<AliAODTrack>(aod->GetTrack(itr));
f15c1f69	277	if(!trk) AliFatal("Not a standard AOD");
8b17ae4a	278	if( !trk ) continue;
8a2486ba	279	if(!(trk->TestFilterBit(1<<0))) continue;
8b17ae4a MB	280
	281	if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
	282	if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
	283	if(trk->GetTPCNcls() < 50)continue;
	284	if(trk->Chi2perNDF() > 4)continue;
	285	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
	286	if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	287	if(TMath::Abs(dca[1]) > 2) continue;
	288
	289	TrackID[nGoodTracks] = trk->GetID();
	290	nGoodTracks++;
	291
	292	if(nGoodTracks > 4) break;
	293	}//Track loop
	294
	295	if(nGoodV0s == 2 && nGoodTracks == 4){
86e51cbb	296	SortArray(TrackID);
	297	SortArray(V0TrackID);
	298	for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
	299
489951cd	300	AliAODv0 *v00 = aod->GetV0(V0Index[0]);
489951cd	301	AliAODv0 *v01 = aod->GetV0(V0Index[1]);
18e1cc2e	302	fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
18e1cc2e	303
8b17ae4a	304	}
8b17ae4a	305
489951cd	306	PostData(2, fListHist);
8b17ae4a MB	307
	308	}
	309
	310	//_____________________________________________________________________________
	311	void AliAnalysisTaskUpcK0sK0s::RunAODtree()
	312	{
	313	//input event
	314	AliAODEvent aod = (AliAODEvent) InputEvent();
	315	if(!aod) return;
	316
	317	//input data
0a918d8d	318	AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
	319	if(!header) AliFatal("Not a standard AOD");
	320
	321	fDataFilnam = header->GetESDFileName();
	322	fEvtNum = header->GetEventNumberESDFile();
18e1cc2e	323	fRunNum = aod->GetRunNumber();
86e51cbb	324
8b17ae4a	325
86e51cbb	326	//Trigger
8b17ae4a MB	327	//Trigger
	328	TString trigger = aod->GetFiredTriggerClasses();
	329
86e51cbb	330	fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
	331	fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
	332	fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
	333
	334	Bool_t isTriggered = kFALSE;
	335	for(Int_t i=0; i<ntrg; i++) {
	336	if( fTrigger[i] ) isTriggered = kTRUE;
	337	}
	338	if(!isTriggered ) return;
8b17ae4a	339
8b17ae4a	340	//trigger inputs
0a918d8d	341	fL0inputs = ((AliVAODHeader*)aod->GetHeader())->GetL0TriggerInputs();
0a918d8d	342	fL1inputs = ((AliVAODHeader*)aod->GetHeader())->GetL1TriggerInputs();
8b17ae4a MB	343
	344	//Event identification
	345	fPerNum = aod ->GetPeriodNumber();
	346	fOrbNum = aod ->GetOrbitNumber();
	347	fBCrossNum = aod ->GetBunchCrossNumber();
	348
	349	//primary vertex
	350	AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
	351	fVtxContrib = fAODVertex->GetNContributors();
	352
	353	//Tracklets
	354	fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
	355
	356	//VZERO, ZDC
	357	AliAODVZERO *fV0data = aod ->GetVZEROData();
	358	AliAODZDC *fZDCdata = aod->GetZDCData();
	359
	360	fV0Adecision = fV0data->GetV0ADecision();
	361	fV0Cdecision = fV0data->GetV0CDecision();
	362	fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
	363	fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
	364
	365	Int_t nGoodV0s=0;
	366	Int_t V0Index[3] = {-1,-1,-1};
	367	Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
	368
	369	Int_t nGoodTracks=0;
	370	Int_t TrackID[5] = {-1,-1,-1,-1,-1};
	371
	372	//V0s loop
	373	for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
	374	AliAODv0 *v0 = aod->GetV0(iV0);
	375	if( !v0 ) continue;
	376	Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
	377	if (lOnFlyStatus) continue;
	378
	379	AliAODTrack pTrack=(AliAODTrack )v0->GetDaughter(0); //0->Positive Daughter
	380	AliAODTrack nTrack=(AliAODTrack )v0->GetDaughter(1); //1->Negative Daughter
	381	if (!pTrack \|\| !nTrack) continue;
	382
	383	if ( pTrack->Charge() == nTrack->Charge())continue;
	384
86e51cbb	385	if(!(pTrack->TestFilterBit(1<<0))) continue;
86e51cbb	386	if(!(nTrack->TestFilterBit(1<<0))) continue;
8b17ae4a MB	387	if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
8b17ae4a MB	388	if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
18e1cc2e	389	if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
18e1cc2e	390	if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a MB	391	if(pTrack->GetTPCNcls() < 50)continue;
	392	if(nTrack->GetTPCNcls() < 50)continue;
	393	if(pTrack->Chi2perNDF() > 4)continue;
	394	if(nTrack->Chi2perNDF() > 4)continue;
	395
	396	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
	397	if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	398	if(TMath::Abs(dca[1]) > 2) continue;
	399	if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	400	if(TMath::Abs(dca[1]) > 2) continue;
	401
	402	V0Index[nGoodV0s] = iV0;
	403	V0TrackID[2*nGoodV0s] = pTrack->GetID();
	404	V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
	405	nGoodV0s++;
	406	if(nGoodV0s > 2) break;
	407	}//V0s loop
	408
	409	//Track loop
	410	for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69	411	AliAODTrack trk = dynamic_cast<AliAODTrack>(aod->GetTrack(itr));
f15c1f69	412	if(!trk) AliFatal("Not a standard AOD");
8b17ae4a	413	if( !trk ) continue;
8a2486ba	414	if(!(trk->TestFilterBit(1<<0))) continue;
8b17ae4a MB	415
	416	if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
	417	if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
	418	if(trk->GetTPCNcls() < 50)continue;
	419	if(trk->Chi2perNDF() > 4)continue;
	420	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
	421	if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
	422	if(TMath::Abs(dca[1]) > 2) continue;
	423
	424	TrackID[nGoodTracks] = trk->GetID();
	425	nGoodTracks++;
	426
	427	if(nGoodTracks > 4) break;
	428	}//Track loop
	429
	430	if(nGoodV0s == 2 && nGoodTracks == 4){
	431	//SortArray(TrackID);
	432	//SortArray(V0TrackID);
489951cd	433	//for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
8b17ae4a MB	434	for(Int_t i=0; i<2; i++){
	435	AliAODv0 *v0 = aod->GetV0(V0Index[i]);
	436	AliAODTrack pTrack=(AliAODTrack )v0->GetDaughter(0); //0->Positive Daughter
	437	AliAODTrack nTrack=(AliAODTrack )v0->GetDaughter(1); //1->Negative Daughter
	438
	439	new((fK0sAODv0s)[i]) AliAODv0(v0);
	440	new((fK0sAODTracks)[2i]) AliAODTrack(*pTrack);
	441	new((fK0sAODTracks)[2i+1]) AliAODTrack(*nTrack);
	442	}
	443	fK0sTree ->Fill();
	444
	445	}
	446
	447	PostData(1, fK0sTree);
8b17ae4a MB	448
	449	}//RunAOD
	450
	451	//_____________________________________________________________________________
86e51cbb	452	void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
8b17ae4a MB	453	for (Int_t i = 3; i > 0; --i) {
	454	for (Int_t j = 0; j < i; ++j) {
	455	if (dArray[j] > dArray[j+1]) {
86e51cbb	456	Int_t dTemp = dArray[j];
8b17ae4a MB	457	dArray[j] = dArray[j+1];
	458	dArray[j+1] = dTemp;
	459	}
	460	}
	461	}
	462	}
	463
	464	//_____________________________________________________________________________
	465	void AliAnalysisTaskUpcK0sK0s::RunESD()
	466	{
	467
8b17ae4a MB	468	}//RunESD
	469
	470	//_____________________________________________________________________________
	471	void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
	472	{
	473
	474	cout<<"Analysis complete."<<endl;
	475	}//Terminate
	476