[u/mrichter/AliRoot.git] / PWG0 / trigger / AliGenTriggerMapSelector.cxx

/* $Id$ */

#include "AliGenTriggerMapSelector.h"

#include <TVector3.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TNtuple.h>
#include <TProfile.h>
#include <TParticle.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TGeoManager.h>
#include <TF1.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TLine.h>

#include <AliLog.h>
#include <AliESD.h>
#include <AliRunLoader.h>
#include <AliStack.h>
#include <AliGenEventHeader.h>
#include <AliHeader.h>

#include <AliITSgeom.h>
#include <AliITSLoader.h>
#include <AliITSdigitSPD.h>
#include <AliITSRecPoint.h>

#include "AliPWG0Helper.h"

//
//

ClassImp(AliGenTriggerMapSelector)

AliGenTriggerMapSelector::AliGenTriggerMapSelector() :
  AliSelectorRL(),
  fChipsFired(0),
  fTracklets(0)
{
  //
  // Constructor. Initialization of pointers
  //
}

AliGenTriggerMapSelector::~AliGenTriggerMapSelector()
{
  //
  // Destructor
  //

  // histograms are in the output list and deleted when the output
  // list is deleted by the TSelector dtor
}

void AliGenTriggerMapSelector::SlaveBegin(TTree *tree)
{
  AliSelectorRL::SlaveBegin(tree);

  fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
  fTracklets = new TNtuple("fTracklets", "", "vertex:chip1:chip2:dirty");
}

void AliGenTriggerMapSelector::Init(TTree* tree)
{
  // read the user objects

  AliSelectorRL::Init(tree);

  // enable only the needed branches
  if (tree)
  {
    tree->SetBranchStatus("*", 0);
    tree->SetBranchStatus("fTriggerMask", 1);
    tree->SetBranchStatus("fSPDVertex*", 1);
    tree->SetBranchStatus("fSPDMult*", 1);
  }
}

