[u/mrichter/AliRoot.git] / FMD / scripts / DrawHitsRecs.C

//____________________________________________________________________
//
// $Id$
//
// Script that contains a class to draw eloss from hits, versus ADC
// counts from digits, using the AliFMDInputHits class in the util library. 
//
// It draws the energy loss versus the p/(mq^2).  It can be overlayed
// with the Bethe-Bloc curve to show how the simulation behaves
// relative to the expected. 
//
// Use the script `Compile.C' to compile this class using ACLic. 
//
#include <TH2D.h>
#include <AliFMDHit.h>
#include <AliFMDDigit.h>
#include <AliFMDRecPoint.h>
#include <AliFMDInput.h>
#include <AliFMDEdepMap.h>
#include <AliFMDFloatMap.h>
#include <iostream>
#include <TStyle.h>
#include <TArrayF.h>
#include <TCanvas.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TTree.h>
#include <AliStack.h>
#include <AliLog.h>

//____________________________________________________________________
/** @class DrawHitsRecs
    @brief Draw hit energy loss versus rec point mult
    @code 
    Root> .L Compile.C
    Root> Compile("DrawHitsRecs.C")
    Root> DrawHitsRecs c
    Root> c.Run();
    @endcode
    @ingroup FMD_script
 */
class DrawHitsRecs : public AliFMDInput
{
private:
  TH2D* fHitEvsAdc;
  TH2D* fHitEvsRecM;  // Histogram 
  TH2D* fHitEvsRecE;  // Histogram 
  TH1D* fDiffE;       // Histogram 
  TH2D* fHitsVsRecM;  // Histogram 
  AliFMDEdepMap  fMap;
  AliFMDFloatMap fEta;
  AliFMDFloatMap fPhi;
  AliFMDFloatMap fMult;
  Bool_t fPrimary;
public:
  //__________________________________________________________________
  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
  {
    TArrayF bins(n+1);
    Float_t dp   = n / TMath::Log10(max / min);
    Float_t pmin = TMath::Log10(min);
    bins[0]      = min;
    for (Int_t i = 1; i < n+1; i++) {
      Float_t p = pmin + i / dp;
      bins[i]   = TMath::Power(10, p);
    }
    return bins;
  }
  //__________________________________________________________________
  DrawHitsRecs(Bool_t primary=kFALSE,
	       Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
	       Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5) 
  { 
    fPrimary = primary;
    AddLoad(kRecPoints);
    AddLoad(kHits);
    AddLoad(kDigits);
    if (fPrimary) AddLoad(kKinematics);
    TArrayF eloss(MakeLogScale(n, emin, emax));
    TArrayF mults(m+1);
    mults[0] = mmin;
    for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;

    fHitEvsAdc  = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
			   n, emin, emax, 1025, -.5, 1024.5);
    fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
    fHitEvsAdc->SetYTitle("ADC");
    
    fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}", 
			   eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
    fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
    fHitEvsRecM->SetYTitle("M_{rec}");

    fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}", 
			    n, emin, emax, n, emin, emax);
    fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
    fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");


    fDiffE = new TH1D("diffE", 
		      "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}", 
		      1100, -1, 1.1);
    fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
    
    Double_t omin = -.5;
    Double_t omax = 7.5;
    Int_t    o    = 8;
    fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
			   o, omin, omax, m, mmin, mmax);
    fHitsVsRecM->SetXTitle("# of Hits");
    fHitsVsRecM->SetYTitle("M_{rec}");
  }
  //__________________________________________________________________
  /** Begining of event
      @param ev Event number
      @return @c false on error */
  Bool_t Begin(Int_t ev) 
  {
    fMap.Reset();
    return AliFMDInput::Begin(ev);
  }
  //__________________________________________________________________
  Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
  {
    // Cache the energy loss 
    if (!hit) return kFALSE;
    UShort_t det = hit->Detector();
    Char_t   rng = hit->Ring();
    UShort_t sec = hit->Sector();
    UShort_t str = hit->Strip();
    if (str > 511) {
      AliWarning(Form("Bad strip number %d in hit", str));
      return kTRUE;
    }
    if (fPrimary) {
      TParticle* kine = fStack->Particle(hit->Track());
      if (!kine) return kTRUE;
      if (!kine->IsPrimary()) return kTRUE;
    }
    
    fMap(det, rng, sec, str).fEdep += hit->Edep();
    fMap(det, rng, sec, str).fN++;
    return kTRUE;
  }

  //__________________________________________________________________
  Bool_t ProcessDigit(AliFMDDigit* digit)
  {
    if (!digit) return kTRUE;

    UShort_t det = digit->Detector();
    Char_t   rng = digit->Ring();
    UShort_t sec = digit->Sector();
    UShort_t str = digit->Strip();
    if (str > 511) {
      AliWarning(Form("Bad strip number %d in digit", str));
      return kTRUE;
    }
    Double_t edep = fMap(det, rng, sec, str).fEdep;
    if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
      
    return kTRUE;
  }

  //__________________________________________________________________
  Bool_t ProcessRecPoint(AliFMDRecPoint* single) 
  {
    if (!single) return kTRUE;
    UShort_t det = single->Detector();
    Char_t   rng = single->Ring();
    UShort_t sec = single->Sector();
    UShort_t str = single->Strip();
    if (str > 511) {
      AliWarning(Form("Bad strip number %d in single", str));
      return kTRUE;
    }
    Double_t edep = fMap(det, rng, sec, str).fEdep;
    Int_t    nhit = fMap(det, rng, sec, str).fN;
    if (edep > 0) {
      fHitEvsRecM->Fill(edep, single->Particles());
      fHitEvsRecE->Fill(edep, single->Edep());
      fDiffE->Fill((single->Edep() - edep) / edep);
    }
    if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
    return kTRUE;
  }
  //__________________________________________________________________
  Bool_t Finish()
  {
    gStyle->SetPalette(1);
    gStyle->SetOptTitle(0);
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderSize(0);
    gStyle->SetPadColor(0);
    gStyle->SetPadBorderSize(0);

    new TCanvas("c0", fHitEvsAdc->GetTitle());
    fHitEvsAdc->SetStats(kFALSE);
    fHitEvsAdc->Draw("COLZ");

    new TCanvas("c1", fHitEvsRecM->GetTitle());
    fHitEvsRecM->SetStats(kFALSE);
    fHitEvsRecM->Draw("COLZ");

    new TCanvas("c2", fHitEvsRecE->GetTitle());
    fHitEvsRecE->SetStats(kFALSE);
    fHitEvsRecE->Draw("COLZ");

    new TCanvas("c3", fDiffE->GetTitle());
    fDiffE->Draw();

    new TCanvas("c4", fHitsVsRecM->GetTitle());
    fHitsVsRecM->SetStats(kFALSE);
    fHitsVsRecM->Draw("COLZ");

    return kTRUE;
  }

  ClassDef(DrawHitsRecs,0);
};

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
bf000c32	1	//____________________________________________________________________
8f6ee336	2	//
	3	// $Id$
	4	//
	5	// Script that contains a class to draw eloss from hits, versus ADC
	6	// counts from digits, using the AliFMDInputHits class in the util library.
	7	//
	8	// It draws the energy loss versus the p/(mq^2). It can be overlayed
	9	// with the Bethe-Bloc curve to show how the simulation behaves
	10	// relative to the expected.
	11	//
	12	// Use the script `Compile.C' to compile this class using ACLic.
	13	//
	14	#include <TH2D.h>
	15	#include <AliFMDHit.h>
8f6ee336	16	#include <AliFMDDigit.h>
15b17c89	17	#include <AliFMDRecPoint.h>
8f6ee336	18	#include <AliFMDInput.h>
	19	#include <AliFMDEdepMap.h>
	20	#include <AliFMDFloatMap.h>
	21	#include <iostream>
	22	#include <TStyle.h>
	23	#include <TArrayF.h>
	24	#include <TCanvas.h>
	25	#include <TParticle.h>
	26	#include <TClonesArray.h>
	27	#include <TTree.h>
	28	#include <AliStack.h>
	29	#include <AliLog.h>
	30
