[u/mrichter/AliRoot.git] / HLT / global / physics / AliHLTCaloHistoInvMass.cxx

//-*- Mode: C++ -*-
/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        * 
 * All rights reserved.                                                   *
 *                                                                        *
 * Primary Authors: Svein Lindal                                          *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

/** 
 * @file   AliHLTCaloHistoInvMass
 * @author Svein Lindal <slindal@fys.uio.no>
 * @date 
 * @brief  Produces plots of invariant mass of two clusters. 
 */

// see header file for class documentation
// or
// refer to README to build package
// or
// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt

#include "AliHLTCaloHistoInvMass.h"
#include "AliHLTCaloClusterDataStruct.h"
#include "TObjArray.h"
#include "AliESDEvent.h"
#include "TRefArray.h"
#include "TH1F.h"
#include "TString.h"
#include "AliESDCaloCluster.h"
#include "TVector3.h"
#include "TLorentzVector.h"

AliHLTCaloHistoInvMass::AliHLTCaloHistoInvMass(TString det) :
  fHistTwoClusterInvMass0(NULL),
  fHistTwoClusterInvMass1(NULL),
  fHistTwoClusterInvMass2(NULL),
  fHistTwoClusterInvMass3(NULL),
  fHistTwoClusterInvMass4(NULL)
{
  // See header file for documentation
  fHistTwoClusterInvMass0 = new TH1F(Form("%s fHistTwoClusterInvMass0", det.Data()), Form("Invariant mass of two clusters in %s, 1.5 GeV < p_{T} < 2.0", det.Data()), 200, 0, 1);
  fHistTwoClusterInvMass0->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
  fHistTwoClusterInvMass0->GetYaxis()->SetTitle("Counts");
  fHistTwoClusterInvMass0->SetMarkerStyle(21);
  fHistArray->AddLast(fHistTwoClusterInvMass0);

  fHistTwoClusterInvMass1 = new TH1F(Form("%s fHistTwoClusterInvMass1", det.Data()), Form("Invariant mass of two clusters in %s, 2.0 GeV < p_{T} < 2.5", det.Data()), 200, 0, 1);
  fHistTwoClusterInvMass1->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
  fHistTwoClusterInvMass1->GetYaxis()->SetTitle("Counts");
  fHistTwoClusterInvMass1->SetMarkerStyle(21);
  fHistArray->AddLast(fHistTwoClusterInvMass1);

  fHistTwoClusterInvMass2 = new TH1F(Form("%s fHistTwoClusterInvMass2", det.Data()), Form("Invariant mass of two clusters in %s, 2.5 GeV < p_{T} < 3.0", det.Data()), 200, 0, 1);
  fHistTwoClusterInvMass2->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
  fHistTwoClusterInvMass2->GetYaxis()->SetTitle("Counts");
  fHistTwoClusterInvMass2->SetMarkerStyle(21);
  fHistArray->AddLast(fHistTwoClusterInvMass2);

  fHistTwoClusterInvMass3 = new TH1F(Form("%s fHistTwoClusterInvMass3", det.Data()), Form("Invariant mass of two clusters in %s, 3.0 GeV < p_{T} < 5.0", det.Data()), 200, 0, 1);
  fHistTwoClusterInvMass3->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
  fHistTwoClusterInvMass3->GetYaxis()->SetTitle("Counts");
  fHistTwoClusterInvMass3->SetMarkerStyle(21);
  fHistArray->AddLast(fHistTwoClusterInvMass3);

  fHistTwoClusterInvMass4 = new TH1F(Form("%s fHistTwoClusterInvMass4", det.Data()), Form("Invariant mass of two clusters in %s p_{T} > 5.0 GeV", det.Data()), 200, 0, 1);
  fHistTwoClusterInvMass4->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
  fHistTwoClusterInvMass4->GetYaxis()->SetTitle("Counts");
  fHistTwoClusterInvMass4->SetMarkerStyle(21);
  fHistArray->AddLast(fHistTwoClusterInvMass4);


}

