[u/mrichter/AliRoot.git] / PWGUD / UPC / AliAnalysisTaskUpcPsi2s.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>

// root headers
#include "TH1I.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TParticle.h"
#include "TObjString.h"
#include "TFile.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 "AliAnalysisTaskUpcPsi2s.h"

ClassImp(AliAnalysisTaskUpcPsi2s);

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

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

//_____________________________________________________________________________
AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s() 
  : AliAnalysisTaskSE(),fType(0),hCounter(0),fJPsiTree(0),fPsi2sTree(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),
    fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0)
{

//Dummy constructor

}//AliAnalysisTaskUpcPsi2s


//_____________________________________________________________________________
AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name) 
  : AliAnalysisTaskSE(name),fType(0),hCounter(0),fJPsiTree(0),fPsi2sTree(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),
    fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0)
{

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

  DefineOutput(1, TTree::Class());
  DefineOutput(2, TTree::Class());
  DefineOutput(3, TH1I::Class());

}//AliAnalysisTaskUpcPsi2s

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

}//Init

//_____________________________________________________________________________
AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s() 
{
  // Destructor
  if(fJPsiTree){
     delete fJPsiTree;
     fJPsiTree = 0x0;
  }
  if(fPsi2sTree){
     delete fPsi2sTree;
     fPsi2sTree = 0x0;
  }
  if(hCounter){
     delete hCounter;
     hCounter = 0x0;
  }


}//~AliAnalysisTaskUpcPsi2s


//_____________________________________________________________________________
void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
{
   hCounter = new TH1I("hCounter", "hCounter", 34000, 1., 34001.);

  //input file
  fDataFilnam = new TObjString();
  fDataFilnam->SetString("");

    //tracks
  fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
  fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
  fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
  fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);

  //output tree with JPsi candidate events
  fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
  fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
  fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
  fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
  
  fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
  fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
  fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
  fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
  fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
  fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
  fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
  fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
  fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
  fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
  fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
  fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
  fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L"); 		       
  if( fType == 0 ) {
    fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
  }
  if( fType == 1 ) {
    fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
  }
 
 //output tree with Psi2s candidate events
  fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
  fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
  fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
  fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
  
  fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
  fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
  fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
  fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
  fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
  fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
  fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
  fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
  fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
  fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
  fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
  fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
  fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");  		       
  if( fType == 0 ) {
    fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
  }
  if( fType == 1 ) {
    fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
  }
  
  PostData(1, fJPsiTree);
  PostData(2, fPsi2sTree);
  PostData(3, hCounter);

}//UserCreateOutputObjects

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

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

  if( fType == 0 ) RunESD();
  if( fType == 1 ) RunAOD();

}//UserExec

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

  //input data
  const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
  fDataFilnam->Clear();
  fDataFilnam->SetString(filnam);
  fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
  fRunNum = aod ->GetRunNumber();

  hCounter->Fill( 1 );

  //Trigger
  TString trigger = aod->GetFiredTriggerClasses();
  
  fTrigger[0]   = trigger.Contains("CINT7-B");
  fTrigger[1]   = trigger.Contains("CCUP4-B"); // CE 
  fTrigger[2]   = trigger.Contains("CCUP4-E"); // CE 
  /*/
  fTrigger[0]  = trigger.Contains("CCUP6-B"); // CE
  fTrigger[1]  = trigger.Contains("CCUP6-ACE");
  fTrigger[2]  = trigger.Contains("CCUP7-B"); // CE
  fTrigger[3]  = trigger.Contains("CCUP7-ACE");/*/

  Bool_t isTRG = kFALSE;
  for(Int_t i=1; i<ntrg; i++) {
    if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
  }
  if( !isTRG ) {PostData(3, hCounter); return;}

  hCounter->Fill( 2 );

  //trigger inputs
  fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
  fL1inputs = 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 nLepton=0;
  // Int_t nPion=0;
  // Int_t nHighPt=0;
  
   Int_t nGoodTracks=0;
  //Two tracks loop
  Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};

  //Track loop
  for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
    AliAODTrack *trk = aod->GetTrack(itr);
    if( !trk ) 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;
     
      TrackIndex[nGoodTracks] = itr;
      nGoodTracks++;
				  
      if(nGoodTracks > 4) break;  
  }//Track loop
  
  if(nGoodTracks == 2){
  	  for(Int_t i=0; i<2; i++){
	  	AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
		
		Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
		trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
				
		trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
		new((*fJPsiAODTracks)[i]) AliAODTrack(*trk); 
		
  		}
  fJPsiTree ->Fill();
  PostData(1, fJPsiTree);
  }
  
  if(nGoodTracks == 4){
  	  for(Int_t i=0; i<4; i++){
	  	AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
		
		Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
		trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
		
		trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
		new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk); 
		
  		}
  fPsi2sTree ->Fill();
  PostData(2, fPsi2sTree);
  }
    
  PostData(3, hCounter);

}//RunAOD

