[u/mrichter/AliRoot.git] / TRD / qaAnalysis / AliTRDqaJPsi.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 iLs" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//
// This class is a part of a package of high level QA monitoring for TRD.
//
// S. Radomski
// radomski@physi.uni-heidelberg.de
// March 2008
//

#include "AliTRDqaJPsi.h"
#include "AliTRDqaAT.h"

#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TChain.h"

#include "AliESDEvent.h" 
#include "AliESDtrack.h"
#include "AliKFParticle.h"

#include "TLorentzVector.h"

//______________________________________________________________________________

AliTRDqaJPsi::AliTRDqaJPsi() 
  : AliAnalysisTask("",""),  
    fChain(0),
    fESD(0),
    fOutputContainer(0),
    fnKFtracks(0)
{
  //
  // default dummy constructor
  //

  for (Int_t i = 0; i < 1000; i++) {
    fSM[i] = 0;
  }
 
}

//______________________________________________________________________________
AliTRDqaJPsi:: AliTRDqaJPsi(const AliTRDqaJPsi & /*trd*/)
  : AliAnalysisTask("",""),  
    fChain(0),
    fESD(0),
    fOutputContainer(0),
    fnKFtracks(0)
{
  //
  // Copy constructor
  //

  for (Int_t i = 0; i < 1000; i++) {
    fSM[i] = 0;
  }

}

//______________________________________________________________________________
AliTRDqaJPsi::AliTRDqaJPsi(const char *name) 
  : AliAnalysisTask(name,""),  
    fChain(0),
    fESD(0),
    fOutputContainer(0),
    fnKFtracks(0)
    
{
  // Constructor.
  // Input slot #0 works with an Ntuple
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TH1 container
  DefineOutput(0,  TObjArray::Class()) ; 
  for (Int_t i = 0; i < 1000; i++) {
    fSM[i] = 0;
  }

}

//______________________________________________________________________________
void AliTRDqaJPsi::ConnectInputData(const Option_t *)
{
  // Initialisation of branch container and histograms 

  //AliInfo(Form("*** Initialization of %s", GetName())) ; 

  fChain = (TChain*)GetInputData(0);
  fESD = new AliESDEvent();
  fESD->ReadFromTree(fChain);
}