AliHLTCaloHistoInvMass::~AliHLTCaloHistoInvMass()
{

  if(fHistTwoClusterInvMass0)
    delete fHistTwoClusterInvMass0;
  fHistTwoClusterInvMass0 = NULL;

  if(fHistTwoClusterInvMass1)
    delete fHistTwoClusterInvMass1;
  fHistTwoClusterInvMass1 = NULL;

  if(fHistTwoClusterInvMass2)
    delete fHistTwoClusterInvMass2;
  fHistTwoClusterInvMass2 = NULL;

  if(fHistTwoClusterInvMass3)
    delete fHistTwoClusterInvMass3;
  fHistTwoClusterInvMass3 = NULL;

  if(fHistTwoClusterInvMass4)
    delete fHistTwoClusterInvMass4;
  fHistTwoClusterInvMass4 = NULL;


}


Int_t AliHLTCaloHistoInvMass::FillHistograms(Int_t nc, vector<AliHLTCaloClusterDataStruct*> &cVec) {
  //See header file for documentation
  
  Float_t cPos[nc][3];
  Float_t cEnergy[nc];

  for(int ic = 0; ic < nc; ic++) {
    AliHLTCaloClusterDataStruct * cluster = cVec.at(ic);
    cluster->GetPosition(cPos[ic]);
    cEnergy[ic] = cluster->E();
  }

  return FillInvariantMassHistograms(nc, cPos, cEnergy);
}


Int_t AliHLTCaloHistoInvMass::FillHistograms(const Int_t nc, TRefArray * clusterArray) {
  //See header file for documentation
  
  Float_t cPos[nc][3];
  Float_t cEnergy[nc];

  for(int ic = 0; ic < nc; ic++) {
    AliESDCaloCluster * cluster = static_cast<AliESDCaloCluster*>(clusterArray->At(ic));
    cluster->GetPosition(cPos[ic]);
    cEnergy[ic] = cluster->E();
  }

  return FillInvariantMassHistograms(nc, cPos, cEnergy);

}


Int_t AliHLTCaloHistoInvMass::FillInvariantMassHistograms(const Int_t nc, Float_t cPos[][3], Float_t  cEnergy[]){ 

  Int_t iResult = 0;

  for(Int_t ic = 0; ic<(nc-1); ic++) { 
     
    //BALLE hardcoded variable
    if(cEnergy[ic] < 0.5)
      continue;

    //Get momentum vector
    TVector3 iVec(cPos[ic]);
    iVec = iVec.Unit();
    iVec = cEnergy[ic] * iVec;
    
    
    for(Int_t jc = ic+1; jc<nc; jc++) { 
     
    //BALLE hardcoded variable
      if(cEnergy[jc] < 0.5)
	continue;

      //Get second momentum vector
      TVector3 jVec(cPos[jc]);
      jVec = jVec.Unit();
      jVec = cEnergy[jc] * jVec;
      
      //Calculate inv mass
      Double_t m2 = 2 *(cEnergy[ic]* cEnergy[jc] - iVec.Dot(jVec));
      Double_t m = 0.0;
      if(m2 > 0) { 
	m = TMath::Sqrt(m2);
      }
      
      //Fill histograms
      
      Float_t sum = cEnergy[ic]+cEnergy[jc];
      if(sum > 1.5)
      {
	 if(sum > 2.0)
	 {
	    if(sum > 2.5)
	    {
	       if(sum > 3.0)
	       {
		  if(sum > 5.0)
		  {
		     fHistTwoClusterInvMass4->Fill(m);
		  }
		  else
		  {
		     fHistTwoClusterInvMass3->Fill(m);
		  }	
	       }
	       else
	       {
		  fHistTwoClusterInvMass2->Fill(m);
	       }
	    }
	    else
	    {
	       fHistTwoClusterInvMass1->Fill(m);
	    }
	 }
	 else
	 {
	    fHistTwoClusterInvMass0->Fill(m);
	 }
      }
    }
  }

  return iResult;

}

// template <class T>
// Int_t AliHLTCaloHistoInvMass::FillHistograms(Int_t nc, vector<T*> clusterVec ) {
//   Float_t cPos[nc][3];
//   Float_t cEnergy[nc];
  
//   for(int ic = 0; ic < nc; ic++) {
//     T * cluster = cVec.at(ic);
//     cluster->GetPosition(cPos[ic]);
//     cEnergy[ic] = cluster->E();
//   }
  