//_____________________________________________________________________________
void AliAnalysisTaskUpcPsi2s::RunESD()
{

  //input event
  AliESDEvent *esd = (AliESDEvent*) InputEvent();
  if(!esd) return;

  //input data
  const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
  fDataFilnam->Clear();
  fDataFilnam->SetString(filnam);
  fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
  fRunNum = esd->GetRunNumber();

  hCounter->Fill( 1 );

  //Trigger
  TString trigger = esd->GetFiredTriggerClasses();
  
  fTrigger[0]   = trigger.Contains("CINT7-B");
  fTrigger[1]   = trigger.Contains("CCUP4-B"); // CE 
  fTrigger[2]   = trigger.Contains("CCUP4-E"); // CE 
 /*/
  fTrigger[0]  = trigger.Contains("CCUP6-B"); // CE
  fTrigger[1]  = trigger.Contains("CCUP6-ACE");
  fTrigger[2]  = trigger.Contains("CCUP7-B"); // CE
  fTrigger[3]  = trigger.Contains("CCUP7-ACE");/*/


  Bool_t isTRG = kFALSE;
  for(Int_t i=1; i<ntrg; i++) {
    if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
  }
  if( !isTRG ) {PostData(3, hCounter); return;}

  hCounter->Fill( 2 );

  //trigger inputs
  fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
  fL1inputs = esd->GetHeader()->GetL1TriggerInputs();

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

  //primary vertex
  AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
  fVtxContrib = fESDVertex->GetNContributors();

  //Tracklets
  fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();

  //VZERO, ZDC
  AliESDVZERO *fV0data = esd->GetVZEROData();
  AliESDZDC *fZDCdata = esd->GetESDZDC();
  
  fV0Adecision = fV0data->GetV0ADecision();
  fV0Cdecision = fV0data->GetV0CDecision();
  fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
  fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];


   //Int_t nLepton=0;
  // Int_t nPion=0;
  // Int_t nHighPt=0;
  
   Int_t nGoodTracks=0;
  //Track loop
  Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
  for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
    AliESDtrack *trk = esd->GetTrack(itr);
    if( !trk ) continue;

      if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
      if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
      if(trk->GetTPCNcls() < 50)continue;
      if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
      Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
      if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
      trk->GetImpactParameters(dca[0],dca[1]);
      if(TMath::Abs(dca[1]) > 2) continue;
      
      TrackIndex[nGoodTracks] = itr;
      nGoodTracks++;
      if(nGoodTracks > 4) break;   
  }//Track loop

  if(nGoodTracks == 2){
  	  for(Int_t i=0; i<2; i++){
	  	AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
		
		AliExternalTrackParam cParam;
      		trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
				
		new((*fJPsiESDTracks)[i]) AliESDtrack(*trk); 
		
  		}
  fJPsiTree ->Fill();
  PostData(1, fJPsiTree);
  }
  
  if(nGoodTracks == 4){
  	  for(Int_t i=0; i<4; i++){
	  	AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
		
		AliExternalTrackParam cParam;
      		trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);

		new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk); 
		
  		}
  fPsi2sTree ->Fill();
  PostData(2, fPsi2sTree);
  }
    
  PostData(3, hCounter);

}//RunESD

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

  cout<<"Analysis complete."<<endl;
}//Terminate
Commit	Line	Data
3d16cd00	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
	28	// aliroot headers
	29	#include "AliAnalysisManager.h"
	30	#include "AliInputEventHandler.h"
	31	#include "AliESDEvent.h"
	32	#include "AliAODEvent.h"
	33	#include "AliMCEvent.h"
	34	#include "AliAODVZERO.h"
	35	#include "AliAODZDC.h"
	36	#include "AliESDVZERO.h"
	37	#include "AliESDZDC.h"
	38	#include "AliPIDResponse.h"
	39	#include "AliAODTrack.h"
	40	#include "AliAODPid.h"
	41	#include "AliAODVertex.h"
	42	#include "AliESDVertex.h"
	43	#include "AliMultiplicity.h"
	44	#include "AliESDtrack.h"
	45	#include "AliESDMuonTrack.h"
	46	#include "AliAODMCParticle.h"
	47	#include "AliMCParticle.h"
	48
	49	// my headers
	50	#include "AliAnalysisTaskUpcPsi2s.h"
	51
	52	ClassImp(AliAnalysisTaskUpcPsi2s);
	53
