Removing dummy if
[u/mrichter/AliRoot.git] / PWGUD / selectors / trigger / AliGenTriggerMapSelector.cxx
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;
233           //Printf("%d --> %d %d", label, moduleIndex, isChip);
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 }
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 }