//   for(Int_t ic = 0; ic<(nc-1); ic++) { 

//     //Get momentum vector
//     TVector3 iVec(cPos[ic]);
//     iVec = iVec.Unit();
//     iVec = cEnergy[ic] * iVec;

    
//     for(Int_t jc = ic+1; jc<nc; jc++) { 

//       //Get second momentum vector
//       TVector3 jVec(cPos[jc]);
//       jVec = jVec.Unit();
//       jVec = cEnergy[jc] * jVec;

//       //Calculate inv mass
//       Double_t m = TMath::Sqrt( 2 *(cEnergy[ic]* cEnergy[jc] - iVec.Dot(jVec) ) );

//       //Fill histogram
//       fHistTwoClusterInvMass->Fill(m);
//     }
//   }
// }
Commit	Line	Data
766aafea	1	//-- Mode: C++ --
	2	/**************************************************************************
	3	* This file is property of and copyright by the ALICE HLT Project *
	4	* All rights reserved. *
	5	* *
ce063697	6	* Primary Authors: Svein Lindal *
766aafea	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	/**
	18	* @file AliHLTCaloHistoInvMass
ce063697	19	* @author Svein Lindal <slindal@fys.uio.no>
766aafea	20	* @date
5f9dec50	21	* @brief Produces plots of invariant mass of two clusters.
766aafea	22	*/
	23
	24	// see header file for class documentation
	25	// or
	26	// refer to README to build package
	27	// or
	28	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	29
	30	#include "AliHLTCaloHistoInvMass.h"
ce063697	31	#include "AliHLTCaloClusterDataStruct.h"
766aafea	32	#include "TObjArray.h"
766aafea	33	#include "AliESDEvent.h"
766aafea	34	#include "TRefArray.h"
	35	#include "TH1F.h"
	36	#include "TString.h"
ce063697	37	#include "AliESDCaloCluster.h"
1a530a79	38	#include "TVector3.h"
1a530a79	39	#include "TLorentzVector.h"
766aafea	40
766aafea	41	AliHLTCaloHistoInvMass::AliHLTCaloHistoInvMass(TString det) :
f00fcfff	42	fHistTwoClusterInvMass0(NULL),
	43	fHistTwoClusterInvMass1(NULL),
	44	fHistTwoClusterInvMass2(NULL),
	45	fHistTwoClusterInvMass3(NULL),
	46	fHistTwoClusterInvMass4(NULL)