//________________________________________________________________________
void AliTRDqaJPsi::CreateOutputObjects()
{
  //
  // Create the objects that contain the analysis output
  //
  
  Int_t c = 0;
  fOutputContainer = new TObjArray(100);
  
  const char *charge[2] = {"Neg", "Pos"};
  
  // build histograms
  
  for(Int_t i=0; i<fgknSteps; i++) {

    fStatus[i] = new TH1D(Form("status_%d", i), "status", 32, -0.5, 31.5);
    fOutputContainer->AddAt(fStatus[i], c++);
    
    fInvMass[i] = new TH1D(Form("mass_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
    fOutputContainer->AddAt(fInvMass[i], c++);
    
    fInvMassVec[i] = new TH1D(Form("massVec_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
    fOutputContainer->AddAt(fInvMassVec[i], c++);
    
    fInvMassDiff[i] = new TH1D(Form("massDiff_%d", i), ";m_{inv} (GeV);", 100, -1, 1);
    fOutputContainer->AddAt(fInvMassDiff[i], c++);
    
    fAngleSM[i] = new TH1D(Form("angleSM_%d", i), ";#delta SM", 19, -0.5, 18.5); 	
    fOutputContainer->AddAt(fAngleSM[i], c++);
    
    fPtAngle[i] = new TH2D(Form("ptAngle_%d", i), ";p_{T} (GeV/c);#delta SM",
			   20, 0, 5, 10, -0.5, 9.5);
    fOutputContainer->AddAt(fPtAngle[i], c++);
    

    for(Int_t j=0; j<2; j++) {
      fnTracks[j*fgknSteps+i] = 
	new TH1D(Form("nTracks%s_%d", charge[j],i), Form("%s;number of tracks",charge[j]), 100, -0.5, 99.5);
      fPt[j*fgknSteps+i] = new TH1D(Form("pt%s_%d", charge[j], i), Form("%s;p_{T} (GeV/c)", charge[j]), 100, 0, 5);
      fPID[j*fgknSteps+i] = new TH1D(Form("pid%s_%d", charge[j], i), ";electron LQ", 100, 0, 1);

      fOutputContainer->AddAt(fnTracks[j*fgknSteps+i], c++);
      fOutputContainer->AddAt(fPt[j*fgknSteps+i], c++); 
      fOutputContainer->AddAt(fPID[j*fgknSteps+i], c++); 
    }
  }

  //TH2D *fnGoodTracks;  

  printf("n hist = %d\n", c);
}
//______________________________________________________________________________
void AliTRDqaJPsi::Exec(Option_t *) 
{
  /*
    Selection steps:
    - Parameters In and Out
    - TRDrefit bit
    - TRDpid and quality
    - pt > 0.8 GeV/c
    - PID > 0.9
   */


  // Process one event
  Long64_t entry = fChain->GetReadEntry() ;
  if (!(entry%100)) Info("Exec", "Entry = %lld", entry);

  // Processing of one event 
   
  if (!fESD) {
    //AliError("fESD is not connected to the input!") ; 
    return ; 
  }
  
  
  Int_t nTracks = fESD->GetNumberOfTracks();
  Int_t cTracks[2*fgknSteps] = {0,0,0,0,0,0,0,0,0,0};
  fnKFtracks = 0;

  // track loop
  for(Int_t i=0; i<nTracks; i++) {
    
    //
    // track selection 
    //
    // param in and Out
    // TRDrefit and TRDPid bit
    //
 
    AliESDtrack *track = fESD->GetTrack(i);
    const AliExternalTrackParam *paramOut = track->GetOuterParam();
    const AliExternalTrackParam *paramIn = track->GetInnerParam();

    // long track ..
    if (!paramIn) continue;
    if (!paramOut) continue;

    Int_t step    = 0;
    Int_t charge  = (track->Charge() > 0) ? 1 : 0;
    UInt_t status = track->GetStatus();
    Double_t pt   = track->Pt();
    Double_t pid  = track->GetTRDpid(AliPID::kElectron);
    
    Double_t esdPid[5];
    track->GetESDpid(esdPid); 

    // create a kalman particle
    Int_t pdg = (charge == 0)? -11 : 11;
    for(Int_t k=0; k<fgknSteps; k++) fInSample[fnKFtracks][k] = 0;  
 
    fVec[fnKFtracks] = CreateVector(track);
    fTracks[fnKFtracks] = new AliKFParticle(*track, pdg);
    fSM[fnKFtracks] = AliTRDqaAT::GetSector(paramOut->GetAlpha());
    fnKFtracks++;

    //AliTRDqaAT::PrintPID(track);

    // apply the cuts

    cTracks[fgknSteps *charge + step]++;
    FillHist(track, step++);

    if (!(status & AliESDtrack::kTRDrefit)) continue;

    cTracks[fgknSteps *charge + step]++;    
    FillHist(track, step++);    

    if (!(status & AliESDtrack::kTRDpid)) continue;
    if (track->GetTRDntracklets() < 6) continue;

    cTracks[fgknSteps *charge + step]++;
    FillHist(track, step++);  

    if (pt < 0.8) continue;

    cTracks[fgknSteps *charge + step]++;
    FillHist(track, step++);      

    if (pid < 0.3) continue; // 
    //if (esdPid[AliPID::kElectron] < 0.5) continue;
    
    cTracks[fgknSteps *charge + step]++;
    FillHist(track, step);  

    for(Int_t k=0; k<2*fgknSteps; k++) fnTracks[k]->Fill(cTracks[k]);
  }

  // calculate invariant mass

  for(Int_t k=0; k<fgknSteps; k++) {
    for(Int_t i=0; i<fnKFtracks; i++) {
      if (!fInSample[i][k]) continue;
      for(Int_t j=i+1; j<fnKFtracks; j++) {
	if (!fInSample[j][k]) continue;
	AliKFParticle jpsi(*(fTracks[i]), *(fTracks[j]));
	TLorentzVector jpsiVec = (*(fVec[i])) + (*fVec[j]);
	fInvMass[k]->Fill(jpsi.GetMass());
	fInvMassVec[k]->Fill(jpsiVec.M());
	fInvMassDiff[k]->Fill(jpsiVec.M() - jpsi.GetMass());
	
	if (jpsi.GetMass() > 2.5 && jpsi.GetMass() < 3.5) {
	  Int_t dSM = TMath::Abs(fSM[i] - fSM[j]);
	  if (dSM > 9) dSM = (18-dSM);
	  fAngleSM[k]->Fill(dSM);
	  fPtAngle[k]->Fill(TMath::Hypot(jpsi.GetPx(), jpsi.GetPy()), dSM);
	}
	
      }
    }
  }

  PostData(0, fOutputContainer);
}

//______________________________________________________________________________
void AliTRDqaJPsi::Terminate(Option_t *)
{
  // save histograms
  fOutputContainer = (TObjArray*)GetOutputData(0);
  
  TFile *file = new TFile("outJPsi.root", "RECREATE");
  fOutputContainer->Write();
  
  file->Flush();
  file->Close();
  delete file;  

  //for(Int_t i=0; i<fOutputContainer->GetEntries(); i++) {
  //  TObject *obj = fOu
  // }
}

//______________________________________________________________________________

void AliTRDqaJPsi::FillHist(AliESDtrack *track, Int_t step) {
  //
  // Fill the histograms
  //

  Int_t charge  = (track->Charge() > 0) ? 1 : 0;
  UInt_t status = track->GetStatus();
  Double_t pt   = track->Pt();
  Double_t pid  = track->GetTRDpid(AliPID::kElectron);
  
  Double_t esdPid[5];
  track->GetESDpid(esdPid);

  Int_t id = charge * fgknSteps + step;
  AliTRDqaAT::FillStatus(fStatus[step], status);
  fPt[id]->Fill(pt);
  fPID[id]->Fill(pid);
  //fPID[id]->Fill(esdPid[AliPID::kElectron]);

  fInSample[fnKFtracks-1][step] = 1;
}

//______________________________________________________________________________

TLorentzVector *AliTRDqaJPsi::CreateVector(AliESDtrack *track) {

  TLorentzVector *vec = new TLorentzVector();
  vec->SetXYZM(track->Px(), track->Py(), track->Pz(), 0);
  return vec;
}

//______________________________________________________________________________
Commit	Line	Data
8e2f611a	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 iLs" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	//
	19	// This class is a part of a package of high level QA monitoring for TRD.
	20	//
	21	// S. Radomski
	22	// radomski@physi.uni-heidelberg.de
	23	// March 2008
	24	//
	25
	26	#include "AliTRDqaJPsi.h"
	27	#include "AliTRDqaAT.h"
	28
	29	#include "TFile.h"
	30	#include "TH1D.h"
	31	#include "TH2D.h"
	32	#include "TChain.h"
	33
	34	#include "AliESDEvent.h"
	35	#include "AliESDtrack.h"
	36	#include "AliKFParticle.h"
	37
	38	#include "TLorentzVector.h"
	39
8e2f611a	40	//______________________________________________________________________________
	41
	42	AliTRDqaJPsi::AliTRDqaJPsi()
	43	: AliAnalysisTask("",""),
	44	fChain(0),
	45	fESD(0),
	46	fOutputContainer(0),
	47	fnKFtracks(0)
	48	{
	49	//
	50	// default dummy constructor
	51	//
02f3bfcc	52
	53	for (Int_t i = 0; i < 1000; i++) {
	54	fSM[i] = 0;
	55	}
8e2f611a	56
8e2f611a	57	}
8e2f611a	58
c336eab0	59	//______________________________________________________________________________
c336eab0	60	AliTRDqaJPsi:: AliTRDqaJPsi(const AliTRDqaJPsi & /trd/)
8e2f611a	61	: AliAnalysisTask("",""),
	62	fChain(0),
	63	fESD(0),
	64	fOutputContainer(0),
	65	fnKFtracks(0)
	66	{
	67	//
c336eab0	68	// Copy constructor
8e2f611a	69	//
8e2f611a	70
a987273c	71	for (Int_t i = 0; i < 1000; i++) {
	72	fSM[i] = 0;
	73	}
	74
8e2f611a	75	}
8e2f611a	76
8e2f611a	77	//______________________________________________________________________________
	78	AliTRDqaJPsi::AliTRDqaJPsi(const char *name)
	79	: AliAnalysisTask(name,""),
	80	fChain(0),
	81	fESD(0),
	82	fOutputContainer(0),
	83	fnKFtracks(0)
	84
	85	{
	86	// Constructor.
	87	// Input slot #0 works with an Ntuple
	88	DefineInput(0, TChain::Class());
	89	// Output slot #0 writes into a TH1 container
	90	DefineOutput(0, TObjArray::Class()) ;
02f3bfcc	91	for (Int_t i = 0; i < 1000; i++) {
	92	fSM[i] = 0;
	93	}
	94
8e2f611a	95	}
	96
	97	//______________________________________________________________________________
	98	void AliTRDqaJPsi::ConnectInputData(const Option_t *)
	99	{
	100	// Initialisation of branch container and histograms
	101
	102	//AliInfo(Form("*** Initialization of %s", GetName())) ;
	103
	104	fChain = (TChain*)GetInputData(0);
	105	fESD = new AliESDEvent();
	106	fESD->ReadFromTree(fChain);
	107	}
	108
	109	//________________________________________________________________________
	110	void AliTRDqaJPsi::CreateOutputObjects()
	111	{
c336eab0	112	//
	113	// Create the objects that contain the analysis output
	114	//
8e2f611a	115
	116	Int_t c = 0;
	117	fOutputContainer = new TObjArray(100);
	118
	119	const char *charge[2] = {"Neg", "Pos"};
	120
	121	// build histograms
	122
c336eab0	123	for(Int_t i=0; i<fgknSteps; i++) {
8e2f611a	124
	125	fStatus[i] = new TH1D(Form("status_%d", i), "status", 32, -0.5, 31.5);
	126	fOutputContainer->AddAt(fStatus[i], c++);
	127
	128	fInvMass[i] = new TH1D(Form("mass_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
	129	fOutputContainer->AddAt(fInvMass[i], c++);
	130
	131	fInvMassVec[i] = new TH1D(Form("massVec_%d", i), ";m_{inv} (GeV);", 100, 0, 5);
	132	fOutputContainer->AddAt(fInvMassVec[i], c++);
	133
	134	fInvMassDiff[i] = new TH1D(Form("massDiff_%d", i), ";m_{inv} (GeV);", 100, -1, 1);
	135	fOutputContainer->AddAt(fInvMassDiff[i], c++);
	136
	137	fAngleSM[i] = new TH1D(Form("angleSM_%d", i), ";#delta SM", 19, -0.5, 18.5);
	138	fOutputContainer->AddAt(fAngleSM[i], c++);
	139
	140	fPtAngle[i] = new TH2D(Form("ptAngle_%d", i), ";p_{T} (GeV/c);#delta SM",
	141	20, 0, 5, 10, -0.5, 9.5);
	142	fOutputContainer->AddAt(fPtAngle[i], c++);
	143
	144
	145	for(Int_t j=0; j<2; j++) {
c336eab0	146	fnTracks[j*fgknSteps+i] =
8e2f611a	147	new TH1D(Form("nTracks%s_%d", charge[j],i), Form("%s;number of tracks",charge[j]), 100, -0.5, 99.5);
c336eab0	148	fPt[j*fgknSteps+i] = new TH1D(Form("pt%s_%d", charge[j], i), Form("%s;p_{T} (GeV/c)", charge[j]), 100, 0, 5);
c336eab0	149	fPID[j*fgknSteps+i] = new TH1D(Form("pid%s_%d", charge[j], i), ";electron LQ", 100, 0, 1);
8e2f611a	150
c336eab0	151	fOutputContainer->AddAt(fnTracks[j*fgknSteps+i], c++);
	152	fOutputContainer->AddAt(fPt[j*fgknSteps+i], c++);
	153	fOutputContainer->AddAt(fPID[j*fgknSteps+i], c++);
8e2f611a	154	}
	155	}
	156
	157	//TH2D *fnGoodTracks;
	158
	159	printf("n hist = %d\n", c);
	160	}
	161	//______________________________________________________________________________
	162	void AliTRDqaJPsi::Exec(Option_t *)
	163	{
	164	/*
	165	Selection steps:
	166	- Parameters In and Out
	167	- TRDrefit bit
	168	- TRDpid and quality
	169	- pt > 0.8 GeV/c
	170	- PID > 0.9
	171	*/
	172
	173
	174	// Process one event
	175	Long64_t entry = fChain->GetReadEntry() ;
20a6e9c2	176	if (!(entry%100)) Info("Exec", "Entry = %lld", entry);
8e2f611a	177
	178	// Processing of one event
	179
	180	if (!fESD) {
	181	//AliError("fESD is not connected to the input!") ;
	182	return ;
	183	}
	184
	185
	186	Int_t nTracks = fESD->GetNumberOfTracks();
c336eab0	187	Int_t cTracks[2*fgknSteps] = {0,0,0,0,0,0,0,0,0,0};
8e2f611a	188	fnKFtracks = 0;
	189
	190	// track loop
	191	for(Int_t i=0; i<nTracks; i++) {
	192
	193	//
	194	// track selection
	195	//
	196	// param in and Out
	197	// TRDrefit and TRDPid bit
	198	//
	199
	200	AliESDtrack *track = fESD->GetTrack(i);
	201	const AliExternalTrackParam *paramOut = track->GetOuterParam();
	202	const AliExternalTrackParam *paramIn = track->GetInnerParam();
	203
	204	// long track ..
	205	if (!paramIn) continue;
	206	if (!paramOut) continue;
	207
	208	Int_t step = 0;
	209	Int_t charge = (track->Charge() > 0) ? 1 : 0;
	210	UInt_t status = track->GetStatus();
	211	Double_t pt = track->Pt();
	212	Double_t pid = track->GetTRDpid(AliPID::kElectron);
	213
	214	Double_t esdPid[5];
	215	track->GetESDpid(esdPid);
	216
	217	// create a kalman particle
	218	Int_t pdg = (charge == 0)? -11 : 11;
c336eab0	219	for(Int_t k=0; k<fgknSteps; k++) fInSample[fnKFtracks][k] = 0;
8e2f611a	220
7bce990c	221	fVec[fnKFtracks] = CreateVector(track);
7bce990c	222	fTracks[fnKFtracks] = new AliKFParticle(*track, pdg);
8e2f611a	223	fSM[fnKFtracks] = AliTRDqaAT::GetSector(paramOut->GetAlpha());
	224	fnKFtracks++;
	225
	226	//AliTRDqaAT::PrintPID(track);
	227
	228	// apply the cuts
	229
c336eab0	230	cTracks[fgknSteps *charge + step]++;
8e2f611a	231	FillHist(track, step++);
	232
	233	if (!(status & AliESDtrack::kTRDrefit)) continue;
	234
c336eab0	235	cTracks[fgknSteps *charge + step]++;
8e2f611a	236	FillHist(track, step++);
	237
	238	if (!(status & AliESDtrack::kTRDpid)) continue;
ed15ef4f	239	if (track->GetTRDntracklets() < 6) continue;
8e2f611a	240
c336eab0	241	cTracks[fgknSteps *charge + step]++;
8e2f611a	242	FillHist(track, step++);
	243
	244	if (pt < 0.8) continue;
	245
c336eab0	246	cTracks[fgknSteps *charge + step]++;
8e2f611a	247	FillHist(track, step++);
	248
	249	if (pid < 0.3) continue; //
	250	//if (esdPid[AliPID::kElectron] < 0.5) continue;
	251
c336eab0	252	cTracks[fgknSteps *charge + step]++;
8e2f611a	253	FillHist(track, step);
8e2f611a	254
c336eab0	255	for(Int_t k=0; k<2*fgknSteps; k++) fnTracks[k]->Fill(cTracks[k]);
8e2f611a	256	}
	257
	258	// calculate invariant mass
	259
c336eab0	260	for(Int_t k=0; k<fgknSteps; k++) {
8e2f611a	261	for(Int_t i=0; i<fnKFtracks; i++) {
	262	if (!fInSample[i][k]) continue;
	263	for(Int_t j=i+1; j<fnKFtracks; j++) {
	264	if (!fInSample[j][k]) continue;
7bce990c	265	AliKFParticle jpsi((fTracks[i]), (fTracks[j]));
7bce990c	266	TLorentzVector jpsiVec = ((fVec[i])) + (fVec[j]);
8e2f611a	267	fInvMass[k]->Fill(jpsi.GetMass());
	268	fInvMassVec[k]->Fill(jpsiVec.M());
	269	fInvMassDiff[k]->Fill(jpsiVec.M() - jpsi.GetMass());
	270
	271	if (jpsi.GetMass() > 2.5 && jpsi.GetMass() < 3.5) {
	272	Int_t dSM = TMath::Abs(fSM[i] - fSM[j]);
	273	if (dSM > 9) dSM = (18-dSM);
	274	fAngleSM[k]->Fill(dSM);
	275	fPtAngle[k]->Fill(TMath::Hypot(jpsi.GetPx(), jpsi.GetPy()), dSM);
	276	}
	277
	278	}
	279	}
	280	}
	281
	282	PostData(0, fOutputContainer);
	283	}
	284
	285	//______________________________________________________________________________
	286	void AliTRDqaJPsi::Terminate(Option_t *)
	287	{
	288	// save histograms
	289	fOutputContainer = (TObjArray*)GetOutputData(0);
	290
	291	TFile *file = new TFile("outJPsi.root", "RECREATE");
	292	fOutputContainer->Write();
	293
	294	file->Flush();
	295	file->Close();
	296	delete file;
	297
	298	//for(Int_t i=0; i<fOutputContainer->GetEntries(); i++) {
	299	// TObject *obj = fOu
	300	// }
	301	}
	302
	303	//______________________________________________________________________________
	304
	305	void AliTRDqaJPsi::FillHist(AliESDtrack *track, Int_t step) {
c336eab0	306	//
	307	// Fill the histograms
	308	//
8e2f611a	309
	310	Int_t charge = (track->Charge() > 0) ? 1 : 0;
	311	UInt_t status = track->GetStatus();
	312	Double_t pt = track->Pt();
	313	Double_t pid = track->GetTRDpid(AliPID::kElectron);
	314
	315	Double_t esdPid[5];
	316	track->GetESDpid(esdPid);
	317
c336eab0	318	Int_t id = charge * fgknSteps + step;
8e2f611a	319	AliTRDqaAT::FillStatus(fStatus[step], status);
	320	fPt[id]->Fill(pt);
	321	fPID[id]->Fill(pid);
	322	//fPID[id]->Fill(esdPid[AliPID::kElectron]);
	323
	324	fInSample[fnKFtracks-1][step] = 1;
	325	}
	326
	327	//______________________________________________________________________________
	328
	329	TLorentzVector AliTRDqaJPsi::CreateVector(AliESDtrack track) {
	330
	331	TLorentzVector *vec = new TLorentzVector();
	332	vec->SetXYZM(track->Px(), track->Py(), track->Pz(), 0);
	333	return vec;
	334	}
	335
	336	//______________________________________________________________________________