Transition PWG0 -> PWGUD
[u/mrichter/AliRoot.git] / PWGUD / selectors / trigger / AliGenTriggerMapSelector.cxx
CommitLineData
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
37ClassImp(AliGenTriggerMapSelector)
38
39AliGenTriggerMapSelector::AliGenTriggerMapSelector() :
40 AliSelectorRL(),
41 fChipsFired(0),
42 fTracklets(0)
43{
44 //
45 // Constructor. Initialization of pointers
46 //
47}
48
49AliGenTriggerMapSelector::~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
59void 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
67void 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
83Bool_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
263Bool_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
280void 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
299void 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
319void 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
329void 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
340void 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}