766aafea	47	{
766aafea	48	// See header file for documentation
15c22d47	49	fHistTwoClusterInvMass0 = new TH1F(Form("%s fHistTwoClusterInvMass0", det.Data()), Form("Invariant mass of two clusters in %s, 1.5 GeV < p_{T} < 2.0", det.Data()), 200, 0, 1);
f00fcfff	50	fHistTwoClusterInvMass0->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
	51	fHistTwoClusterInvMass0->GetYaxis()->SetTitle("Counts");
	52	fHistTwoClusterInvMass0->SetMarkerStyle(21);
	53	fHistArray->AddLast(fHistTwoClusterInvMass0);
	54
15c22d47	55	fHistTwoClusterInvMass1 = new TH1F(Form("%s fHistTwoClusterInvMass1", det.Data()), Form("Invariant mass of two clusters in %s, 2.0 GeV < p_{T} < 2.5", det.Data()), 200, 0, 1);
f00fcfff	56	fHistTwoClusterInvMass1->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
	57	fHistTwoClusterInvMass1->GetYaxis()->SetTitle("Counts");
	58	fHistTwoClusterInvMass1->SetMarkerStyle(21);
	59	fHistArray->AddLast(fHistTwoClusterInvMass1);
	60
15c22d47	61	fHistTwoClusterInvMass2 = new TH1F(Form("%s fHistTwoClusterInvMass2", det.Data()), Form("Invariant mass of two clusters in %s, 2.5 GeV < p_{T} < 3.0", det.Data()), 200, 0, 1);
6a6e8a01	62	fHistTwoClusterInvMass2->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
	63	fHistTwoClusterInvMass2->GetYaxis()->SetTitle("Counts");
	64	fHistTwoClusterInvMass2->SetMarkerStyle(21);
	65	fHistArray->AddLast(fHistTwoClusterInvMass2);
	66
15c22d47	67	fHistTwoClusterInvMass3 = new TH1F(Form("%s fHistTwoClusterInvMass3", det.Data()), Form("Invariant mass of two clusters in %s, 3.0 GeV < p_{T} < 5.0", det.Data()), 200, 0, 1);
f00fcfff	68	fHistTwoClusterInvMass3->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
	69	fHistTwoClusterInvMass3->GetYaxis()->SetTitle("Counts");
	70	fHistTwoClusterInvMass3->SetMarkerStyle(21);
	71	fHistArray->AddLast(fHistTwoClusterInvMass3);
	72
15c22d47	73	fHistTwoClusterInvMass4 = new TH1F(Form("%s fHistTwoClusterInvMass4", det.Data()), Form("Invariant mass of two clusters in %s p_{T} > 5.0 GeV", det.Data()), 200, 0, 1);
f00fcfff	74	fHistTwoClusterInvMass4->GetXaxis()->SetTitle("m_{#gamma#gamma} GeV");
	75	fHistTwoClusterInvMass4->GetYaxis()->SetTitle("Counts");
	76	fHistTwoClusterInvMass4->SetMarkerStyle(21);
	77	fHistArray->AddLast(fHistTwoClusterInvMass4);
	78
6a6e8a01	79
766aafea	80	}
	81
	82	AliHLTCaloHistoInvMass::~AliHLTCaloHistoInvMass()
	83	{
6a6e8a01	84
f00fcfff	85	if(fHistTwoClusterInvMass0)
	86	delete fHistTwoClusterInvMass0;
	87	fHistTwoClusterInvMass0 = NULL;
	88
	89	if(fHistTwoClusterInvMass1)
	90	delete fHistTwoClusterInvMass1;
	91	fHistTwoClusterInvMass1 = NULL;
6a6e8a01	92
	93	if(fHistTwoClusterInvMass2)
	94	delete fHistTwoClusterInvMass2;
	95	fHistTwoClusterInvMass2 = NULL;
	96
f00fcfff	97	if(fHistTwoClusterInvMass3)
	98	delete fHistTwoClusterInvMass3;
	99	fHistTwoClusterInvMass3 = NULL;
	100
	101	if(fHistTwoClusterInvMass4)
	102	delete fHistTwoClusterInvMass4;
	103	fHistTwoClusterInvMass4 = NULL;
	104
6a6e8a01	105
ce063697	106	}
766aafea	107
ce063697	108
	109	Int_t AliHLTCaloHistoInvMass::FillHistograms(Int_t nc, vector<AliHLTCaloClusterDataStruct*> &cVec) {
	110	//See header file for documentation
	111
	112	Float_t cPos[nc][3];
	113	Float_t cEnergy[nc];
	114
	115	for(int ic = 0; ic < nc; ic++) {
	116	AliHLTCaloClusterDataStruct * cluster = cVec.at(ic);
	117	cluster->GetPosition(cPos[ic]);
	118	cEnergy[ic] = cluster->E();
	119	}
	120
6a6e8a01	121	return FillInvariantMassHistograms(nc, cPos, cEnergy);
766aafea	122	}
766aafea	123
6a6e8a01	124
6a6e8a01	125
695f6531	126	Int_t AliHLTCaloHistoInvMass::FillHistograms(const Int_t nc, TRefArray * clusterArray) {
ce063697	127	//See header file for documentation
	128
	129	Float_t cPos[nc][3];
	130	Float_t cEnergy[nc];
	131
	132	for(int ic = 0; ic < nc; ic++) {
	133	AliESDCaloCluster * cluster = static_cast<AliESDCaloCluster*>(clusterArray->At(ic));
	134	cluster->GetPosition(cPos[ic]);
	135	cEnergy[ic] = cluster->E();
	136	}
	137
6a6e8a01	138	return FillInvariantMassHistograms(nc, cPos, cEnergy);
	139
	140	}
	141
	142
