]>
Commit | Line | Data |
---|---|---|
0fa651e5 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //------------------------------------------------------------------------- | |
17 | // Implementation of class AliAnalysisTaskPileup | |
18 | // | |
19 | // The class re-weights the number of analysed events with the correction | |
20 | // factor to account for the event pileup. | |
21 | // | |
22 | // The final idea will be to include all of the methods for the pileup | |
23 | // calculation. | |
24 | // | |
25 | // For the moment, the implemented method is based on a work | |
26 | // by L. Aphecetche. | |
27 | // The probability of having at least 1 collision is given by: | |
28 | // MB L0 trigger / max collision rate | |
29 | // The max collision rate is given by: | |
30 | // number of bunches X revolution frequency ~ CBEAMB L0 trigger | |
31 | // Assuming a poissonian distribution of the collision, with: | |
32 | // mu = mean number of collisions | |
33 | // the probability of having 0 collisions is: | |
34 | // exp{-mu} = 1 - L0b_MB / max collision rate | |
35 | // Since the MB trigger is measuring the probability of having at least | |
36 | // 1 collisions, the number should be corrected by the factor: | |
37 | // CF = mu / ( 1 - exp{-mu} ) | |
38 | // The class weights the number of events with this correction factor | |
39 | // | |
40 | // CAVEATS: | |
41 | // - so far the class needs to access the OCDB | |
42 | // (hopefully it will be changed in the future) | |
43 | // Hence the following libraries need to be included (in addition to | |
44 | // the standard analysis ones): | |
45 | // libRAWDatabase.so libCDB.so libSTEER.so libPWG3base.so | |
46 | //------------------------------------------------------------------------- | |
47 | ||
48 | #include <Riostream.h> | |
49 | ||
50 | // ROOT includes | |
51 | #include "TH1.h" | |
52 | #include "TH2.h" | |
53 | #include "TCanvas.h" | |
54 | #include "TROOT.h" | |
55 | #include "TString.h" | |
56 | #include "TObjArray.h" | |
57 | #include "TObjString.h" | |
58 | #include "TMath.h" | |
59 | ||
60 | // STEER includes | |
61 | #include "AliESDEvent.h" | |
62 | #include "AliESDInputHandler.h" | |
63 | #include "AliAODEvent.h" | |
64 | #include "AliAODInputHandler.h" | |
65 | #include "AliTriggerScalersESD.h" | |
66 | #include "AliLog.h" | |
67 | ||
68 | // ANALYSIS includes | |
69 | #include "AliAnalysisTaskSE.h" | |
70 | #include "AliAnalysisManager.h" | |
71 | #include "AliAnalysisDataSlot.h" | |
72 | #include "AliAnalysisDataContainer.h" | |
73 | ||
74 | #include "AliCounterCollection.h" | |
75 | ||
76 | #include "AliAnalysisTaskPileup.h" | |
77 | ||
78 | #if defined(READOCDB) | |
79 | #include "AliTriggerRunScalers.h" | |
80 | #include "AliTriggerConfiguration.h" | |
81 | #include "AliTriggerClass.h" | |
82 | #include "AliCDBManager.h" | |
83 | #include "AliCDBEntry.h" | |
84 | #include "AliCDBPath.h" | |
85 | #endif | |
86 | ||
87 | ||
88 | ClassImp(AliAnalysisTaskPileup) | |
89 | ||
90 | ||
91 | //________________________________________________________________________ | |
92 | AliAnalysisTaskPileup::AliAnalysisTaskPileup(const char *name) : | |
93 | AliAnalysisTaskSE(name), | |
94 | fEventCounters(0x0), | |
95 | fHistoEventsList(0x0), | |
96 | fTriggerClasses(0x0), | |
97 | fTriggerClassIndex(new TArrayI(50)) | |
98 | #if defined(READOCDB) | |
99 | , fTriggerRunScalers(0x0), | |
100 | fDefaultStorage(new TString()) | |
101 | #endif | |
102 | ||
103 | { | |
104 | /// Constructor | |
105 | ||
106 | DefineOutput(1,AliCounterCollection::Class()); | |
107 | DefineOutput(2,TObjArray::Class()); | |
108 | } | |
109 | ||
110 | //________________________________________________________________________ | |
111 | AliAnalysisTaskPileup::~AliAnalysisTaskPileup() | |
112 | { | |
113 | /// Destructor | |
114 | delete fEventCounters; | |
115 | delete fHistoEventsList; | |
116 | delete fTriggerClasses; | |
117 | delete fTriggerClassIndex; | |
118 | ||
119 | #if defined(READOCDB) | |
120 | delete fTriggerRunScalers; | |
121 | delete fDefaultStorage; | |
122 | #endif | |
123 | ||
124 | } | |
125 | ||
126 | ||
127 | //___________________________________________________________________________ | |
128 | void AliAnalysisTaskPileup::NotifyRun() | |
129 | { | |
130 | /// Notify run | |
131 | ||
132 | #if defined(READOCDB) | |
133 | if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) { | |
134 | AliCDBManager::Instance()->SetDefaultStorage(fDefaultStorage->Data()); | |
135 | } | |
136 | ||
137 | AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber()); | |
138 | ||
139 | AliCDBEntry *entry = 0x0; | |
140 | entry = AliCDBManager::Instance()->Get("GRP/CTP/Config"); | |
141 | if ( entry ) { | |
142 | AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject(); | |
143 | const TObjArray& classesArray = trigConf->GetClasses(); | |
144 | Int_t nclasses = classesArray.GetEntriesFast(); | |
145 | ||
146 | if ( fTriggerClasses ) delete fTriggerClasses; | |
147 | ||
148 | fTriggerClasses = new TObjArray(nclasses); | |
149 | fTriggerClasses->SetOwner(); | |
150 | ||
151 | Int_t currActive = -1; | |
152 | for( Int_t iclass=0; iclass < nclasses; iclass++ ) { | |
153 | AliTriggerClass* trclass = (AliTriggerClass*)classesArray.At(iclass); | |
154 | if (trclass && trclass->GetMask()>0) { | |
155 | currActive++; | |
156 | Int_t currPos = currActive; | |
157 | ||
158 | // Store the CBEAMB class at the first position | |
159 | TString trigName = trclass->GetName(); | |
160 | if ( trigName.Contains("CBEAMB") && ! trigName.Contains("WU") ) | |
161 | currPos = 0; | |
162 | else if ( ! fTriggerClasses->At(0) ) | |
163 | currPos++; | |
164 | ||
165 | Int_t trindex = TMath::Nint(TMath::Log2(trclass->GetMask())); | |
166 | TObjString* objString = new TObjString(trigName.Data()); | |
167 | fTriggerClasses->AddAtAndExpand(objString, currPos); | |
168 | (*fTriggerClassIndex)[currPos] = trindex; | |
169 | AliDebug(3, Form("Current class %s index %i position %i", trigName.Data(), trindex, currPos)); | |
170 | } // is class active | |
171 | } // loop on trigger classes | |
172 | } // if entry | |
173 | ||
174 | entry = AliCDBManager::Instance()->Get("GRP/CTP/Scalers"); | |
175 | if (entry) { | |
176 | AliInfo("Found an AliTriggerRunScalers in GRP/CTP/Scalers, reading it"); | |
177 | fTriggerRunScalers = dynamic_cast<AliTriggerRunScalers*> (entry->GetObject()); | |
178 | entry->SetOwner(0); | |
179 | if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow"); | |
180 | } | |
181 | #endif | |
182 | ||
183 | } | |
184 | ||
185 | ||
186 | //___________________________________________________________________________ | |
187 | void AliAnalysisTaskPileup::UserCreateOutputObjects() | |
188 | { | |
189 | /// Create histograms and counters | |
190 | ||
191 | // The framework has problems if the name of the object | |
192 | // and the one of container differ | |
193 | // To be complaint, get the name from container and set it | |
194 | TString containerName = GetOutputSlot(1)->GetContainer()->GetName(); | |
195 | ||
196 | // initialize event counters | |
197 | fEventCounters = new AliCounterCollection(containerName.Data()); | |
198 | fEventCounters->AddRubric("event", "any/correctedL0"); | |
199 | fEventCounters->AddRubric("trigger", 1000000); | |
200 | fEventCounters->AddRubric("selection", "any/hasVtxContrib/nonPileupSPD"); | |
201 | fEventCounters->Init(kTRUE); | |
202 | ||
203 | // Post data at least once per task to ensure data synchronisation (required for merging) | |
204 | PostData(1, fEventCounters); | |
205 | } | |
206 | ||
207 | //________________________________________________________________________ | |
208 | void AliAnalysisTaskPileup::UserExec(Option_t *) | |
209 | { | |
210 | /// Called for each event | |
211 | ||
212 | AliAODEvent* aodEvent = 0x0; | |
213 | ||
214 | AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent()); | |
215 | if ( ! esdEvent ) { | |
216 | aodEvent = dynamic_cast<AliAODEvent*> (InputEvent()); | |
217 | } | |
218 | ||
219 | if ( ! aodEvent && ! esdEvent ) { | |
220 | AliError ("AOD or ESD event not found. Nothing done!"); | |
221 | return; | |
222 | } | |
223 | ||
224 | // check physics selection | |
225 | //Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected()); | |
226 | ||
227 | #if defined(READOCDB) | |
228 | ||
229 | TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses(); | |
230 | ||
231 | if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events | |
232 | ||
233 | AliDebug(2, Form("\nEvent %lli\n", Entry())); | |
234 | ||
235 | Int_t nPoints = fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast(); | |
236 | ||
237 | AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber()); | |
238 | Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp); | |
239 | if ( position < 0 ) { | |
240 | AliWarning("Position out of range: put to 1"); | |
241 | position = 1; | |
242 | } | |
243 | if ( position == 0 ) position++; // Don't trust the first one | |
244 | else if ( position + 1 >= nPoints ) position--; | |
245 | AliDebug(2, Form("position %i\n", position)); | |
246 | AliTriggerScalersRecordESD* trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position); | |
247 | AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0; | |
248 | ||
249 | // Sometimes scalers are filled very close to each others | |
250 | // in this case skip to the next entry | |
251 | for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) { | |
252 | trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos); | |
253 | Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() ); | |
254 | AliDebug(2, Form("Pos %i TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime)); | |
255 | ||
256 | if ( deltaTime > 1 ) | |
257 | break; | |
258 | } | |
259 | ||
260 | #endif | |
261 | ||
262 | ULong64_t trigMask = 0; | |
263 | Double_t correctFactor = 1.; | |
264 | ||
265 | Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors(); | |
266 | Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8); | |
267 | ||
268 | TString selKey[3] = {"any","hasVtxContrib","nonPileupSPD"}; | |
269 | Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )}; | |
270 | ||
271 | //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT | |
272 | ||
273 | TString trigName = ""; | |
274 | TString eventType = ""; | |
275 | ||
276 | Int_t nTriggerClasses = fTriggerClasses->GetEntries(); | |
277 | Int_t classIndex = -1; | |
278 | Double_t deltaScalersBeam = 0., deltaScalers = 0.; | |
279 | Bool_t isFiredOnce = kFALSE; | |
280 | for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) { | |
281 | ||
282 | Double_t correctFactorL0 = 1.; | |
283 | ||
284 | Bool_t isClassFired = kFALSE; | |
285 | ||
286 | if ( itrig < nTriggerClasses ) { | |
287 | ||
288 | // Check if current mask contains trigger | |
289 | trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString(); | |
290 | classIndex = (*fTriggerClassIndex)[itrig]; | |
291 | trigMask = ( 1ull << classIndex ); | |
292 | isClassFired = ( trigMask & InputEvent()->GetTriggerMask() ); | |
293 | ||
294 | if ( isClassFired || itrig == 0 ) { | |
295 | // Get scalers | |
296 | const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1); | |
297 | const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1); | |
298 | deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB(); | |
299 | ||
300 | if ( itrig == 0 ) | |
301 | deltaScalersBeam = deltaScalers; | |
302 | else if ( isClassFired ) { | |
303 | correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam); | |
304 | AliDebug(2, Form("Scalers: %s %.0f %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0)); | |
305 | } | |
306 | } | |
307 | } | |
308 | else { | |
309 | trigName = "any"; | |
310 | classIndex = -1; | |
311 | isClassFired = isFiredOnce; | |
312 | } | |
313 | ||
314 | if ( ! isClassFired ) continue; | |
315 | isFiredOnce = kTRUE; | |
316 | ||
317 | //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT | |
318 | //if ( classIndex > 1 ) printf("Index: trigger %i scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT | |
319 | ||
320 | AliDebug(2, Form("Fired trig %s\n", trigName.Data())); | |
321 | ||
322 | for ( Int_t ievType=0; ievType<2; ievType++ ){ | |
323 | switch ( ievType ) { | |
324 | #if defined(READOCDB) | |
325 | case kHeventsCorrectL0: | |
326 | correctFactor = correctFactorL0; | |
327 | eventType = "correctedL0"; | |
328 | break; | |
329 | #endif | |
330 | default: | |
331 | correctFactor = 1.; | |
332 | eventType = "any"; | |
333 | } | |
334 | ||
335 | for ( Int_t isel=0; isel<3; isel++ ) { | |
336 | if ( ! fillSel[isel] ) continue; | |
337 | fEventCounters->Count(Form("event:%s/trigger:%s/selection:%s",eventType.Data(),trigName.Data(), selKey[isel].Data()),correctFactor); | |
338 | } // loop on vertex selection | |
339 | } // loop on event type | |
340 | } // loop on trigger classes | |
341 | ||
342 | // Post final data. It will be written to a file with option "RECREATE" | |
343 | PostData(1, fEventCounters); | |
344 | } | |
345 | ||
346 | ||
347 | //________________________________________________________________________ | |
348 | void AliAnalysisTaskPileup::Terminate(Option_t *) | |
349 | { | |
350 | // | |
351 | /// Save final histograms | |
352 | /// and draw result to the screen | |
353 | // | |
354 | ||
355 | fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1)); | |
356 | if ( ! fEventCounters ) return; | |
357 | ||
358 | if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(2); | |
359 | fHistoEventsList->SetOwner(); | |
360 | ||
361 | TH2D* histo = 0x0; | |
362 | histo = fEventCounters->Draw("trigger","selection","event:any"); | |
363 | if ( histo ) { | |
364 | histo->SetName("hEvents"); | |
365 | histo->SetTitle("Events per trigger"); | |
366 | fHistoEventsList->AddAtAndExpand(histo, kHevents); | |
367 | } | |
368 | ||
369 | histo = fEventCounters->Draw("trigger","selection","event:correctedL0"); | |
370 | if ( histo ) { | |
371 | histo->SetName("hEventsCorrectL0"); | |
372 | histo->SetTitle("L0 corrected events"); | |
373 | fHistoEventsList->AddAtAndExpand(histo, kHeventsCorrectL0); | |
374 | } | |
375 | ||
376 | PostData(2, fHistoEventsList); | |
377 | ||
378 | if ( gROOT->IsBatch() ) | |
379 | return; | |
380 | ||
381 | TH2D* histoPileupL0 = (TH2D*)fHistoEventsList->At(kHeventsCorrectL0)->Clone("hPileupL0"); | |
382 | histoPileupL0->Divide((TH2D*)fHistoEventsList->At(kHevents)); | |
383 | ||
384 | TCanvas *can = new TCanvas("can1_Pileup","Pileup",10,10,310,310); | |
385 | can->SetFillColor(10); can->SetHighLightColor(10); | |
386 | can->SetLeftMargin(0.15); can->SetBottomMargin(0.15); | |
387 | histoPileupL0->DrawCopy("text"); | |
388 | } | |
389 | ||
390 | ||
391 | //________________________________________________________________________ | |
392 | Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB) | |
393 | { | |
394 | /// Get the correction factor for L0 calculation | |
395 | ||
396 | if ( nCBEAMB == 0. ) | |
397 | return 1.; | |
398 | ||
399 | Double_t ratio = nCINT1B / nCBEAMB; | |
400 | ||
401 | if ( ratio >= 1. || ratio == 0. ) | |
402 | return 1.; | |
403 | ||
404 | Double_t mu = -TMath::Log(1-ratio); | |
405 | ||
406 | return mu / ( 1. - TMath::Exp(-mu) ); | |
407 | ||
408 | } | |
409 |