Bool_t AliGenTriggerMapSelector::Process(Long64_t entry)
{
  //
  // processing
  //

  if (AliSelectorRL::Process(entry) == kFALSE)
    return kFALSE;

  // Check prerequisites
  if (!fESD)
  {
    AliDebug(AliLog::kError, "ESD branch not available");
    return kFALSE;
  }

  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);

  if (!eventTriggered)
  {
    AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
    return kTRUE;
  }

  AliHeader* header = GetHeader();
  if (!header)
  {
    AliDebug(AliLog::kError, "Header not available");
    return kFALSE;
  }

  // get the MC vertex
  AliGenEventHeader* genHeader = header->GenEventHeader();
  TArrayF vtxMC(3);
  genHeader->PrimaryVertex(vtxMC);
  Float_t zVertex = vtxMC[2];

  AliStack* stack = GetStack();
  if (!stack)
  {
    AliDebug(AliLog::kError, "Stack not available");
    return kFALSE;
  }
  Int_t nPrim  = stack->GetNprimary();

  AliRunLoader* runLoader = GetRunLoader();
  if (!runLoader)
  {
    AliDebug(AliLog::kError, "runloader not available");
    return kFALSE;
  }

  // TDirectory::TContext restores the current directory is restored when the scope ends.
  // This helps around ROOT bug #26025 and is good behaviour anyway
  TDirectory::TContext context(0);
  AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
  loader->LoadDigits("READ");
  TTree* treeD = loader->TreeD();
  if (!treeD)
  {
    AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
    return kFALSE;
  }

  treeD->SetBranchStatus("*", 0);
  treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
  treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);

  TClonesArray* digits = 0;
  treeD->SetBranchAddress("ITSDigitsSPD", &digits);
  if (digits);
    digits->Clear();

  const Int_t startSPD = 0; //geom->GetStartSPD();
  const Int_t lastSPD  = 239; //geom->GetLastSPD();

  const AliMultiplicity* mult = fESD->GetMultiplicity();
  if (!mult)
  {
    AliDebug(AliLog::kError, "AliMultiplicity not available");
    return kFALSE;
  }

  // loop over tracklets
  for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
  {
    Float_t theta = mult->GetTheta(i);
    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));

    Int_t label = mult->GetLabel(i);
    if (label >= nPrim)
    {
      AliDebug(AliLog::kDebug, Form("Skipping particle %d because it is not a primary (label %d)", i, label));
      continue;
    }
    if (label < 0)
    {
      AliDebug(AliLog::kDebug, Form("Skipping particle %d because its label is negative (label %d)", i, label));
      continue;
    }

    Int_t chip[2];
    Int_t nClusters[2];

    for (Int_t j=0; j<2; ++j)
    {
      chip[j] = -1;
      nClusters[j] = 0;
    }

    // now find clusters that belong to this label
    // loop over modules (ladders)
    for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
    {
      Int_t currentLayer = 0;
      if (moduleIndex >= 80)
        currentLayer = 1;

      treeD->GetEvent(moduleIndex);

      // get number of digits in this module
      Int_t ndigitsInModule = digits->GetEntriesFast();

      // loop over digits in this module
      for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
      {
        AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);

        // check if digit belongs to this tracklet
        Bool_t belongs = kFALSE;
        for (Int_t iTrack = 0; iTrack < dp->GetNTracks(); ++iTrack)
        {
          //Printf("%d %d %d %d: %d =? %d", i, moduleIndex, iDig, iTrack, label, dp->GetTrack(iTrack));
          if (dp->GetTrack(iTrack) == label)
            belongs = kTRUE;
        }

        if (!belongs)
          continue;

        Int_t column = dp->GetCoord1();
        Int_t isChip = column / 32;

        //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);

        fChipsFired->Fill(moduleIndex, isChip);

        if (chip[currentLayer] == -1)
        {
          chip[currentLayer] = moduleIndex * 5 + isChip;
          //Printf("%d --> %d %d", label, moduleIndex, isChip);
          nClusters[currentLayer]++;
        }
        else
        {
          if (chip[currentLayer] != moduleIndex * 5 + isChip)
            nClusters[currentLayer]++;
        }
      }
    }

    AliDebug(AliLog::kDebug, Form("Tracklet %d fired by chip %d and %d", i, chip[0], chip[1]));

    Bool_t dirty = kFALSE;
    if (nClusters[0] == 0 || nClusters[1] == 0)
    {
      AliDebug(AliLog::kDebug, Form("Not enough clusters found for tracklet %d", i));
      continue;
    }
    if (nClusters[0] > 1 || nClusters[1] > 1)
    {
      AliDebug(AliLog::kDebug, Form("Too many clusters found for tracklet %d; L1: %d; L2: %d", i, nClusters[0], nClusters[1]));
      dirty = kTRUE;
    }
    fTracklets->Fill(zVertex, chip[0], chip[1], dirty);
  }

  return kTRUE;
}

Bool_t AliGenTriggerMapSelector::Notify()
{
  AliRunLoader* runLoader = GetRunLoader();

  if (runLoader)
  {
    AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
    if (loader)
    {
      loader->UnloadDigits();
      loader->UnloadRecPoints();
    }
  }

  return AliSelectorRL::Notify();
}

void AliGenTriggerMapSelector::SlaveTerminate()
{
  // The SlaveTerminate() function is called after all entries or objects
  // have been processed. When running with PROOF SlaveTerminate() is called
  // on each slave server.

  AliSelectorRL::SlaveTerminate();

  // Add the histograms to the output on each slave server
  if (!fOutput)
  {
    AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
    return;
  }

  fOutput->Add(fChipsFired);
  fOutput->Add(fTracklets);
}

void AliGenTriggerMapSelector::Terminate()
{
  // 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.

  AliSelectorRL::Terminate();

  fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
  fTracklets = dynamic_cast<TNtuple*> (fOutput->FindObject("fTracklets"));

  if (!fTracklets || !fChipsFired)
  {
    AliError("Histograms not available");
    return;
  }

  WriteHistograms();
}

void AliGenTriggerMapSelector::WriteHistograms(const char* filename)
{
  TFile* file = TFile::Open(filename, "RECREATE");

  fChipsFired->Write();
  fTracklets->Write();

  file->Close();
}

void AliGenTriggerMapSelector::ReadHistograms(const char* filename)
{
  TFile* file = TFile::Open(filename);

  if (!file)
    return;

  fTracklets  = dynamic_cast<TNtuple*> (file->Get("fTracklets"));
  fChipsFired  = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
}