695f6531	143	Int_t AliHLTCaloHistoInvMass::FillInvariantMassHistograms(const Int_t nc, Float_t cPos[][3], Float_t cEnergy[]){
6a6e8a01	144
	145	Int_t iResult = 0;
	146
1a530a79	147	for(Int_t ic = 0; ic<(nc-1); ic++) {
6a6e8a01	148
6a6e8a01	149	//BALLE hardcoded variable
b54ade95	150	if(cEnergy[ic] < 0.5)
6a6e8a01	151	continue;
4ec58661	152
4ec58661	153	//Get momentum vector
1a530a79	154	TVector3 iVec(cPos[ic]);
1a530a79	155	iVec = iVec.Unit();
1a530a79	156	iVec = cEnergy[ic] * iVec;
6a6e8a01	157
1a530a79	158
1a530a79	159	for(Int_t jc = ic+1; jc<nc; jc++) {
6a6e8a01	160
6a6e8a01	161	//BALLE hardcoded variable
b54ade95	162	if(cEnergy[jc] < 0.5)
6a6e8a01	163	continue;
1a530a79	164
4ec58661	165	//Get second momentum vector
4ec58661	166	TVector3 jVec(cPos[jc]);
1a530a79	167	jVec = jVec.Unit();
1a530a79	168	jVec = cEnergy[jc] * jVec;
6a6e8a01	169
4ec58661	170	//Calculate inv mass
695f6531	171	Double_t m2 = 2 (cEnergy[ic] cEnergy[jc] - iVec.Dot(jVec));
	172	Double_t m = 0.0;
	173	if(m2 > 0) {
	174	m = TMath::Sqrt(m2);
	175	}
6a6e8a01	176
f00fcfff	177	//Fill histograms
6a6e8a01	178
695f6531	179	Float_t sum = cEnergy[ic]+cEnergy[jc];
b54ade95	180	if(sum > 1.5)
f00fcfff	181	{
b54ade95	182	if(sum > 2.0)
f00fcfff	183	{
b54ade95	184	if(sum > 2.5)
f00fcfff	185	{
b54ade95	186	if(sum > 3.0)
f00fcfff	187	{
b54ade95	188	if(sum > 5.0)
	189	{
	190	fHistTwoClusterInvMass4->Fill(m);
	191	}
	192	else
	193	{
	194	fHistTwoClusterInvMass3->Fill(m);
	195	}
f00fcfff	196	}
	197	else
	198	{
b54ade95	199	fHistTwoClusterInvMass2->Fill(m);
f00fcfff	200	}
	201	}
	202	else
	203	{
b54ade95	204	fHistTwoClusterInvMass1->Fill(m);
f00fcfff	205	}
	206	}
	207	else
	208	{
b54ade95	209	fHistTwoClusterInvMass0->Fill(m);
f00fcfff	210	}
f00fcfff	211	}
ce063697	212	}
	213	}
	214
6a6e8a01	215	return iResult;
7b5c97a1	216
766aafea	217	}
6a6e8a01	218
	219	// template <class T>
	220	// Int_t AliHLTCaloHistoInvMass::FillHistograms(Int_t nc, vector<T*> clusterVec ) {
	221	// Float_t cPos[nc][3];
	222	// Float_t cEnergy[nc];
	223
	224	// for(int ic = 0; ic < nc; ic++) {
	225	// T * cluster = cVec.at(ic);
	226	// cluster->GetPosition(cPos[ic]);
	227	// cEnergy[ic] = cluster->E();
	228	// }
	229
	230
	231	// for(Int_t ic = 0; ic<(nc-1); ic++) {
	232
	233	// //Get momentum vector
	234	// TVector3 iVec(cPos[ic]);
	235	// iVec = iVec.Unit();
	236	// iVec = cEnergy[ic] * iVec;
	237
	238
	239	// for(Int_t jc = ic+1; jc<nc; jc++) {
	240
	241	// //Get second momentum vector
	242	// TVector3 jVec(cPos[jc]);
	243	// jVec = jVec.Unit();
	244	// jVec = cEnergy[jc] * jVec;
	245
	246	// //Calculate inv mass
	247	// Double_t m = TMath::Sqrt( 2 (cEnergy[ic] cEnergy[jc] - iVec.Dot(jVec) ) );
	248
	249	// //Fill histogram
	250	// fHistTwoClusterInvMass->Fill(m);
	251	// }
	252	// }
	253	// }