[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / HMPID / AliHMPIDAnalysisTask.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/

//==============================================================================
// AliHMPIDAnalysysTask - Class representing a basic analysis tool of HMPID data
// A set of histograms is created.
//==============================================================================
//
// By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
// to perform the analysis on local data, on data on alien using local machine
// and on CAF.

#ifndef AliHMPIDAnalysisTASK_CXX
#define AliHMPIDAnalysisTASK_CXX


#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "AliPIDResponse.h"
#include "AliPID.h"
#include "AliVEvent.h"
#include "AliVParticle.h"
#include "AliVTrack.h"
#include "AliESDtrackCuts.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliPID.h"
#include "AliLog.h"
#include "AliHMPIDAnalysisTask.h"

ClassImp(AliHMPIDAnalysisTask)

//__________________________________________________________________________
AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
  fESD(0x0),fMC(0x0),fUseMC(kTRUE),
  fHmpHistList(0x0),
  fHmpNevents(0x0),
  fZvertex(0x0),
  fPIDResponse(0x0),
  fTrackCuts(0x0),
  fTrackFilter(0x0),
  fTree(0x0)
{
  //
  //Default ctor
  //
  for (Int_t i=0; i<69; i++) fVar[i]=0;
}

//___________________________________________________________________________
AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
  AliAnalysisTaskSE(name),
  fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  fHmpHistList(0x0),
  fHmpNevents(0x0),
  fZvertex(0x0),
  fPIDResponse(0x0),
  fTrackCuts(0x0),
  fTrackFilter(0x0),
  fTree(0x0)
{
  //
  // Constructor. Initialization of Inputs and Outputs
  //
  for (Int_t i=0; i<69; i++) fVar[i]=0;

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

//___________________________________________________________________________
AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
{
  //
  // Assignment operator
  //
  if (this!=&c) {
    AliAnalysisTaskSE::operator=(c);
    fESD             = c.fESD;
    fMC              = c.fMC;
    fUseMC           = c.fUseMC;
    fHmpHistList     = c.fHmpHistList;
    fHmpNevents      = c.fHmpNevents;
    fZvertex         = c.fZvertex;
    fPIDResponse     = c.fPIDResponse;
    fTrackCuts       = c.fTrackCuts;
    fTrackFilter     = c.fTrackFilter;
    fTree            = c.fTree;
    for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
  }
  return *this;
}

//___________________________________________________________________________
AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
  AliAnalysisTaskSE(c),
  fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
  fHmpHistList(c.fHmpHistList),
  fHmpNevents(c.fHmpNevents),
  fZvertex(c.fZvertex),
  fPIDResponse(c.fPIDResponse),  
  fTrackCuts(c.fTrackCuts),
  fTrackFilter(c.fTrackFilter),
  fTree(c.fTree)
{
  //
  // Copy Constructor
  //
  for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
}
 
//___________________________________________________________________________
AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
  //
  //destructor
  //
  Info("~AliHMPIDAnalysisTask","Calling Destructor");
  if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
}

//___________________________________________________________________________
void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
{
  // Connect ESD here

  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if (!esdH) {
    AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
  } else
    fESD = esdH->GetEvent();
    fPIDResponse = esdH->GetPIDResponse();

  if (fUseMC){
    // Connect MC
    AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!mcH) {
      AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
      fUseMC = kFALSE;
    } else
      fMC = mcH->MCEvent();
      if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
  }
  
  fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");
  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fTrackCuts->SetMinNClustersTPC(80);
  fTrackCuts->SetMaxChi2PerClusterTPC(4);
  fTrackCuts->SetMaxDCAToVertexXY(3);
  fTrackCuts->SetMaxDCAToVertexZ(3);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetRequireITSRefit(kTRUE);
  fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
  
  fTrackFilter = new AliAnalysisFilter("trackFilter");
  fTrackFilter->AddCuts(fTrackCuts);   
    
}
//***************************************************************************************************************************************************************************
void AliHMPIDAnalysisTask::UserExec(Option_t *)
{
  
  fHmpNevents->Fill(0);

  const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
  
  if(!vertex || !vertex->GetStatus() || vertex->GetNContributors()<1) {
    
    // SPD vertex
    vertex = fESD->GetPrimaryVertexSPD();
    if(!vertex) {
      PostData(1,fHmpHistList);
      PostData(2,fTree);
      return;
    }
    if(!vertex->GetStatus()){
      PostData(1,fHmpHistList);
      PostData(2,fTree);
      return;
    }  
    if(vertex->GetNContributors()<1){
      PostData(1,fHmpHistList);
      PostData(2,fTree);
      return;
    }    
  }
  
  Double_t vtxPos[3] = {999, 999, 999};
  if(vertex) vertex->GetXYZ(vtxPos);

  fHmpNevents->Fill(1);
  
  fZvertex->Fill(vtxPos[2]);
      
  if(!fPIDResponse){
      PostData(1,fHmpHistList);
      PostData(2,fTree);
      return;
   }
      
  Double_t probsHMP[5], probsTPC[5], probsTOF[5];
  Float_t  nSigmasTOF[5], nSigmasTPC[5];
      
  AliESDtrack *track = 0;

  Double_t ktol = 0.001;
  Double_t r[3], rin[3], rout[3];

  //
  // Main loop function, executed on Event basis
  //
  for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {

    track = fESD->GetTrack(iTrack);
    if(!track) continue;
    
    //if(!fTrackFilter->IsSelected(track)) continue; 
    //if(!track->GetStatus() || !AliESDtrack::kTOFout || !AliESDtrack::kTIME || track->GetTOFsignal()<12000 || track->GetTOFsignal()>200000 || track->GetIntegratedLength()<365) continue; 
     
    track->GetXYZ(r);
    track->GetInnerXYZ(rin);
    track->GetOuterXYZ(rout);

    if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
    if(track->GetHMPIDcluIdx() < 0) continue;

    Int_t q, nph;
    Float_t x, y;
    Float_t xpc, ypc, th, ph;
    track->GetHMPIDmip(x,y,q,nph);
    track->GetHMPIDtrk(xpc,ypc,th,ph);

    Int_t ITSrefit = track->IsOn(AliESDtrack::kITSrefit);
    Int_t TPCrefit = track->IsOn(AliESDtrack::kTPCrefit);
    
    if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;

    Double_t pHmp[3] = {0}, pHmp3 = 0;
    if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);

    Float_t b[2];
    Float_t bCov[3];
    track->GetImpactParameters(b,bCov);    
    
    track->GetHMPIDpid(probsHMP);
    fPIDResponse->ComputeTOFProbability(track,AliPID::kSPECIES,probsTOF);
    fPIDResponse->ComputeTPCProbability(track,AliPID::kSPECIES,probsTPC);
    
    for(Int_t ispecie=0; ispecie<AliPID::kSPECIES; ispecie++){
      nSigmasTOF[ispecie] = fPIDResponse->NumberOfSigmasTOF((AliVTrack*)track,(AliPID::EParticleType)ispecie);
      nSigmasTPC[ispecie] = fPIDResponse->NumberOfSigmasTPC((AliVTrack*)track,(AliPID::EParticleType)ispecie);      
    }
        
    fVar[0] = track->GetHMPIDcluIdx()/1000000;
    fVar[1] = pHmp3;
    fVar[2] = (Float_t)track->P();
    fVar[3] = xpc;
    fVar[4] = ypc;
    fVar[5] = x;
    fVar[6] = y;
    fVar[7] = (Float_t)track->GetHMPIDsignal();
    fVar[8] = q;
    fVar[9] = th;
    fVar[10] = ph;
    fVar[11] = (Float_t)track->GetSign();
    fVar[12] = (Float_t)nph;
    fVar[13] = (Float_t)track->GetNcls(1);
    fVar[14] = (Float_t)probsHMP[0];
    fVar[15] = (Float_t)probsHMP[1];
    fVar[16] = (Float_t)probsHMP[2];
    fVar[17] = (Float_t)probsHMP[3];
    fVar[18] = (Float_t)probsHMP[4];
    fVar[19] = (Float_t)track->GetTOFsignal();
    fVar[20] = (Float_t)track->GetKinkIndex(0);
    fVar[21] = (Float_t)track->Xv();
    fVar[22] = (Float_t)track->Yv();
    fVar[23] = (Float_t)track->Zv();
    fVar[24] = (Float_t)track->GetTPCchi2();
    fVar[25] = b[0];
    fVar[26] = b[1];
    fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
    fVar[28] = vtxPos[0];
    fVar[29] = vtxPos[1];
    fVar[30] = vtxPos[2];
    fVar[31] = (Float_t)ITSrefit;
    fVar[32] = (Float_t)TPCrefit;
    fVar[33] = (Float_t)track->Eta();
    fVar[34] = (Float_t)r[0];
    fVar[35] = (Float_t)r[1];
    fVar[36] = (Float_t)r[2];
    fVar[37] = (Float_t)rout[0];           
    fVar[38] = (Float_t)rout[1];
    fVar[39] = (Float_t)rout[2];
    fVar[40] = (Float_t)pHmp[0];
    fVar[41] = (Float_t)pHmp[1];
    fVar[42] = (Float_t)pHmp[2];
    fVar[43] = (Float_t)track->Px();
    fVar[44] = (Float_t)track->Py();
    fVar[45] = (Float_t)track->Pz();
    fVar[46] = (Float_t)track->GetTPCClusterInfo(2,1);
    fVar[47] = (Float_t)track->GetTPCNclsF();
    fVar[48] = nSigmasTOF[0];
    fVar[49] = nSigmasTOF[1];
    fVar[50] = nSigmasTOF[2];
    fVar[51] = nSigmasTOF[3];
    fVar[52] = nSigmasTOF[4];
    fVar[53] = nSigmasTPC[0];
    fVar[54] = nSigmasTPC[1];
    fVar[55] = nSigmasTPC[2];
    fVar[56] = nSigmasTPC[3];
    fVar[57] = nSigmasTPC[4];
    fVar[58] = (Float_t)probsTOF[0];
    fVar[59] = (Float_t)probsTOF[1];
    fVar[60] = (Float_t)probsTOF[2];
    fVar[61] = (Float_t)probsTOF[3];
    fVar[62] = (Float_t)probsTOF[4];
    fVar[63] = (Float_t)probsTPC[0];
    fVar[64] = (Float_t)probsTPC[1];
    fVar[65] = (Float_t)probsTPC[2];
    fVar[66] = (Float_t)probsTPC[3];
    fVar[67] = (Float_t)probsTPC[4];
    fVar[68] = (Float_t)track->GetHMPIDchi2();
       
    fTree->Fill();
   
  }//track loop

  /* PostData(0) is taken care of by AliAnalysisTaskSE */
  PostData(1,fHmpHistList);
  PostData(2,fTree);
}
//___________________________________________________________________________
void AliHMPIDAnalysisTask::Terminate(Option_t*)
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

  Info("Terminate"," ");

  if (!fUseMC) return;

  fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));

  if (!fHmpHistList) {
    AliError("Histogram List is not available");
    return;
  }


  AliAnalysisTaskSE::Terminate();

}
//___________________________________________________________________________
void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
  //
  //HERE ONE CAN CREATE OUTPUT OBJECTS
  //TO BE SET BEFORE THE EXECUTION OF THE TASK
  //

  //slot #1