b8f92f9d	54	using std::cout;
	55	using std::endl;
	56
3d16cd00	57	//trees for UPC analysis,
	58	// michal.broz@cern.ch
	59
	60	//_____________________________________________________________________________
	61	AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
	62	: AliAnalysisTaskSE(),fType(0),hCounter(0),fJPsiTree(0),fPsi2sTree(0),
	63	fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
	64	fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
	65	fDataFilnam(0),fRecoPass(0),fEvtNum(0),
	66	fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0)
	67	{
	68
	69	//Dummy constructor
	70
	71	}//AliAnalysisTaskUpcPsi2s
	72
	73
	74	//_____________________________________________________________________________
	75	AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
	76	: AliAnalysisTaskSE(name),fType(0),hCounter(0),fJPsiTree(0),fPsi2sTree(0),
	77	fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
	78	fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
	79	fDataFilnam(0),fRecoPass(0),fEvtNum(0),
	80	fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0)
	81	{
	82
	83	// Constructor
	84	if( strstr(name,"ESD") ) fType = 0;
	85	if( strstr(name,"AOD") ) fType = 1;
	86
	87	Init();
	88
	89	DefineOutput(1, TTree::Class());
	90	DefineOutput(2, TTree::Class());
	91	DefineOutput(3, TH1I::Class());
	92
	93	}//AliAnalysisTaskUpcPsi2s
	94
	95	//_____________________________________________________________________________
	96	void AliAnalysisTaskUpcPsi2s::Init()
	97	{
	98
	99	for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
	100
	101	}//Init
	102
	103	//_____________________________________________________________________________
	104	AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
	105	{
	106	// Destructor
	107	if(fJPsiTree){
	108	delete fJPsiTree;
	109	fJPsiTree = 0x0;
	110	}
	111	if(fPsi2sTree){
	112	delete fPsi2sTree;
	113	fPsi2sTree = 0x0;
	114	}
	115	if(hCounter){
	116	delete hCounter;
	117	hCounter = 0x0;
	118	}
	119
	120
121	}//~AliAnalysisTaskUpcPsi2s
122
123
124	//_____________________________________________________________________________
125	void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
126	{
127	hCounter = new TH1I("hCounter", "hCounter", 34000, 1., 34001.);
128
129	//input file
130	fDataFilnam = new TObjString();
131	fDataFilnam->SetString("");
132
133	//tracks
134	fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
135	fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
136	fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
137	fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
138
139	//output tree with JPsi candidate events
140	fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
141	fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
142	fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
143	fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
144
145	fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
146	fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
147	fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
148	fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
149	fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
150	fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
151	fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
152	fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
153	fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
154	fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
155	fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
156	fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
157	fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
158	if( fType == 0 ) {
159	fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
160	}
161	if( fType == 1 ) {
162	fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
163	}
164
165	//output tree with Psi2s candidate events
166	fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
167	fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
168	fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
169	fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
170
171	fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
172	fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
173	fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
174	fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
175	fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
176	fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
177	fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
178	fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
179	fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
180	fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
181	fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
182	fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
183	fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
184	if( fType == 0 ) {
185	fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
186	}
187	if( fType == 1 ) {
188	fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
189	}
190
191	PostData(1, fJPsiTree);
192	PostData(2, fPsi2sTree);
193	PostData(3, hCounter);
194
195	}//UserCreateOutputObjects
196
197	//_____________________________________________________________________________
198	void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
199	{
200
201	//cout<<"#################### Next event ##################"<<endl;
202
203	if( fType == 0 ) RunESD();
204	if( fType == 1 ) RunAOD();
205
206	}//UserExec
207
208	//_____________________________________________________________________________
209	void AliAnalysisTaskUpcPsi2s::RunAOD()
210	{
211	//input event
212	AliAODEvent aod = (AliAODEvent) InputEvent();
213	if(!aod) return;
214
215	//input data
216	const char filnam = ((TTree) GetInputData(0))->GetCurrentFile()->GetName();
217	fDataFilnam->Clear();
218	fDataFilnam->SetString(filnam);
219	fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
220	fRunNum = aod ->GetRunNumber();
221
222	hCounter->Fill( 1 );
223
224	//Trigger
225	TString trigger = aod->GetFiredTriggerClasses();
226
227	fTrigger[0] = trigger.Contains("CINT7-B");
228	fTrigger[1] = trigger.Contains("CCUP4-B"); // CE
229	fTrigger[2] = trigger.Contains("CCUP4-E"); // CE
230	/*/
231	fTrigger[0] = trigger.Contains("CCUP6-B"); // CE
232	fTrigger[1] = trigger.Contains("CCUP6-ACE");
233	fTrigger[2] = trigger.Contains("CCUP7-B"); // CE
234	fTrigger[3] = trigger.Contains("CCUP7-ACE");/*/
235
236	Bool_t isTRG = kFALSE;
237	for(Int_t i=1; i<ntrg; i++) {
238	if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
239	}
240	if( !isTRG ) {PostData(3, hCounter); return;}
241
242	hCounter->Fill( 2 );
243
244	//trigger inputs
245	fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
246	fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
247
248	//Event identification
249	fPerNum = aod ->GetPeriodNumber();
250	fOrbNum = aod ->GetOrbitNumber();
251	fBCrossNum = aod ->GetBunchCrossNumber();
252
253	//primary vertex
254	AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
255	fVtxContrib = fAODVertex->GetNContributors();
256
257	//Tracklets
258	fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
259
260	//VZERO, ZDC
261	AliAODVZERO *fV0data = aod ->GetVZEROData();
262	AliAODZDC *fZDCdata = aod->GetZDCData();
263
264	fV0Adecision = fV0data->GetV0ADecision();
265	fV0Cdecision = fV0data->GetV0CDecision();
266	fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
267	fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
268
269	//Int_t nLepton=0;
270	// Int_t nPion=0;
271	// Int_t nHighPt=0;
272
273	Int_t nGoodTracks=0;
274	//Two tracks loop
275	Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
276
277	//Track loop
278	for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
279	AliAODTrack *trk = aod->GetTrack(itr);
280	if( !trk ) continue;
281
282	if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
283	if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
284	if(trk->GetTPCNcls() < 50)continue;
285	if(trk->Chi2perNDF() > 4)continue;
286	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
287	if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
288	if(TMath::Abs(dca[1]) > 2) continue;
289
290	TrackIndex[nGoodTracks] = itr;
291	nGoodTracks++;
292
293	if(nGoodTracks > 4) break;
294	}//Track loop
295
296	if(nGoodTracks == 2){
297	for(Int_t i=0; i<2; i++){
298	AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
299
300	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
301	trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
302
303	trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
304	new((fJPsiAODTracks)[i]) AliAODTrack(trk);
305
306	}
307	fJPsiTree ->Fill();
308	PostData(1, fJPsiTree);
309	}
310
311	if(nGoodTracks == 4){
312	for(Int_t i=0; i<4; i++){
313	AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
314
315	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
316	trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
317
318	trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZatDCA();
319	new((fPsi2sAODTracks)[i]) AliAODTrack(trk);
320
321	}
322	fPsi2sTree ->Fill();
323	PostData(2, fPsi2sTree);
324	}
325
326	PostData(3, hCounter);
327
328	}//RunAOD
329
330	//_____________________________________________________________________________
331	void AliAnalysisTaskUpcPsi2s::RunESD()
332	{
333
334	//input event
335	AliESDEvent esd = (AliESDEvent) InputEvent();
336	if(!esd) return;
337
338	//input data
339	const char filnam = ((TTree) GetInputData(0))->GetCurrentFile()->GetName();
340	fDataFilnam->Clear();
341	fDataFilnam->SetString(filnam);
342	fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
343	fRunNum = esd->GetRunNumber();
344
345	hCounter->Fill( 1 );
346
347	//Trigger
348	TString trigger = esd->GetFiredTriggerClasses();
349
350	fTrigger[0] = trigger.Contains("CINT7-B");
351	fTrigger[1] = trigger.Contains("CCUP4-B"); // CE
352	fTrigger[2] = trigger.Contains("CCUP4-E"); // CE
353	/*/
354	fTrigger[0] = trigger.Contains("CCUP6-B"); // CE
355	fTrigger[1] = trigger.Contains("CCUP6-ACE");
356	fTrigger[2] = trigger.Contains("CCUP7-B"); // CE
357	fTrigger[3] = trigger.Contains("CCUP7-ACE");/*/
358
359
360	Bool_t isTRG = kFALSE;
361	for(Int_t i=1; i<ntrg; i++) {
362	if( fTrigger[i] ) {isTRG = kTRUE; hCounter->Fill( fRunNum - 167806 + 1 + i*2000 );}
363	}
364	if( !isTRG ) {PostData(3, hCounter); return;}
365
366	hCounter->Fill( 2 );
367
368	//trigger inputs
369	fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
370	fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
371
372	//Event identification
373	fPerNum = esd->GetPeriodNumber();
374	fOrbNum = esd->GetOrbitNumber();
375	fBCrossNum = esd->GetBunchCrossNumber();
376
377	//primary vertex
378	AliESDVertex fESDVertex = (AliESDVertex) esd->GetPrimaryVertex();
379	fVtxContrib = fESDVertex->GetNContributors();
380
381	//Tracklets
382	fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
383
384	//VZERO, ZDC
385	AliESDVZERO *fV0data = esd->GetVZEROData();
386	AliESDZDC *fZDCdata = esd->GetESDZDC();
387
388	fV0Adecision = fV0data->GetV0ADecision();
389	fV0Cdecision = fV0data->GetV0CDecision();
390	fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
391	fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
392
393
394	//Int_t nLepton=0;
395	// Int_t nPion=0;
396	// Int_t nHighPt=0;
397
398	Int_t nGoodTracks=0;
399	//Track loop
400	Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
401	for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
402	AliESDtrack *trk = esd->GetTrack(itr);
403	if( !trk ) continue;
404
405	if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
406	if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
407	if(trk->GetTPCNcls() < 50)continue;
408	if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
409	Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
410	if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
411	trk->GetImpactParameters(dca[0],dca[1]);
412	if(TMath::Abs(dca[1]) > 2) continue;
413
414	TrackIndex[nGoodTracks] = itr;
415	nGoodTracks++;
416	if(nGoodTracks > 4) break;
417	}//Track loop
418
419	if(nGoodTracks == 2){
420	for(Int_t i=0; i<2; i++){
421	AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
422
423	AliExternalTrackParam cParam;
424	trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
425
426	new((fJPsiESDTracks)[i]) AliESDtrack(trk);
427
428	}
429	fJPsiTree ->Fill();
430	PostData(1, fJPsiTree);
431	}
432
433	if(nGoodTracks == 4){
434	for(Int_t i=0; i<4; i++){
435	AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
436
437	AliExternalTrackParam cParam;
438	trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
439
440	new((fPsi2sESDTracks)[i]) AliESDtrack(trk);
441
442	}
443	fPsi2sTree ->Fill();
444	PostData(2, fPsi2sTree);
445	}
446
447	PostData(3, hCounter);
448
449	}//RunESD
450
451	//_____________________________________________________________________________
452	void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
453	{
454
455	cout<<"Analysis complete."<<endl;
456	}//Terminate
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486