void AliGenTriggerMapSelector::GenerateTriggerMap(Bool_t clean)
{
  Int_t nEntries = fTracklets->GetEntries();
  //nEntries = 10;

  TH2F* hist = new TH2F("hist", ";layer1;layer2", 400, -0.5, 399.5, 800, 399.5, 1199.5);

  Long64_t nSkipped = 0;
  for (Int_t i = 0; i < nEntries; ++i)
  {
    fTracklets->GetEntry(i);
    Float_t* vars = fTracklets->GetArgs();

    if (clean && (Bool_t) vars[3] != kFALSE)
    {
      nSkipped++;
      continue;
    }

    hist->Fill(vars[1], vars[2]);
  }
  if (clean)
    Printf("Skipped %lld dirty entries", nSkipped);

  hist->Draw("COLZ");

  for (Int_t i = 1; i <= hist->GetNbinsX(); ++i)
    for (Int_t j = 1; j <= hist->GetNbinsY(); ++j)
      if (hist->GetBinContent(i, j) > 0)
        Printf("(%d,%d) ", TMath::Nint(hist->GetXaxis()->GetBinCenter(i)), TMath::Nint(hist->GetYaxis()->GetBinCenter(j)));
}
Commit	Line	Data
a6ce57c4	1	/* $Id$ */
	2
	3	#include "AliGenTriggerMapSelector.h"
	4
	5	#include <TVector3.h>
	6	#include <TFile.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TNtuple.h>
	10	#include <TProfile.h>
	11	#include <TParticle.h>
	12	#include <TCanvas.h>
	13	#include <TTree.h>
	14	#include <TGeoManager.h>
	15	#include <TF1.h>
	16	#include <TGraph.h>
	17	#include <TLegend.h>
	18	#include <TLine.h>
	19
	20	#include <AliLog.h>
	21	#include <AliESD.h>
	22	#include <AliRunLoader.h>
	23	#include <AliStack.h>
	24	#include <AliGenEventHeader.h>
	25	#include <AliHeader.h>
	26
	27	#include <AliITSgeom.h>
	28	#include <AliITSLoader.h>
	29	#include <AliITSdigitSPD.h>
	30	#include <AliITSRecPoint.h>
	31
	32	#include "AliPWG0Helper.h"
	33
	34	//
	35	//
	36
	37	ClassImp(AliGenTriggerMapSelector)
	38
	39	AliGenTriggerMapSelector::AliGenTriggerMapSelector() :
	40	AliSelectorRL(),
	41	fChipsFired(0),
	42	fTracklets(0)
	43	{
	44	//
	45	// Constructor. Initialization of pointers
	46	//
	47	}
	48
	49	AliGenTriggerMapSelector::~AliGenTriggerMapSelector()
	50	{
	51	//
	52	// Destructor
	53	//
	54
	55	// histograms are in the output list and deleted when the output
	56	// list is deleted by the TSelector dtor
	57	}
	58
	59	void AliGenTriggerMapSelector::SlaveBegin(TTree *tree)
	60	{
	61	AliSelectorRL::SlaveBegin(tree);
	62
	63	fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
	64	fTracklets = new TNtuple("fTracklets", "", "vertex:chip1:chip2:dirty");
65	}
66
67	void AliGenTriggerMapSelector::Init(TTree* tree)
68	{
69	// read the user objects
70
71	AliSelectorRL::Init(tree);
72
73	// enable only the needed branches
74	if (tree)
75	{
76	tree->SetBranchStatus("*", 0);
77	tree->SetBranchStatus("fTriggerMask", 1);
78	tree->SetBranchStatus("fSPDVertex*", 1);
79	tree->SetBranchStatus("fSPDMult*", 1);
80	}
81	}
82
83	Bool_t AliGenTriggerMapSelector::Process(Long64_t entry)
84	{
85	//
86	// processing
87	//
88
89	if (AliSelectorRL::Process(entry) == kFALSE)
90	return kFALSE;
91
92	// Check prerequisites
93	if (!fESD)
94	{
95	AliDebug(AliLog::kError, "ESD branch not available");
96	return kFALSE;
97	}
98
99	Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
100
101	if (!eventTriggered)
102	{
103	AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
104	return kTRUE;
105	}
106
107	AliHeader* header = GetHeader();
108	if (!header)
109	{
110	AliDebug(AliLog::kError, "Header not available");
111	return kFALSE;
112	}
113
114	// get the MC vertex
115	AliGenEventHeader* genHeader = header->GenEventHeader();
116	TArrayF vtxMC(3);
117	genHeader->PrimaryVertex(vtxMC);
118	Float_t zVertex = vtxMC[2];
119
120	AliStack* stack = GetStack();
121	if (!stack)
122	{
123	AliDebug(AliLog::kError, "Stack not available");
124	return kFALSE;
125	}
126	Int_t nPrim = stack->GetNprimary();
127
128	AliRunLoader* runLoader = GetRunLoader();
129	if (!runLoader)
130	{
131	AliDebug(AliLog::kError, "runloader not available");
132	return kFALSE;
133	}
134
135	// TDirectory::TContext restores the current directory is restored when the scope ends.
136	// This helps around ROOT bug #26025 and is good behaviour anyway
137	TDirectory::TContext context(0);
138	AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
139	loader->LoadDigits("READ");
140	TTree* treeD = loader->TreeD();
141	if (!treeD)
142	{
143	AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
144	return kFALSE;
145	}
146
147	treeD->SetBranchStatus("*", 0);
148	treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
149	treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
150
151	TClonesArray* digits = 0;
152	treeD->SetBranchAddress("ITSDigitsSPD", &digits);
153	if (digits);
154	digits->Clear();
155
156	const Int_t startSPD = 0; //geom->GetStartSPD();
157	const Int_t lastSPD = 239; //geom->GetLastSPD();
158
159	const AliMultiplicity* mult = fESD->GetMultiplicity();
160	if (!mult)
161	{
162	AliDebug(AliLog::kError, "AliMultiplicity not available");
163	return kFALSE;
164	}
165
166	// loop over tracklets
167	for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
168	{
169	Float_t theta = mult->GetTheta(i);
170	Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
171
172	Int_t label = mult->GetLabel(i);
173	if (label >= nPrim)
174	{
175	AliDebug(AliLog::kDebug, Form("Skipping particle %d because it is not a primary (label %d)", i, label));
176	continue;
177	}
178	if (label < 0)
179	{
180	AliDebug(AliLog::kDebug, Form("Skipping particle %d because its label is negative (label %d)", i, label));
181	continue;
182	}
183
184	Int_t chip[2];
185	Int_t nClusters[2];
186
187	for (Int_t j=0; j<2; ++j)
188	{
189	chip[j] = -1;
190	nClusters[j] = 0;
191	}
192
193	// now find clusters that belong to this label
194	// loop over modules (ladders)
195	for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
196	{
197	Int_t currentLayer = 0;
198	if (moduleIndex >= 80)
199	currentLayer = 1;
200
201	treeD->GetEvent(moduleIndex);
202
203	// get number of digits in this module
204	Int_t ndigitsInModule = digits->GetEntriesFast();
205
206	// loop over digits in this module
207	for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
208	{
209	AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
210
211	// check if digit belongs to this tracklet
212	Bool_t belongs = kFALSE;
213	for (Int_t iTrack = 0; iTrack < dp->GetNTracks(); ++iTrack)
214	{
215	//Printf("%d %d %d %d: %d =? %d", i, moduleIndex, iDig, iTrack, label, dp->GetTrack(iTrack));
216	if (dp->GetTrack(iTrack) == label)
217	belongs = kTRUE;
218	}
219
220	if (!belongs)
221	continue;
222
223	Int_t column = dp->GetCoord1();
224	Int_t isChip = column / 32;
225
226	//printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
227
228	fChipsFired->Fill(moduleIndex, isChip);
229
230	if (chip[currentLayer] == -1)
231	{
232	chip[currentLayer] = moduleIndex * 5 + isChip;
9cb768f2	233	//Printf("%d --> %d %d", label, moduleIndex, isChip);
a6ce57c4	234	nClusters[currentLayer]++;
	235	}
	236	else
	237	{
	238	if (chip[currentLayer] != moduleIndex * 5 + isChip)
	239	nClusters[currentLayer]++;
	240	}
	241	}
	242	}
	243
	244	AliDebug(AliLog::kDebug, Form("Tracklet %d fired by chip %d and %d", i, chip[0], chip[1]));
	245
	246	Bool_t dirty = kFALSE;
	247	if (nClusters[0] == 0 \|\| nClusters[1] == 0)
	248	{
	249	AliDebug(AliLog::kDebug, Form("Not enough clusters found for tracklet %d", i));
	250	continue;
	251	}
	252	if (nClusters[0] > 1 \|\| nClusters[1] > 1)
	253	{
	254	AliDebug(AliLog::kDebug, Form("Too many clusters found for tracklet %d; L1: %d; L2: %d", i, nClusters[0], nClusters[1]));
	255	dirty = kTRUE;
	256	}
	257	fTracklets->Fill(zVertex, chip[0], chip[1], dirty);
	258	}
	259
	260	return kTRUE;
	261	}
	262
	263	Bool_t AliGenTriggerMapSelector::Notify()
	264	{
	265	AliRunLoader* runLoader = GetRunLoader();
	266
	267	if (runLoader)
	268	{
	269	AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
	270	if (loader)
	271	{
	272	loader->UnloadDigits();
	273	loader->UnloadRecPoints();
	274	}
	275	}
	276
	277	return AliSelectorRL::Notify();
	278	}
	279
	280	void AliGenTriggerMapSelector::SlaveTerminate()
	281	{
	282	// The SlaveTerminate() function is called after all entries or objects
	283	// have been processed. When running with PROOF SlaveTerminate() is called
	284	// on each slave server.
	285
	286	AliSelectorRL::SlaveTerminate();
	287
	288	// Add the histograms to the output on each slave server
	289	if (!fOutput)
	290	{
	291	AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
	292	return;
	293	}
	294
	295	fOutput->Add(fChipsFired);
	296	fOutput->Add(fTracklets);
	297	}