//   OpenFile(1);
   fHmpHistList = new TList();
   fHmpHistList->SetOwner();

   fHmpNevents = new TH1F("fHmpNevents","Number of events",2,0,2);
   fHmpHistList->Add(fHmpNevents);

   fZvertex = new TH1F("fZvertex","Z primary vertex distribution",4000,-20,20);
   fHmpHistList->Add(fZvertex);
   
//   OpenFile(2);
   fTree = new TTree("Tree","Tree with data");
   fTree->Branch("Chamber",&fVar[0]);
   fTree->Branch("pHmp3",&fVar[1]);
   fTree->Branch("P",&fVar[2]);
   fTree->Branch("Xpc",&fVar[3]);
   fTree->Branch("Ypc",&fVar[4]);
   fTree->Branch("X",&fVar[5]);
   fTree->Branch("Y",&fVar[6]);
   fTree->Branch("HMPIDsignal",&fVar[7]);
   fTree->Branch("Charge",&fVar[8]);
   fTree->Branch("Theta",&fVar[9]);
   fTree->Branch("Phi",&fVar[10]);
   fTree->Branch("Sign",&fVar[11]);
   fTree->Branch("NumPhotons",&fVar[12]);
   fTree->Branch("NumTPCclust",&fVar[13]);
   fTree->Branch("ProbHMP0",&fVar[14]);
   fTree->Branch("ProbHMP1",&fVar[15]);
   fTree->Branch("ProbHMP2",&fVar[16]);
   fTree->Branch("ProbHMP3",&fVar[17]);
   fTree->Branch("ProbHMP4",&fVar[18]);
   fTree->Branch("TOFsignal",&fVar[19]);
   fTree->Branch("KinkIndex",&fVar[20]);
   fTree->Branch("Xv",&fVar[21]);
   fTree->Branch("Yv",&fVar[22]);
   fTree->Branch("Zv",&fVar[23]);
   fTree->Branch("TPCchi2",&fVar[24]);
   fTree->Branch("b0",&fVar[25]);
   fTree->Branch("b1",&fVar[26]);
   fTree->Branch("ClustSize",&fVar[27]);
   fTree->Branch("PrimVertexX",&fVar[28]);
   fTree->Branch("PrimVertexY",&fVar[29]);
   fTree->Branch("PrimVertexZ",&fVar[30]);
   fTree->Branch("ITSrefit",&fVar[31]);
   fTree->Branch("TPCrefit",&fVar[32]);   
   fTree->Branch("Eta",&fVar[33]);
   fTree->Branch("xTrack",&fVar[34]);
   fTree->Branch("yTrack",&fVar[35]);
   fTree->Branch("zTrack",&fVar[36]);
   fTree->Branch("xOuterTrack",&fVar[37]);
   fTree->Branch("yOuterTrack",&fVar[38]);
   fTree->Branch("zOuterTrack",&fVar[39]);
   fTree->Branch("pHmpX",&fVar[40]);
   fTree->Branch("pHmpY",&fVar[41]);
   fTree->Branch("pHmpZ",&fVar[42]);
   fTree->Branch("Px",&fVar[43]);
   fTree->Branch("Py",&fVar[44]);
   fTree->Branch("Pz",&fVar[45]);
   fTree->Branch("nCrossedRowsTPC",&fVar[46]);
   fTree->Branch("findableClustTPC",&fVar[47]);
   fTree->Branch("NSigmasTOF0",&fVar[48]);   
   fTree->Branch("NSigmasTOF1",&fVar[49]);   
   fTree->Branch("NSigmasTOF2",&fVar[50]);   
   fTree->Branch("NSigmasTOF3",&fVar[51]);   
   fTree->Branch("NSigmasTOF4",&fVar[52]);   
   fTree->Branch("NSigmasTPC0",&fVar[53]);   
   fTree->Branch("NSigmasTPC1",&fVar[54]);   
   fTree->Branch("NSigmasTPC2",&fVar[55]);   
   fTree->Branch("NSigmasTPC3",&fVar[56]);   
   fTree->Branch("NSigmasTPC4",&fVar[57]);   
   fTree->Branch("ProbTOF0",&fVar[58]);
   fTree->Branch("ProbTOF1",&fVar[59]);
   fTree->Branch("ProbTOF2",&fVar[60]);
   fTree->Branch("ProbTOF3",&fVar[61]);
   fTree->Branch("ProbTOF4",&fVar[62]);
   fTree->Branch("ProbTPC0",&fVar[63]);
   fTree->Branch("ProbTPC1",&fVar[64]);
   fTree->Branch("ProbTPC2",&fVar[65]);
   fTree->Branch("ProbTPC3",&fVar[66]);
   fTree->Branch("ProbTPC4",&fVar[67]);
   fTree->Branch("HmpSigma",&fVar[68]);
   
   PostData(1,fHmpHistList);
   PostData(2,fTree);
}

