]>
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 | } |