298
299	void AliGenTriggerMapSelector::Terminate()
300	{
301	// The Terminate() function is the last function to be called during
302	// a query. It always runs on the client, it can be used to present
303	// the results graphically or save the results to file.
304
305	AliSelectorRL::Terminate();
306
307	fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
308	fTracklets = dynamic_cast<TNtuple*> (fOutput->FindObject("fTracklets"));
309
310	if (!fTracklets \|\| !fChipsFired)
311	{
312	AliError("Histograms not available");
313	return;
314	}
315
316	WriteHistograms();
317	}
318
319	void AliGenTriggerMapSelector::WriteHistograms(const char* filename)
320	{
321	TFile* file = TFile::Open(filename, "RECREATE");
322
323	fChipsFired->Write();
324	fTracklets->Write();
325
326	file->Close();
327	}
328
329	void AliGenTriggerMapSelector::ReadHistograms(const char* filename)
330	{
331	TFile* file = TFile::Open(filename);
332
333	if (!file)
334	return;
335
336	fTracklets = dynamic_cast<TNtuple*> (file->Get("fTracklets"));
337	fChipsFired = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
338	}
9cb768f2	339
	340	void AliGenTriggerMapSelector::GenerateTriggerMap(Bool_t clean)
	341	{
	342	Int_t nEntries = fTracklets->GetEntries();
	343	//nEntries = 10;
	344
	345	TH2F* hist = new TH2F("hist", ";layer1;layer2", 400, -0.5, 399.5, 800, 399.5, 1199.5);
	346
	347	Long64_t nSkipped = 0;
	348	for (Int_t i = 0; i < nEntries; ++i)
	349	{
	350	fTracklets->GetEntry(i);
	351	Float_t* vars = fTracklets->GetArgs();
	352
	353	if (clean && (Bool_t) vars[3] != kFALSE)
	354	{
	355	nSkipped++;
	356	continue;
	357	}
	358
	359	hist->Fill(vars[1], vars[2]);
	360	}
	361	if (clean)
	362	Printf("Skipped %lld dirty entries", nSkipped);
	363
	364	hist->Draw("COLZ");
	365
	366	for (Int_t i = 1; i <= hist->GetNbinsX(); ++i)
	367	for (Int_t j = 1; j <= hist->GetNbinsY(); ++j)
	368	if (hist->GetBinContent(i, j) > 0)
	369	Printf("(%d,%d) ", TMath::Nint(hist->GetXaxis()->GetBinCenter(i)), TMath::Nint(hist->GetYaxis()->GetBinCenter(j)));
	370	}