//____________________________________________________________________________________________________________________________________
Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
{
 return abs(x - y) <= tolerance ;
}
   
#endif
Commit	Line	Data
a65a7e70	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	//==============================================================================
	17	// AliHMPIDAnalysysTask - Class representing a basic analysis tool of HMPID data
	18	// A set of histograms is created.
	19	//==============================================================================
	20	//
	21	// By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
	22	// to perform the analysis on local data, on data on alien using local machine
	23	// and on CAF.
	24
	25	#ifndef AliHMPIDAnalysisTASK_CXX
	26	#define AliHMPIDAnalysisTASK_CXX
	27
	28
	29	#include "TH1.h"
	30	#include "TH2.h"
	31	#include "TFile.h"
	32	#include "TCanvas.h"
	33	#include "TGraphErrors.h"
	34	#include "AliPIDResponse.h"
	35	#include "AliPID.h"
	36	#include "AliVEvent.h"
	37	#include "AliVParticle.h"
	38	#include "AliVTrack.h"
	39	#include "AliESDtrackCuts.h"
	40	#include "AliAnalysisFilter.h"
	41	#include "AliAnalysisManager.h"
	42	#include "AliESDInputHandler.h"
	43	#include "AliMCEventHandler.h"
	44	#include "AliMCEvent.h"
	45	#include "AliESDtrack.h"
	46	#include "AliPID.h"
	47	#include "AliLog.h"
	48	#include "AliHMPIDAnalysisTask.h"
	49
	50	ClassImp(AliHMPIDAnalysisTask)
	51
	52	//__________________________________________________________________________
	53	AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
	54	fESD(0x0),fMC(0x0),fUseMC(kTRUE),
	55	fHmpHistList(0x0),
	56	fHmpNevents(0x0),
	57	fZvertex(0x0),
	58	fPIDResponse(0x0),
	59	fTrackCuts(0x0),
	60	fTrackFilter(0x0),
	61	fTree(0x0)
	62	{
	63	//
	64	//Default ctor
65	//
66	for (Int_t i=0; i<69; i++) fVar[i]=0;
67	}
68
69	//___________________________________________________________________________
70	AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
71	AliAnalysisTaskSE(name),
72	fESD(0x0), fMC(0x0), fUseMC(kTRUE),
73	fHmpHistList(0x0),
74	fHmpNevents(0x0),
75	fZvertex(0x0),
76	fPIDResponse(0x0),
77	fTrackCuts(0x0),
78	fTrackFilter(0x0),
79	fTree(0x0)
80	{
81	//
82	// Constructor. Initialization of Inputs and Outputs
83	//
84	for (Int_t i=0; i<69; i++) fVar[i]=0;
85
86	DefineOutput(1,TList::Class());
87	DefineOutput(2,TTree::Class());
88	}
89
90	//___________________________________________________________________________
91	AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
92	{
93	//
94	// Assignment operator
95	//
96	if (this!=&c) {
97	AliAnalysisTaskSE::operator=(c);
98	fESD = c.fESD;
99	fMC = c.fMC;
100	fUseMC = c.fUseMC;
101	fHmpHistList = c.fHmpHistList;
102	fHmpNevents = c.fHmpNevents;
103	fZvertex = c.fZvertex;
104	fPIDResponse = c.fPIDResponse;
105	fTrackCuts = c.fTrackCuts;
106	fTrackFilter = c.fTrackFilter;
107	fTree = c.fTree;
108	for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
109	}
110	return *this;
111	}
112
113	//___________________________________________________________________________
114	AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
115	AliAnalysisTaskSE(c),
116	fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
117	fHmpHistList(c.fHmpHistList),
118	fHmpNevents(c.fHmpNevents),
119	fZvertex(c.fZvertex),
120	fPIDResponse(c.fPIDResponse),
121	fTrackCuts(c.fTrackCuts),
122	fTrackFilter(c.fTrackFilter),
123	fTree(c.fTree)
124	{
125	//
126	// Copy Constructor
127	//
128	for (Int_t i=0; i<69; i++) fVar[i]=c.fVar[i];
129	}
130
131	//___________________________________________________________________________
132	AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
133	//
134	//destructor
135	//
136	Info("~AliHMPIDAnalysisTask","Calling Destructor");
137	if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
138	}
139
140	//___________________________________________________________________________
141	void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
142	{
143	// Connect ESD here
144
145	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
146	if (!esdH) {
147	AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
148	} else
149	fESD = esdH->GetEvent();
150	fPIDResponse = esdH->GetPIDResponse();
151
152	if (fUseMC){
153	// Connect MC
154	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
155	if (!mcH) {
156	AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
157	fUseMC = kFALSE;
158	} else
159	fMC = mcH->MCEvent();
160	if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
161	}
162
163	fTrackCuts = new AliESDtrackCuts("fTrackCuts", "Standard");
164	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
165	fTrackCuts->SetMinNClustersTPC(80);
166	fTrackCuts->SetMaxChi2PerClusterTPC(4);
167	fTrackCuts->SetMaxDCAToVertexXY(3);
168	fTrackCuts->SetMaxDCAToVertexZ(3);
169	fTrackCuts->SetRequireTPCRefit(kTRUE);
170	fTrackCuts->SetRequireITSRefit(kTRUE);
171	fTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
172
173	fTrackFilter = new AliAnalysisFilter("trackFilter");
174	fTrackFilter->AddCuts(fTrackCuts);
175
176	}
177	//***************************************************************************************************************************************************************************
178	void AliHMPIDAnalysisTask::UserExec(Option_t *)
179	{
180
181	fHmpNevents->Fill(0);
182
183	const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
184
185	if(!vertex \|\| !vertex->GetStatus() \|\| vertex->GetNContributors()<1) {
186
187	// SPD vertex
188	vertex = fESD->GetPrimaryVertexSPD();
189	if(!vertex) {
190	PostData(1,fHmpHistList);
191	PostData(2,fTree);
192	return;
193	}
194	if(!vertex->GetStatus()){
195	PostData(1,fHmpHistList);
196	PostData(2,fTree);
197	return;
198	}
199	if(vertex->GetNContributors()<1){
200	PostData(1,fHmpHistList);
201	PostData(2,fTree);
202	return;
203	}
204	}
205
206	Double_t vtxPos[3] = {999, 999, 999};
207	if(vertex) vertex->GetXYZ(vtxPos);
208
209	fHmpNevents->Fill(1);
210
211	fZvertex->Fill(vtxPos[2]);
212
213	if(!fPIDResponse){
214	PostData(1,fHmpHistList);
215	PostData(2,fTree);
216	return;
217	}
218
219	Double_t probsHMP[5], probsTPC[5], probsTOF[5];
220	Float_t nSigmasTOF[5], nSigmasTPC[5];
221
222	AliESDtrack *track = 0;
223
224	Double_t ktol = 0.001;
225	Double_t r[3], rin[3], rout[3];
226
227	//
228	// Main loop function, executed on Event basis
229	//
230	for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
231
232	track = fESD->GetTrack(iTrack);
233	if(!track) continue;
234
235	//if(!fTrackFilter->IsSelected(track)) continue;
236	//if(!track->GetStatus() \|\| !AliESDtrack::kTOFout \|\| !AliESDtrack::kTIME \|\| track->GetTOFsignal()<12000 \|\| track->GetTOFsignal()>200000 \|\| track->GetIntegratedLength()<365) continue;
237
238	track->GetXYZ(r);
239	track->GetInnerXYZ(rin);
240	track->GetOuterXYZ(rout);
241
242	if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
243	if(track->GetHMPIDcluIdx() < 0) continue;
244
245	Int_t q, nph;
246	Float_t x, y;
247	Float_t xpc, ypc, th, ph;
248	track->GetHMPIDmip(x,y,q,nph);
249	track->GetHMPIDtrk(xpc,ypc,th,ph);
250
251	Int_t ITSrefit = track->IsOn(AliESDtrack::kITSrefit);
252	Int_t TPCrefit = track->IsOn(AliESDtrack::kTPCrefit);
253
254	if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
255
256	Double_t pHmp[3] = {0}, pHmp3 = 0;
257	if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]pHmp[0]+pHmp[1]pHmp[1]+pHmp[2]*pHmp[2]);
258
259	Float_t b[2];
260	Float_t bCov[3];
261	track->GetImpactParameters(b,bCov);
262
263	track->GetHMPIDpid(probsHMP);
264	fPIDResponse->ComputeTOFProbability(track,AliPID::kSPECIES,probsTOF);
265	fPIDResponse->ComputeTPCProbability(track,AliPID::kSPECIES,probsTPC);
266
267	for(Int_t ispecie=0; ispecie<AliPID::kSPECIES; ispecie++){
268	nSigmasTOF[ispecie] = fPIDResponse->NumberOfSigmasTOF((AliVTrack*)track,(AliPID::EParticleType)ispecie);
269	nSigmasTPC[ispecie] = fPIDResponse->NumberOfSigmasTPC((AliVTrack*)track,(AliPID::EParticleType)ispecie);
270	}
271
272	fVar[0] = track->GetHMPIDcluIdx()/1000000;
273	fVar[1] = pHmp3;
274	fVar[2] = (Float_t)track->P();
275	fVar[3] = xpc;
276	fVar[4] = ypc;
277	fVar[5] = x;
278	fVar[6] = y;
279	fVar[7] = (Float_t)track->GetHMPIDsignal();
280	fVar[8] = q;
281	fVar[9] = th;
282	fVar[10] = ph;
283	fVar[11] = (Float_t)track->GetSign();
284	fVar[12] = (Float_t)nph;
285	fVar[13] = (Float_t)track->GetNcls(1);
286	fVar[14] = (Float_t)probsHMP[0];
287	fVar[15] = (Float_t)probsHMP[1];
288	fVar[16] = (Float_t)probsHMP[2];
289	fVar[17] = (Float_t)probsHMP[3];
290	fVar[18] = (Float_t)probsHMP[4];
291	fVar[19] = (Float_t)track->GetTOFsignal();
292	fVar[20] = (Float_t)track->GetKinkIndex(0);
293	fVar[21] = (Float_t)track->Xv();
294	fVar[22] = (Float_t)track->Yv();
295	fVar[23] = (Float_t)track->Zv();
296	fVar[24] = (Float_t)track->GetTPCchi2();
297	fVar[25] = b[0];
298	fVar[26] = b[1];
299	fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
300	fVar[28] = vtxPos[0];
301	fVar[29] = vtxPos[1];
302	fVar[30] = vtxPos[2];
303	fVar[31] = (Float_t)ITSrefit;
304	fVar[32] = (Float_t)TPCrefit;
305	fVar[33] = (Float_t)track->Eta();
306	fVar[34] = (Float_t)r[0];
307	fVar[35] = (Float_t)r[1];
308	fVar[36] = (Float_t)r[2];
309	fVar[37] = (Float_t)rout[0];
310	fVar[38] = (Float_t)rout[1];
311	fVar[39] = (Float_t)rout[2];
312	fVar[40] = (Float_t)pHmp[0];
313	fVar[41] = (Float_t)pHmp[1];
314	fVar[42] = (Float_t)pHmp[2];
315	fVar[43] = (Float_t)track->Px();
316	fVar[44] = (Float_t)track->Py();
317	fVar[45] = (Float_t)track->Pz();
318	fVar[46] = (Float_t)track->GetTPCClusterInfo(2,1);
319	fVar[47] = (Float_t)track->GetTPCNclsF();
320	fVar[48] = nSigmasTOF[0];
321	fVar[49] = nSigmasTOF[1];
322	fVar[50] = nSigmasTOF[2];
323	fVar[51] = nSigmasTOF[3];
324	fVar[52] = nSigmasTOF[4];
325	fVar[53] = nSigmasTPC[0];
326	fVar[54] = nSigmasTPC[1];
327	fVar[55] = nSigmasTPC[2];
328	fVar[56] = nSigmasTPC[3];
329	fVar[57] = nSigmasTPC[4];
330	fVar[58] = (Float_t)probsTOF[0];
331	fVar[59] = (Float_t)probsTOF[1];
332	fVar[60] = (Float_t)probsTOF[2];
333	fVar[61] = (Float_t)probsTOF[3];
334	fVar[62] = (Float_t)probsTOF[4];
335	fVar[63] = (Float_t)probsTPC[0];
336	fVar[64] = (Float_t)probsTPC[1];
337	fVar[65] = (Float_t)probsTPC[2];
338	fVar[66] = (Float_t)probsTPC[3];
339	fVar[67] = (Float_t)probsTPC[4];
340	fVar[68] = (Float_t)track->GetHMPIDchi2();
341
342	fTree->Fill();
343
344	}//track loop
345
346	/* PostData(0) is taken care of by AliAnalysisTaskSE */
347	PostData(1,fHmpHistList);
348	PostData(2,fTree);
349	}
350	//___________________________________________________________________________
351	void AliHMPIDAnalysisTask::Terminate(Option_t*)
352	{
353	// The Terminate() function is the last function to be called during
354	// a query. It always runs on the client, it can be used to present
355	// the results graphically or save the results to file.
356
357	Info("Terminate"," ");
358
359	if (!fUseMC) return;
360
361	fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
362
363	if (!fHmpHistList) {
364	AliError("Histogram List is not available");
365	return;
366	}
367
368
369	AliAnalysisTaskSE::Terminate();
370
371	}
372	//___________________________________________________________________________
373	void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
374	//
375	//HERE ONE CAN CREATE OUTPUT OBJECTS
376	//TO BE SET BEFORE THE EXECUTION OF THE TASK
377	//
378
379	//slot #1
380	// OpenFile(1);
381	fHmpHistList = new TList();
382	fHmpHistList->SetOwner();
383
384	fHmpNevents = new TH1F("fHmpNevents","Number of events",2,0,2);
385	fHmpHistList->Add(fHmpNevents);
386
387	fZvertex = new TH1F("fZvertex","Z primary vertex distribution",4000,-20,20);
388	fHmpHistList->Add(fZvertex);
389
390	// OpenFile(2);
391	fTree = new TTree("Tree","Tree with data");
392	fTree->Branch("Chamber",&fVar[0]);
393	fTree->Branch("pHmp3",&fVar[1]);
394	fTree->Branch("P",&fVar[2]);
395	fTree->Branch("Xpc",&fVar[3]);
396	fTree->Branch("Ypc",&fVar[4]);
397	fTree->Branch("X",&fVar[5]);
398	fTree->Branch("Y",&fVar[6]);
399	fTree->Branch("HMPIDsignal",&fVar[7]);
400	fTree->Branch("Charge",&fVar[8]);
401	fTree->Branch("Theta",&fVar[9]);
402	fTree->Branch("Phi",&fVar[10]);
403	fTree->Branch("Sign",&fVar[11]);
404	fTree->Branch("NumPhotons",&fVar[12]);
405	fTree->Branch("NumTPCclust",&fVar[13]);
406	fTree->Branch("ProbHMP0",&fVar[14]);
407	fTree->Branch("ProbHMP1",&fVar[15]);
408	fTree->Branch("ProbHMP2",&fVar[16]);
409	fTree->Branch("ProbHMP3",&fVar[17]);
410	fTree->Branch("ProbHMP4",&fVar[18]);
411	fTree->Branch("TOFsignal",&fVar[19]);
412	fTree->Branch("KinkIndex",&fVar[20]);
413	fTree->Branch("Xv",&fVar[21]);
414	fTree->Branch("Yv",&fVar[22]);
415	fTree->Branch("Zv",&fVar[23]);
416	fTree->Branch("TPCchi2",&fVar[24]);
417	fTree->Branch("b0",&fVar[25]);
418	fTree->Branch("b1",&fVar[26]);
419	fTree->Branch("ClustSize",&fVar[27]);
420	fTree->Branch("PrimVertexX",&fVar[28]);
421	fTree->Branch("PrimVertexY",&fVar[29]);
422	fTree->Branch("PrimVertexZ",&fVar[30]);
423	fTree->Branch("ITSrefit",&fVar[31]);
424	fTree->Branch("TPCrefit",&fVar[32]);
425	fTree->Branch("Eta",&fVar[33]);
426	fTree->Branch("xTrack",&fVar[34]);
427	fTree->Branch("yTrack",&fVar[35]);
428	fTree->Branch("zTrack",&fVar[36]);
429	fTree->Branch("xOuterTrack",&fVar[37]);
430	fTree->Branch("yOuterTrack",&fVar[38]);
431	fTree->Branch("zOuterTrack",&fVar[39]);
432	fTree->Branch("pHmpX",&fVar[40]);
433	fTree->Branch("pHmpY",&fVar[41]);
434	fTree->Branch("pHmpZ",&fVar[42]);
435	fTree->Branch("Px",&fVar[43]);
436	fTree->Branch("Py",&fVar[44]);
437	fTree->Branch("Pz",&fVar[45]);
438	fTree->Branch("nCrossedRowsTPC",&fVar[46]);
439	fTree->Branch("findableClustTPC",&fVar[47]);
440	fTree->Branch("NSigmasTOF0",&fVar[48]);
441	fTree->Branch("NSigmasTOF1",&fVar[49]);
442	fTree->Branch("NSigmasTOF2",&fVar[50]);
443	fTree->Branch("NSigmasTOF3",&fVar[51]);
444	fTree->Branch("NSigmasTOF4",&fVar[52]);
445	fTree->Branch("NSigmasTPC0",&fVar[53]);
446	fTree->Branch("NSigmasTPC1",&fVar[54]);
447	fTree->Branch("NSigmasTPC2",&fVar[55]);
448	fTree->Branch("NSigmasTPC3",&fVar[56]);
449	fTree->Branch("NSigmasTPC4",&fVar[57]);
450	fTree->Branch("ProbTOF0",&fVar[58]);
451	fTree->Branch("ProbTOF1",&fVar[59]);
452	fTree->Branch("ProbTOF2",&fVar[60]);
453	fTree->Branch("ProbTOF3",&fVar[61]);
454	fTree->Branch("ProbTOF4",&fVar[62]);
455	fTree->Branch("ProbTPC0",&fVar[63]);
456	fTree->Branch("ProbTPC1",&fVar[64]);
457	fTree->Branch("ProbTPC2",&fVar[65]);
458	fTree->Branch("ProbTPC3",&fVar[66]);
459	fTree->Branch("ProbTPC4",&fVar[67]);
460	fTree->Branch("HmpSigma",&fVar[68]);
461
462	PostData(1,fHmpHistList);
463	PostData(2,fTree);
464	}
465
466	//____________________________________________________________________________________________________________________________________
467	Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
468	{
469	return abs(x - y) <= tolerance ;
470	}
471
472	#endif