bf000c32	31	//____________________________________________________________________
9b48326f	32	/** @class DrawHitsRecs
	33	@brief Draw hit energy loss versus rec point mult
	34	@code
	35	Root> .L Compile.C
	36	Root> Compile("DrawHitsRecs.C")
	37	Root> DrawHitsRecs c
	38	Root> c.Run();
	39	@endcode
	40	@ingroup FMD_script
	41	*/
15b17c89	42	class DrawHitsRecs : public AliFMDInput
8f6ee336	43	{
8f6ee336	44	private:
15b17c89	45	TH2D* fHitEvsAdc;
	46	TH2D* fHitEvsRecM; // Histogram
	47	TH2D* fHitEvsRecE; // Histogram
	48	TH1D* fDiffE; // Histogram
	49	TH2D* fHitsVsRecM; // Histogram
	50	AliFMDEdepMap fMap;
8f6ee336	51	AliFMDFloatMap fEta;
	52	AliFMDFloatMap fPhi;
	53	AliFMDFloatMap fMult;
	54	Bool_t fPrimary;
	55	public:
bf000c32	56	//__________________________________________________________________
8f6ee336	57	TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
	58	{
	59	TArrayF bins(n+1);
	60	Float_t dp = n / TMath::Log10(max / min);
	61	Float_t pmin = TMath::Log10(min);
	62	bins[0] = min;
	63	for (Int_t i = 1; i < n+1; i++) {
	64	Float_t p = pmin + i / dp;
	65	bins[i] = TMath::Power(10, p);
	66	}
	67	return bins;
	68	}
bf000c32	69	//__________________________________________________________________
8f6ee336	70	DrawHitsRecs(Bool_t primary=kFALSE,
	71	Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
	72	Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
	73	{
	74	fPrimary = primary;
	75	AddLoad(kRecPoints);
8ec606c2	76	AddLoad(kHits);
15b17c89	77	AddLoad(kDigits);
8f6ee336	78	if (fPrimary) AddLoad(kKinematics);
	79	TArrayF eloss(MakeLogScale(n, emin, emax));
	80	TArrayF mults(m+1);
	81	mults[0] = mmin;
	82	for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
15b17c89	83
	84	fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
	85	n, emin, emax, 1025, -.5, 1024.5);
	86	fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
	87	fHitEvsAdc->SetYTitle("ADC");
	88
	89	fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
8f6ee336	90	eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
15b17c89	91	fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
	92	fHitEvsRecM->SetYTitle("M_{rec}");
	93
	94	fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
	95	n, emin, emax, n, emin, emax);
	96	fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
	97	fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
8f6ee336	98
15b17c89	99
	100	fDiffE = new TH1D("diffE",
	101	"#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
	102	1100, -1, 1.1);
	103	fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
	104
8f6ee336	105	Double_t omin = -.5;
	106	Double_t omax = 7.5;
	107	Int_t o = 8;
15b17c89	108	fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
8f6ee336	109	o, omin, omax, m, mmin, mmax);
15b17c89	110	fHitsVsRecM->SetXTitle("# of Hits");
15b17c89	111	fHitsVsRecM->SetYTitle("M_{rec}");
8f6ee336	112	}
bf000c32	113	//__________________________________________________________________
9b48326f	114	/** Begining of event
	115	@param ev Event number
	116	@return @c false on error */
