[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAnalysisTaskdEdxSSDQA.cxx

//SSD dEdx gain calibration task
//to run in QA train
//autor Marek Chojnacki
//Marek.Chojnacki@cern.ch
// last line
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TFile.h"

#include "AliAnalysisTaskdEdxSSDQA.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"


#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "Riostream.h"

#include "AliESDfriend.h"
#include "AliESDfriendTrack.h"
#include "AliTrackPointArray.h"
#include "AliCDBManager.h"
#include "AliGeomManager.h"
#include "TMath.h"
#include "TGeoMatrix.h"

using namespace std;
ClassImp(AliAnalysisTaskdEdxSSDQA)

//________________________________________________________________________
AliAnalysisTaskdEdxSSDQA::AliAnalysisTaskdEdxSSDQA(const char *name) 
: AliAnalysisTaskSE(name), fHist1(0),fHist2(0),fHist3(0),fHist4(0),fHist5(0),fHist6(0),fListOfHistos(0),fPcut(0.0),fdothecorrection(0)
{
  // Constructor
	// Define input and output slots here
  // Input slot #1 works with a TChain
 // fcorrections=new Float_t[1698*12];
  for(int i=0;i<1698*12;i++)
  	fcorrections[i]=1.0;
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskdEdxSSDQA::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
  	
	fListOfHistos = new TList();
	fListOfHistos->SetOwner();
	fHist1 =new TH2F("QAChargeRatio","QAChargeRatio;Module;CR",1698,-0.5,1697.5,80,-1.0,1.0);
	fListOfHistos->Add(fHist1);
	fHist2=new TH2F("QACharge","QACharge;Module;Q",1698,-0.5,1697.5,150,0,300);
	fListOfHistos->Add(fHist2);
	fHist3=new TH3F("ChargeRatiovCharge","CRvQ;module;Q;CR",1698,-0.5,1697.5,80,0,1600,50,-1.0,1.0);
	fListOfHistos->Add(fHist3);
	fHist4=new TH2F("QAChargedinchips","chargedinchips;chip;Q",1698*12,-0.5,1698*12-0.5,150,0,300);
	fListOfHistos->Add(fHist4);
	if(fdothecorrection)
	{
		fHist5=new TH2F("QAChargeCorrected","chargedCorrected;Module;QCorrected",1698,-0.5,1697.5,150,0,300);
		fListOfHistos->Add(fHist5);
	}
	fHist6=new TH2F("QNvQP",";QN;QP",160,0,1600,1600,0,1600);
	fListOfHistos->Add(fHist6);
	
	PostData(1,  fListOfHistos);
}
//________________________________________________________________________
void AliAnalysisTaskdEdxSSDQA::LocalInit() 
{
	
	Printf("end of LocalInit");
}
//______________________________________________________________________
AliAnalysisTaskdEdxSSDQA::~AliAnalysisTaskdEdxSSDQA()
{
}

//________________________________________________________________________

void AliAnalysisTaskdEdxSSDQA::UserExec(Option_t *) 
{
    // Main loop
    // Called for each event
    AliESDEvent* esd = dynamic_cast<AliESDEvent*> (InputEvent());
    if(!esd) 
    {
	Printf("ERROR: Input ESD Event not available");
	PostData(1,  fListOfHistos);
	return;
    }
    
    Int_t nTracks=esd->GetNumberOfTracks();
    
    
	if(!ESDfriend())
	{
		Printf("problem with friend");
		PostData(1,  fListOfHistos);
		return;
	}
	//Printf("Event nTracks %d", nTracks);	
	AliTrackPointArray*  trackar=0x0;
	Bool_t l5;
	Bool_t l6;
	Int_t npoints;
	AliTrackPoint point;
	Int_t nlayer=0;
	Int_t id=0;
	Float_t chargeratio=0.0;
	Float_t charge=0.0;
	for(int itr=0;itr<nTracks;itr++)
	{
		AliESDtrack* track= esd->GetTrack(itr);

		if(TMath::Abs(track->Eta())>0.9)
			continue;
		if(track->GetP()>10.0)
			continue;	
		if (track->IsOn(AliESDtrack::kITSrefit)&&track->IsOn(AliESDtrack::kTPCrefit))
		{
			 l5=track->HasPointOnITSLayer(4);
			 l6=track->HasPointOnITSLayer(5);
			if (!(l5||l6))//only tracks with SSD point
				continue;
			Double_t tmpQESD[4]={-1.0,-1.0,-1.0,-1.0};
			track->GetITSdEdxSamples(tmpQESD);
			trackar=(AliTrackPointArray*)track->GetTrackPointArray(); 
			if(!trackar)
				continue;
			npoints=trackar->GetNPoints();			
			for(int itnp=0;itnp<npoints;itnp++)		
			{									
				if(trackar->GetPoint(point,itnp))
				{
					nlayer=AliGeomManager::VolUIDToLayer(point.GetVolumeID(), id);//layer number					
				}	
				else
					continue;		
				if(nlayer==5||nlayer==6)
				{
					TGeoHMatrix* geomatrix=AliGeomManager::GetMatrix(point.GetVolumeID());
					//geomatrix->Print();
					if(point.GetCharge()>0.0&&point.IsExtra()==kFALSE)//do not use additional clusters
					{		
						Double_t local[3]={0.0,0.0,0.0};
						Double_t global[3]={0.0,0.0,0.0};
						global[0]=point.GetX();
						global[1]=point.GetY();
						global[2]=point.GetZ();
						geomatrix->MasterToLocal(global,local);
						 chargeratio=point.GetChargeRatio();	
						 charge=point.GetCharge();
						 Float_t fQNnotcorr=charge*(1.0+chargeratio);
						 Float_t fQPnotcorr=charge*(1.0-chargeratio);						 
						 fHist6->Fill(fQNnotcorr,fQPnotcorr);
						// cout<<point.GetCharge()<<" "<<point.GetChargeRatio()<<" "<<nlayer<<" "<<id<<endl;
						if(nlayer==5&&tmpQESD[2]>0.0)
						{
							fHist1->Fill(id,chargeratio);
							if(track->GetP()>fPcut)
							{
								Float_t fQP=tmpQESD[2]*(1.0-chargeratio);
								Float_t fQN=tmpQESD[2]*(1.0+chargeratio);
								Int_t fPchip=Pstrip5(10000.0*local[0],10000.0*local[2])/128;
								Int_t fNchip=Nstrip5(10000.0*local[0],10000.0*local[2])/128;
								
								fHist2->Fill(id,tmpQESD[2]);
								fHist4->Fill(id*12+fPchip,fQP);	
								fHist4->Fill(id*12+6+fNchip,fQN);
								if(fdothecorrection)
									fHist5->Fill(id,0.5*(fQP*fcorrections[12*id+fPchip]+fQN*fcorrections[12*id+6+fNchip]));	
							}	
							fHist3->Fill(id,charge,chargeratio);
						}	
						if(nlayer==6&&tmpQESD[3]>0.0)
						{
							fHist1->Fill(id+748,chargeratio);
							if(track->GetP()>fPcut)
							{
							
								Float_t fQP=tmpQESD[3]*(1.0-chargeratio);
								Float_t fQN=tmpQESD[3]*(1.0+chargeratio);
								Int_t fPchip=Pstrip6(10000.0*local[0],10000.0*local[2])/128;
								Int_t fNchip=Nstrip6(10000.0*local[0],10000.0*local[2])/128;
							
								fHist2->Fill(id+748,tmpQESD[3]);
								fHist4->Fill((id+748)*12+fPchip,fQP);	
								fHist4->Fill((id+748)*12+6+fNchip,fQN);
								if(fdothecorrection)
									fHist5->Fill((id+748),0.5*(fQP*fcorrections[12*(id+748)+fPchip]+fQN*fcorrections[12*(id+748)+6+fNchip]));			
								
							}	
							fHist3->Fill(id+748,charge,chargeratio);	
						}
							
					}	
				}
			}	
		}
	}
	//Printf("Event nTracks end %d", nTracks);
	// Post output data.
	PostData(1,  fListOfHistos);
}      

//________________________________________________________________________
void AliAnalysisTaskdEdxSSDQA::Terminate(Option_t *) 
{
	//terminate
	Printf("end of Terminate");
}
//_____________________________________________________________________________
Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip5(Float_t x,Float_t z) const
{
	// P strip for layer 5 from local cooridnates 
	Float_t value;
	x=x+36432.5-223.9;
	z=z-(76000.0/7.0);
	value=-(3.0/38000.0)*z+x/95.0;
	return TMath::Nint(value);
}
//___________________________________________________________________________
Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip5(Float_t x,Float_t z)  const
{
	// N strip for layer 5 from local cooridnates 
	Float_t value;
	x=x+36432.5-223.9;
	z=z-(76000.0/7.0);
	value=x/95.0+(11.0/38000.0)*z;
	return TMath::Nint(value);
}
//___________________________________________________________________________
Int_t  AliAnalysisTaskdEdxSSDQA::Pstrip6(Float_t x,Float_t z)  const
{
	// P strip for layer 6 from local cooridnates 
	Float_t value;
	x=x-36432.5-223.9;
	z=z+(76000.0/7.0);
	value=-(3.0/38000.0)*z-x/95.0;
	return TMath::Nint(value);
}
//__________________________________________________________________________________
Int_t  AliAnalysisTaskdEdxSSDQA::Nstrip6(Float_t x,Float_t z)  const
{
	// N strip for layer 6 from local cooridnates 
	Float_t value;
	x=x-36432.5-223.9;
	z=z+(76000.0/7.0);
	value=x*(-1.0/95.0)+(11.0/38000.0)*z;
	return TMath::Nint(value);
}
//___________________________________________________________________________________
void  AliAnalysisTaskdEdxSSDQA::SetDoChipCorretions(const char* filename)
{
//Upload corrections for each chip only for test 
	cout<<filename<<endl;
	ifstream infile(filename);
	if(!infile.good())
		return;
	Float_t value=0.0;
	fdothecorrection=1;
	for (int i=0;i<1698;i++)
	{
		for (int j=0;j<6;j++)
		{
			infile>>value;
			cout<<value<<" ";
			if(value>0.0)
				fcorrections[i*12+j]=84.0/value;
			else
				fcorrections[i*12+j]=1.0;	
			infile>>value;
			infile>>value;
			cout<<value<<" ";
			if(value>0.0)
				fcorrections[i*12+j+6]=84.0/value;
			else
				fcorrections[i*12+j+6]=1.0;
			infile>>value;
		}
		cout<<endl;		
	}	
}
Commit	Line	Data
17ae74df	1	//SSD dEdx gain calibration task
	2	//to run in QA train
	3	//autor Marek Chojnacki
	4	//Marek.Chojnacki@cern.ch
	5	// last line
1625e300	6	#include "TChain.h"
	7	#include "TTree.h"
	8	#include "TH1F.h"
	9	#include "TH2F.h"
17ae74df	10	#include "TH3F.h"
1625e300	11	#include "TFile.h"
	12
	13	#include "AliAnalysisTaskdEdxSSDQA.h"
	14	#include "AliAnalysisManager.h"
	15
	16	#include "AliESDEvent.h"
	17	#include "AliESDInputHandler.h"
	18
	19
	20	#include "AliESDtrack.h"
	21	#include "AliExternalTrackParam.h"
	22	#include "Riostream.h"
	23
	24	#include "AliESDfriend.h"
	25	#include "AliESDfriendTrack.h"
	26	#include "AliTrackPointArray.h"
	27	#include "AliCDBManager.h"
	28	#include "AliGeomManager.h"
	29	#include "TMath.h"
17ae74df	30	#include "TGeoMatrix.h"
1625e300	31
	32	using namespace std;
	33	ClassImp(AliAnalysisTaskdEdxSSDQA)
	34
	35	//________________________________________________________________________
	36	AliAnalysisTaskdEdxSSDQA::AliAnalysisTaskdEdxSSDQA(const char *name)
17ae74df	37	: AliAnalysisTaskSE(name), fHist1(0),fHist2(0),fHist3(0),fHist4(0),fHist5(0),fHist6(0),fListOfHistos(0),fPcut(0.0),fdothecorrection(0)
1625e300	38	{
	39	// Constructor
	40	// Define input and output slots here
	41	// Input slot #1 works with a TChain
17ae74df	42	// fcorrections=new Float_t[1698*12];
	43	for(int i=0;i<1698*12;i++)
	44	fcorrections[i]=1.0;
1625e300	45	DefineOutput(1, TList::Class());
	46	}
	47
	48	//________________________________________________________________________
767feb9b	49	void AliAnalysisTaskdEdxSSDQA::UserCreateOutputObjects()
1625e300	50	{
	51	// Create histograms
	52	// Called once
	53
	54	fListOfHistos = new TList();
7335fa13	55	fListOfHistos->SetOwner();
1625e300	56	fHist1 =new TH2F("QAChargeRatio","QAChargeRatio;Module;CR",1698,-0.5,1697.5,80,-1.0,1.0);
	57	fListOfHistos->Add(fHist1);
	58	fHist2=new TH2F("QACharge","QACharge;Module;Q",1698,-0.5,1697.5,150,0,300);
	59	fListOfHistos->Add(fHist2);
17ae74df	60	fHist3=new TH3F("ChargeRatiovCharge","CRvQ;module;Q;CR",1698,-0.5,1697.5,80,0,1600,50,-1.0,1.0);
15b086e2	61	fListOfHistos->Add(fHist3);
17ae74df	62	fHist4=new TH2F("QAChargedinchips","chargedinchips;chip;Q",169812,-0.5,169812-0.5,150,0,300);
	63	fListOfHistos->Add(fHist4);
	64	if(fdothecorrection)
	65	{
	66	fHist5=new TH2F("QAChargeCorrected","chargedCorrected;Module;QCorrected",1698,-0.5,1697.5,150,0,300);
	67	fListOfHistos->Add(fHist5);
	68	}
	69	fHist6=new TH2F("QNvQP",";QN;QP",160,0,1600,1600,0,1600);
	70	fListOfHistos->Add(fHist6);
	71
	72	PostData(1, fListOfHistos);
1625e300	73	}
	74	//________________________________________________________________________
	75	void AliAnalysisTaskdEdxSSDQA::LocalInit()
	76	{
	77
	78	Printf("end of LocalInit");
	79	}
17ae74df	80	//______________________________________________________________________
	81	AliAnalysisTaskdEdxSSDQA::~AliAnalysisTaskdEdxSSDQA()
	82	{
	83	}
1625e300	84
	85	//________________________________________________________________________
	86
767feb9b	87	void AliAnalysisTaskdEdxSSDQA::UserExec(Option_t *)
1625e300	88	{
	89	// Main loop
	90	// Called for each event
	91	AliESDEvent* esd = dynamic_cast<AliESDEvent*> (InputEvent());
	92	if(!esd)
	93	{
	94	Printf("ERROR: Input ESD Event not available");
17ae74df	95	PostData(1, fListOfHistos);
1625e300	96	return;
	97	}
	98
	99	Int_t nTracks=esd->GetNumberOfTracks();
	100
	101
	102
	103	if(!ESDfriend())
	104	{
	105	Printf("problem with friend");
17ae74df	106	PostData(1, fListOfHistos);
1625e300	107	return;
1625e300	108	}
17ae74df	109	//Printf("Event nTracks %d", nTracks);
1625e300	110	AliTrackPointArray* trackar=0x0;
	111	Bool_t l5;
	112	Bool_t l6;
	113	Int_t npoints;
	114	AliTrackPoint point;
	115	Int_t nlayer=0;
	116	Int_t id=0;
15b086e2	117	Float_t chargeratio=0.0;
15b086e2	118	Float_t charge=0.0;
1625e300	119	for(int itr=0;itr<nTracks;itr++)
	120	{
	121	AliESDtrack* track= esd->GetTrack(itr);
	122
	123	if(TMath::Abs(track->Eta())>0.9)
	124	continue;
	125	if(track->GetP()>10.0)
	126	continue;
	127	if (track->IsOn(AliESDtrack::kITSrefit)&&track->IsOn(AliESDtrack::kTPCrefit))
	128	{
	129	l5=track->HasPointOnITSLayer(4);
	130	l6=track->HasPointOnITSLayer(5);
	131	if (!(l5\|\|l6))//only tracks with SSD point
	132	continue;
	133	Double_t tmpQESD[4]={-1.0,-1.0,-1.0,-1.0};
	134	track->GetITSdEdxSamples(tmpQESD);
	135	trackar=(AliTrackPointArray*)track->GetTrackPointArray();
	136	if(!trackar)
	137	continue;
	138	npoints=trackar->GetNPoints();
	139	for(int itnp=0;itnp<npoints;itnp++)
	140	{
	141	if(trackar->GetPoint(point,itnp))
17ae74df	142	{
	143	nlayer=AliGeomManager::VolUIDToLayer(point.GetVolumeID(), id);//layer number
	144	}
1625e300	145	else
	146	continue;
	147	if(nlayer==5\|\|nlayer==6)
	148	{
17ae74df	149	TGeoHMatrix* geomatrix=AliGeomManager::GetMatrix(point.GetVolumeID());
17ae74df	150	//geomatrix->Print();
1625e300	151	if(point.GetCharge()>0.0&&point.IsExtra()==kFALSE)//do not use additional clusters
1625e300	152	{
17ae74df	153	Double_t local[3]={0.0,0.0,0.0};
	154	Double_t global[3]={0.0,0.0,0.0};
	155	global[0]=point.GetX();
	156	global[1]=point.GetY();
	157	global[2]=point.GetZ();
	158	geomatrix->MasterToLocal(global,local);
15b086e2	159	chargeratio=point.GetChargeRatio();
17ae74df	160	charge=point.GetCharge();
	161	Float_t fQNnotcorr=charge*(1.0+chargeratio);
	162	Float_t fQPnotcorr=charge*(1.0-chargeratio);
	163	fHist6->Fill(fQNnotcorr,fQPnotcorr);
	164	// cout<<point.GetCharge()<<" "<<point.GetChargeRatio()<<" "<<nlayer<<" "<<id<<endl;
1625e300	165	if(nlayer==5&&tmpQESD[2]>0.0)
	166	{
	167	fHist1->Fill(id,chargeratio);
	168	if(track->GetP()>fPcut)
17ae74df	169	{
	170	Float_t fQP=tmpQESD[2]*(1.0-chargeratio);
	171	Float_t fQN=tmpQESD[2]*(1.0+chargeratio);
	172	Int_t fPchip=Pstrip5(10000.0local[0],10000.0local[2])/128;
	173	Int_t fNchip=Nstrip5(10000.0local[0],10000.0local[2])/128;
	174
1625e300	175	fHist2->Fill(id,tmpQESD[2]);
17ae74df	176	fHist4->Fill(id*12+fPchip,fQP);
	177	fHist4->Fill(id*12+6+fNchip,fQN);
	178	if(fdothecorrection)
	179	fHist5->Fill(id,0.5(fQPfcorrections[12id+fPchip]+fQNfcorrections[12*id+6+fNchip]));
	180	}
	181	fHist3->Fill(id,charge,chargeratio);
1625e300	182	}
	183	if(nlayer==6&&tmpQESD[3]>0.0)
	184	{
	185	fHist1->Fill(id+748,chargeratio);
	186	if(track->GetP()>fPcut)
17ae74df	187	{
	188
	189	Float_t fQP=tmpQESD[3]*(1.0-chargeratio);
	190	Float_t fQN=tmpQESD[3]*(1.0+chargeratio);
	191	Int_t fPchip=Pstrip6(10000.0local[0],10000.0local[2])/128;
	192	Int_t fNchip=Nstrip6(10000.0local[0],10000.0local[2])/128;
	193
1625e300	194	fHist2->Fill(id+748,tmpQESD[3]);
17ae74df	195	fHist4->Fill((id+748)*12+fPchip,fQP);
	196	fHist4->Fill((id+748)*12+6+fNchip,fQN);
	197	if(fdothecorrection)
	198	fHist5->Fill((id+748),0.5(fQPfcorrections[12(id+748)+fPchip]+fQNfcorrections[12*(id+748)+6+fNchip]));
	199
	200	}
	201	fHist3->Fill(id+748,charge,chargeratio);
15b086e2	202	}
17ae74df	203
1625e300	204	}
	205	}
	206	}
	207	}
	208	}
17ae74df	209	//Printf("Event nTracks end %d", nTracks);
1625e300	210	// Post output data.
	211	PostData(1, fListOfHistos);
	212	}
	213
	214	//________________________________________________________________________
	215	void AliAnalysisTaskdEdxSSDQA::Terminate(Option_t *)
	216	{
17ae74df	217	//terminate
1625e300	218	Printf("end of Terminate");
	219	}
	220	//_____________________________________________________________________________
17ae74df	221	Int_t AliAnalysisTaskdEdxSSDQA::Pstrip5(Float_t x,Float_t z) const
	222	{
	223	// P strip for layer 5 from local cooridnates
	224	Float_t value;
	225	x=x+36432.5-223.9;
	226	z=z-(76000.0/7.0);
	227	value=-(3.0/38000.0)*z+x/95.0;
	228	return TMath::Nint(value);
	229	}
	230	//___________________________________________________________________________
	231	Int_t AliAnalysisTaskdEdxSSDQA::Nstrip5(Float_t x,Float_t z) const
	232	{
	233	// N strip for layer 5 from local cooridnates
	234	Float_t value;
	235	x=x+36432.5-223.9;
	236	z=z-(76000.0/7.0);
	237	value=x/95.0+(11.0/38000.0)*z;
	238	return TMath::Nint(value);
	239	}
	240	//___________________________________________________________________________
	241	Int_t AliAnalysisTaskdEdxSSDQA::Pstrip6(Float_t x,Float_t z) const
	242	{
	243	// P strip for layer 6 from local cooridnates
	244	Float_t value;
	245	x=x-36432.5-223.9;
	246	z=z+(76000.0/7.0);
	247	value=-(3.0/38000.0)*z-x/95.0;
	248	return TMath::Nint(value);
	249	}
	250	//__________________________________________________________________________________
	251	Int_t AliAnalysisTaskdEdxSSDQA::Nstrip6(Float_t x,Float_t z) const
	252	{
	253	// N strip for layer 6 from local cooridnates
	254	Float_t value;
	255	x=x-36432.5-223.9;
	256	z=z+(76000.0/7.0);
	257	value=x(-1.0/95.0)+(11.0/38000.0)z;
	258	return TMath::Nint(value);
	259	}
	260	//___________________________________________________________________________________
	261	void AliAnalysisTaskdEdxSSDQA::SetDoChipCorretions(const char* filename)
	262	{
	263	//Upload corrections for each chip only for test
	264	cout<<filename<<endl;
	265	ifstream infile(filename);
	266	if(!infile.good())
	267	return;
	268	Float_t value=0.0;
	269	fdothecorrection=1;
	270	for (int i=0;i<1698;i++)
	271	{
	272	for (int j=0;j<6;j++)
	273	{
	274	infile>>value;
	275	cout<<value<<" ";
	276	if(value>0.0)
	277	fcorrections[i*12+j]=84.0/value;
	278	else
	279	fcorrections[i*12+j]=1.0;
	280	infile>>value;
	281	infile>>value;
	282	cout<<value<<" ";
	283	if(value>0.0)
	284	fcorrections[i*12+j+6]=84.0/value;
285	else
286	fcorrections[i*12+j+6]=1.0;
287	infile>>value;
288	}
289	cout<<endl;
290	}
291	}
1625e300	292