[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);

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