8f6ee336	117	Bool_t Begin(Int_t ev)
	118	{
	119	fMap.Reset();
15b17c89	120	return AliFMDInput::Begin(ev);
8f6ee336	121	}
bf000c32	122	//__________________________________________________________________
8f6ee336	123	Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
	124	{
	125	// Cache the energy loss
	126	if (!hit) return kFALSE;
	127	UShort_t det = hit->Detector();
	128	Char_t rng = hit->Ring();
	129	UShort_t sec = hit->Sector();
	130	UShort_t str = hit->Strip();
	131	if (str > 511) {
	132	AliWarning(Form("Bad strip number %d in hit", str));
	133	return kTRUE;
	134	}
	135	if (fPrimary) {
	136	TParticle* kine = fStack->Particle(hit->Track());
	137	if (!kine) return kTRUE;
	138	if (!kine->IsPrimary()) return kTRUE;
	139	}
	140
	141	fMap(det, rng, sec, str).fEdep += hit->Edep();
	142	fMap(det, rng, sec, str).fN++;
	143	return kTRUE;
	144	}
15b17c89	145
	146	//__________________________________________________________________
	147	Bool_t ProcessDigit(AliFMDDigit* digit)
	148	{
	149	if (!digit) return kTRUE;
	150
	151	UShort_t det = digit->Detector();
	152	Char_t rng = digit->Ring();
	153	UShort_t sec = digit->Sector();
	154	UShort_t str = digit->Strip();
	155	if (str > 511) {
	156	AliWarning(Form("Bad strip number %d in digit", str));
	157	return kTRUE;
	158	}
	159	Double_t edep = fMap(det, rng, sec, str).fEdep;
	160	if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
	161
	162	return kTRUE;
	163	}
	164
bf000c32	165	//__________________________________________________________________
bf000c32	166	Bool_t ProcessRecPoint(AliFMDRecPoint* single)
8f6ee336	167	{
15b17c89	168	if (!single) return kTRUE;
bf000c32	169	UShort_t det = single->Detector();
	170	Char_t rng = single->Ring();
	171	UShort_t sec = single->Sector();
	172	UShort_t str = single->Strip();
	173	if (str > 511) {
	174	AliWarning(Form("Bad strip number %d in single", str));
15b17c89	175	return kTRUE;
	176	}
	177	Double_t edep = fMap(det, rng, sec, str).fEdep;
	178	Int_t nhit = fMap(det, rng, sec, str).fN;
	179	if (edep > 0) {
	180	fHitEvsRecM->Fill(edep, single->Particles());
	181	fHitEvsRecE->Fill(edep, single->Edep());
	182	fDiffE->Fill((single->Edep() - edep) / edep);
8f6ee336	183	}
15b17c89	184	if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
8f6ee336	185	return kTRUE;
8f6ee336	186	}
bf000c32	187	//__________________________________________________________________
8f6ee336	188	Bool_t Finish()
	189	{
	190	gStyle->SetPalette(1);
	191	gStyle->SetOptTitle(0);
	192	gStyle->SetCanvasColor(0);
	193	gStyle->SetCanvasBorderSize(0);
	194	gStyle->SetPadColor(0);
	195	gStyle->SetPadBorderSize(0);
	196
15b17c89	197	new TCanvas("c0", fHitEvsAdc->GetTitle());
	198	fHitEvsAdc->SetStats(kFALSE);
	199	fHitEvsAdc->Draw("COLZ");
	200
	201	new TCanvas("c1", fHitEvsRecM->GetTitle());
	202	fHitEvsRecM->SetStats(kFALSE);
	203	fHitEvsRecM->Draw("COLZ");
	204
	205	new TCanvas("c2", fHitEvsRecE->GetTitle());
	206	fHitEvsRecE->SetStats(kFALSE);
	207	fHitEvsRecE->Draw("COLZ");
	208
	209	new TCanvas("c3", fDiffE->GetTitle());
	210	fDiffE->Draw();
8f6ee336	211
15b17c89	212	new TCanvas("c4", fHitsVsRecM->GetTitle());
	213	fHitsVsRecM->SetStats(kFALSE);
	214	fHitsVsRecM->Draw("COLZ");
8f6ee336	215
	216	return kTRUE;
	217	}
	218
	219	ClassDef(DrawHitsRecs,0);
	220	};
	221
	222	//____________________________________________________________________
	223	//
	224	// EOF
	225	//