]>
Commit | Line | Data |
---|---|---|
9608df41 | 1 | /* $Id$ */ |
2 | ||
3 | #include "AliHighMultiplicitySelector.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> | |
dd367a14 | 19 | #include <TMath.h> |
b3b260d1 | 20 | #include <TLegend.h> |
21 | #include <TText.h> | |
9608df41 | 22 | |
23 | #include <AliLog.h> | |
24 | #include <AliESD.h> | |
25 | #include <AliRunLoader.h> | |
26 | #include <AliStack.h> | |
27 | ||
28 | #include <AliITSgeom.h> | |
29 | #include <AliITSLoader.h> | |
30 | #include <AliITSdigitSPD.h> | |
31 | #include <AliITSRecPoint.h> | |
32 | ||
33 | #include "AliPWG0Helper.h" | |
34 | ||
35 | // | |
a9aed06c | 36 | // Selector that produces plots needed for the high multiplicity analysis with the |
37 | // pixel detector | |
38 | // Needs cluster file | |
9608df41 | 39 | // |
40 | ||
41 | ClassImp(AliHighMultiplicitySelector) | |
42 | ||
43 | AliHighMultiplicitySelector::AliHighMultiplicitySelector() : | |
44 | AliSelectorRL(), | |
45 | fChipsLayer1(0), | |
46 | fChipsLayer2(0), | |
47 | fL1vsL2(0), | |
48 | fMvsL1(0), | |
49 | fMvsL2(0), | |
50 | fChipsFired(0), | |
51 | fPrimaryL1(0), | |
52 | fPrimaryL2(0), | |
53 | fClusterZL1(0), | |
54 | fClusterZL2(0), | |
55 | fClvsL1(0), | |
56 | fClvsL2(0), | |
a9aed06c | 57 | fCentralRegion(kFALSE) |
9608df41 | 58 | { |
59 | // | |
60 | // Constructor. Initialization of pointers | |
61 | // | |
62 | } | |
63 | ||
64 | AliHighMultiplicitySelector::~AliHighMultiplicitySelector() | |
65 | { | |
66 | // | |
67 | // Destructor | |
68 | // | |
69 | ||
70 | // histograms are in the output list and deleted when the output | |
71 | // list is deleted by the TSelector dtor | |
72 | } | |
73 | ||
74 | void AliHighMultiplicitySelector::SlaveBegin(TTree *tree) | |
75 | { | |
a9aed06c | 76 | // create histograms |
77 | ||
9608df41 | 78 | AliSelectorRL::SlaveBegin(tree); |
79 | ||
80 | fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5); | |
81 | ||
82 | fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters"); | |
83 | fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters"); | |
84 | ||
85 | fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20); | |
86 | fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20); | |
87 | } | |
88 | ||
89 | void AliHighMultiplicitySelector::Init(TTree* tree) | |
90 | { | |
91 | // read the user objects | |
92 | ||
93 | AliSelectorRL::Init(tree); | |
94 | ||
95 | // enable only the needed branches | |
96 | if (tree) | |
97 | { | |
98 | tree->SetBranchStatus("*", 0); | |
99 | tree->SetBranchStatus("fTriggerMask", 1); | |
100 | ||
101 | /*if (fTree->GetCurrentFile()) | |
102 | { | |
103 | TString fileName(fTree->GetCurrentFile()->GetName()); | |
104 | fileName.ReplaceAll("AliESDs", "geometry"); | |
105 | ||
106 | // load geometry | |
107 | TGeoManager::Import(fileName); | |
108 | }*/ | |
109 | } | |
110 | } | |
111 | ||
112 | Bool_t AliHighMultiplicitySelector::Process(Long64_t entry) | |
113 | { | |
114 | // | |
115 | // processing | |
116 | // | |
117 | ||
118 | if (AliSelectorRL::Process(entry) == kFALSE) | |
119 | return kFALSE; | |
120 | ||
121 | // Check prerequisites | |
122 | if (!fESD) | |
123 | { | |
124 | AliDebug(AliLog::kError, "ESD branch not available"); | |
125 | return kFALSE; | |
126 | } | |
127 | ||
69b09e3b | 128 | Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1); |
9608df41 | 129 | |
130 | if (!eventTriggered) | |
131 | { | |
132 | AliDebug(AliLog::kDebug, "Event not triggered. Skipping."); | |
133 | return kTRUE; | |
134 | } | |
135 | ||
136 | AliStack* stack = GetStack(); | |
137 | if (!stack) | |
138 | { | |
139 | AliDebug(AliLog::kError, "Stack not available"); | |
140 | return kFALSE; | |
141 | } | |
142 | ||
143 | Int_t nPrim = stack->GetNprimary(); | |
144 | Int_t multiplicity21 = 0; | |
145 | Int_t multiplicity16 = 0; | |
146 | ||
147 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
148 | { | |
149 | AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc)); | |
150 | ||
151 | TParticle* particle = stack->Particle(iMc); | |
152 | ||
153 | if (!particle) | |
154 | { | |
155 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc)); | |
156 | continue; | |
157 | } | |
158 | ||
159 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
160 | continue; | |
161 | ||
a9aed06c | 162 | if (fCentralRegion) |
9608df41 | 163 | { |
164 | if (TMath::Abs(particle->Eta()) < 1.05) | |
165 | multiplicity21++; | |
166 | if (TMath::Abs(particle->Eta()) < 0.8) | |
167 | multiplicity16++; | |
168 | } | |
169 | else | |
170 | { | |
171 | if (TMath::Abs(particle->Eta()) < 2.1) | |
172 | multiplicity21++; | |
173 | if (TMath::Abs(particle->Eta()) < 1.6) | |
174 | multiplicity16++; | |
175 | } | |
176 | }// end of mc particle | |
177 | ||
178 | AliRunLoader* runLoader = GetRunLoader(); | |
179 | if (!runLoader) | |
180 | { | |
181 | AliDebug(AliLog::kError, "runloader not available"); | |
182 | return kFALSE; | |
183 | } | |
184 | ||
185 | // TDirectory::TContext restores the current directory is restored when the scope ends. | |
186 | // This helps around ROOT bug #26025 and is good behaviour anyway | |
187 | TDirectory::TContext context(0); | |
188 | AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" ); | |
189 | loader->LoadDigits("READ"); | |
190 | TTree* treeD = loader->TreeD(); | |
191 | if (!treeD) | |
192 | { | |
193 | AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS"); | |
194 | return kFALSE; | |
195 | } | |
196 | ||
197 | treeD->SetBranchStatus("*", 0); | |
198 | treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1); | |
199 | treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1); | |
200 | ||
201 | TClonesArray* digits = 0; | |
202 | treeD->SetBranchAddress("ITSDigitsSPD", &digits); | |
a6e0ebfe | 203 | if (digits) |
9608df41 | 204 | digits->Clear(); |
205 | ||
206 | // each value for both layers | |
207 | Int_t totalNumberOfFO[2]; | |
208 | Int_t chipsHitByPrimaries[2]; | |
209 | //Int_t chipsHitBySecondaries[2]; | |
210 | ||
211 | for (Int_t i=0; i<2; ++i) | |
212 | { | |
213 | totalNumberOfFO[i] = 0; | |
214 | chipsHitByPrimaries[i] = 0; | |
215 | //chipsHitBySecondaries[i] = 0; | |
216 | } | |
217 | ||
218 | Int_t startSPD = 0; //geom->GetStartSPD(); | |
219 | Int_t lastSPD = 239; //geom->GetLastSPD(); | |
220 | ||
221 | //printf("%d %d\n", startSPD, lastSPD); | |
222 | // for (Int_t l=0; l<240; ++l) { AliITSgeomTGeo::GetModuleId(l, i, j, k); printf("%d --> %d\n", l, i); } | |
223 | ||
224 | // loop over modules (ladders) | |
225 | for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++) | |
226 | { | |
a9aed06c | 227 | if (fCentralRegion) |
9608df41 | 228 | { |
229 | if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3) | |
230 | continue; | |
231 | } | |
232 | ||
233 | Int_t currentLayer = 0; | |
234 | if (moduleIndex >= 80) | |
235 | currentLayer = 1; | |
236 | ||
237 | treeD->GetEvent(moduleIndex); | |
238 | ||
239 | // get number of digits in this module | |
240 | Int_t ndigitsInModule = digits->GetEntriesFast(); | |
241 | ||
242 | // get number of digits in each chip of the module | |
243 | Int_t ndigitsInChip[5]; | |
244 | Bool_t hitByPrimary[5]; | |
245 | for( Int_t iChip=0; iChip<5; iChip++) | |
246 | { | |
247 | ndigitsInChip[iChip]=0; | |
248 | hitByPrimary[iChip] = kFALSE; | |
249 | } | |
250 | ||
251 | // loop over digits in this module | |
252 | for (Int_t iDig=0; iDig<ndigitsInModule; iDig++) | |
253 | { | |
254 | AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig); | |
255 | Int_t column = dp->GetCoord1(); | |
256 | Int_t isChip = column / 32; | |
257 | ||
258 | //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip); | |
259 | ||
260 | fChipsFired->Fill(moduleIndex, isChip); | |
261 | ||
262 | ndigitsInChip[isChip]++; | |
263 | ||
264 | Bool_t debug = kFALSE; | |
265 | ||
266 | // find out which particle caused this chip to fire | |
267 | // if we find at least one primary we consider this chip being fired by a primary | |
268 | for (Int_t track=0; track<10; ++track) | |
269 | { | |
270 | Int_t label = dp->GetTrack(track); | |
271 | ||
272 | if (label < 0) | |
273 | continue; | |
274 | ||
275 | if (debug) | |
276 | printf("track %d contains label %d\n", track, label); | |
277 | ||
278 | TParticle* particle = stack->Particle(label); | |
279 | ||
280 | if (!particle) | |
281 | { | |
282 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label)); | |
283 | continue; | |
284 | } | |
285 | ||
286 | if (debug) | |
287 | { | |
288 | particle->Print(); | |
289 | printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother()); | |
290 | } | |
291 | ||
292 | // TODO delta electrons should be traced back to their mother. this is e.g. solved in AliITSClusterFinderV2::CheckLabels2 | |
293 | while (particle->P() < 0.02 && particle->GetStatusCode() == 0 && particle->GetFirstMother() >= 0) | |
294 | { | |
295 | label = particle->GetFirstMother(); | |
296 | particle = stack->Particle(label); | |
297 | ||
298 | if (!particle) | |
299 | break; | |
300 | ||
301 | if (debug) | |
302 | { | |
303 | printf("-->\n"); | |
304 | printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother()); | |
305 | particle->Print(); | |
306 | } | |
307 | } | |
308 | ||
309 | if (!particle) | |
310 | { | |
311 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label)); | |
312 | continue; | |
313 | } | |
314 | ||
315 | if (label > nPrim) | |
316 | continue; | |
317 | ||
318 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
319 | continue; | |
320 | ||
321 | if (debug) | |
322 | printf("This was a primary (or delta-electron of a primary)!\n"); | |
323 | ||
324 | hitByPrimary[isChip] = kTRUE; | |
325 | } | |
326 | } | |
327 | ||
328 | // get number of FOs in the module | |
329 | for (Int_t ifChip=0; ifChip<5; ifChip++) | |
330 | if( ndigitsInChip[ifChip] >= 1 ) | |
331 | { | |
332 | totalNumberOfFO[currentLayer]++; | |
333 | if (hitByPrimary[ifChip]) | |
334 | { | |
335 | chipsHitByPrimaries[currentLayer]++; | |
336 | } | |
337 | //else | |
338 | // chipsHitBySecondaries[currentLayer]++; | |
339 | } | |
340 | } | |
341 | ||
342 | //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2); | |
343 | ||
344 | // now find clusters | |
345 | Int_t clustersLayer[2]; | |
346 | clustersLayer[0] = 0; | |
347 | clustersLayer[1] = 0; | |
348 | ||
349 | loader->LoadRecPoints("READ"); | |
350 | TTree* treeR = loader->TreeR(); | |
351 | if (!treeR) | |
352 | { | |
353 | AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS"); | |
354 | return kFALSE; | |
355 | } | |
356 | ||
357 | // TODO branches! | |
358 | //treeR->SetBranchStatus("*", 0); | |
359 | ||
360 | TClonesArray* itsClusters = 0; | |
361 | treeR->SetBranchAddress("ITSRecPoints", &itsClusters); | |
362 | ||
363 | Int_t nTreeEntries = treeR->GetEntries(); | |
364 | for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry) | |
365 | { | |
366 | if (!treeR->GetEvent(iEntry)) | |
367 | continue; | |
368 | ||
369 | Int_t nClusters = itsClusters->GetEntriesFast(); | |
370 | ||
371 | while(nClusters--) | |
372 | { | |
373 | AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters); | |
374 | ||
375 | if (cluster->GetLayer() == 0) | |
376 | { | |
377 | clustersLayer[0]++; | |
378 | fClusterZL1->Fill(cluster->GetZ()); | |
379 | } | |
380 | else if (cluster->GetLayer() == 1) | |
381 | { | |
382 | clustersLayer[1]++; | |
383 | fClusterZL2->Fill(cluster->GetZ()); | |
384 | } | |
385 | } | |
386 | } | |
387 | ||
388 | fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]); | |
389 | fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]); | |
390 | ||
391 | return kTRUE; | |
392 | } | |
393 | ||
394 | Bool_t AliHighMultiplicitySelector::Notify() | |
395 | { | |
a9aed06c | 396 | // get next ITS runloader |
397 | ||
9608df41 | 398 | AliRunLoader* runLoader = GetRunLoader(); |
399 | ||
400 | if (runLoader) | |
401 | { | |
402 | AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" ); | |
403 | if (loader) | |
404 | { | |
405 | loader->UnloadDigits(); | |
406 | loader->UnloadRecPoints(); | |
407 | } | |
408 | } | |
409 | ||
410 | return AliSelectorRL::Notify(); | |
411 | } | |
412 | ||
413 | void AliHighMultiplicitySelector::SlaveTerminate() | |
414 | { | |
415 | // The SlaveTerminate() function is called after all entries or objects | |
416 | // have been processed. When running with PROOF SlaveTerminate() is called | |
417 | // on each slave server. | |
418 | ||
419 | AliSelectorRL::SlaveTerminate(); | |
420 | ||
421 | // Add the histograms to the output on each slave server | |
422 | if (!fOutput) | |
423 | { | |
424 | AliDebug(AliLog::kError, "ERROR: Output list not initialized."); | |
425 | return; | |
426 | } | |
427 | ||
428 | fOutput->Add(fChipsFired); | |
429 | fOutput->Add(fPrimaryL1); | |
430 | fOutput->Add(fPrimaryL2); | |
431 | fOutput->Add(fClusterZL1); | |
432 | fOutput->Add(fClusterZL2); | |
433 | } | |
434 | ||
435 | void AliHighMultiplicitySelector::Terminate() | |
436 | { | |
437 | // The Terminate() function is the last function to be called during | |
438 | // a query. It always runs on the client, it can be used to present | |
439 | // the results graphically or save the results to file. | |
440 | ||
441 | AliSelectorRL::Terminate(); | |
442 | ||
443 | fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired")); | |
444 | fPrimaryL1 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL1")); | |
445 | fPrimaryL2 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL2")); | |
446 | fClusterZL1 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL1")); | |
447 | fClusterZL2 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL2")); | |
448 | ||
449 | if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2) | |
450 | { | |
451 | AliError("Histograms not available"); | |
452 | return; | |
453 | } | |
454 | ||
455 | WriteHistograms(); | |
456 | } | |
457 | ||
458 | void AliHighMultiplicitySelector::WriteHistograms(const char* filename) | |
459 | { | |
a9aed06c | 460 | // write the histograms to a file |
461 | ||
9608df41 | 462 | TFile* file = TFile::Open(filename, "RECREATE"); |
463 | ||
464 | fChipsFired->Write(); | |
465 | fPrimaryL1->Write(); | |
466 | fPrimaryL2->Write(); | |
467 | fClusterZL1->Write(); | |
468 | fClusterZL2->Write(); | |
469 | ||
470 | file->Close(); | |
471 | } | |
472 | ||
473 | void AliHighMultiplicitySelector::ReadHistograms(const char* filename) | |
474 | { | |
a9aed06c | 475 | // read the data from a file and fill histograms |
476 | ||
9608df41 | 477 | TFile* file = TFile::Open(filename); |
478 | ||
479 | if (!file) | |
480 | return; | |
481 | ||
482 | fPrimaryL1 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL1")); | |
483 | fPrimaryL2 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL2")); | |
484 | fChipsFired = dynamic_cast<TH2F*> (file->Get("fChipsFired")); | |
485 | fClusterZL1 = dynamic_cast<TH1F*> (file->Get("fClusterZL1")); | |
486 | fClusterZL2 = dynamic_cast<TH1F*> (file->Get("fClusterZL2")); | |
487 | ||
488 | #define MULT 1001, -0.5, 1000.5 | |
489 | #define BINNING_LAYER1 401, -0.5, 400.5 | |
490 | #define BINNING_LAYER2 801, -0.5, 800.5 | |
491 | ||
492 | fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1); | |
493 | fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2); | |
494 | ||
495 | fL1vsL2 = new TH2F("fL1vsL2", ";Fired Chips Layer 1;Fired Chips Layer 2", BINNING_LAYER1, BINNING_LAYER2); | |
496 | fMvsL1 = new TH2F("fMvsL1", ";true multiplicity;fired chips layer1", MULT, BINNING_LAYER1); | |
497 | fMvsL2 = new TH2F("fMvsL2", ";true multiplicity;fired chips layer2", MULT, BINNING_LAYER2); | |
498 | ||
499 | fClvsL1 = new TH2F("fClvsL1", ";clusters layer1;fired chips layer1", MULT, BINNING_LAYER1); | |
500 | fClvsL2 = new TH2F("fClvsL2", ";clusters layer2;fired chips layer2", MULT, BINNING_LAYER2); | |
501 | ||
502 | for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++) | |
503 | { | |
504 | fPrimaryL1->GetEvent(i); | |
505 | fPrimaryL2->GetEvent(i); | |
506 | ||
507 | Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0]; | |
508 | Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0]; | |
509 | ||
510 | Int_t totalNumberOfFO[2]; | |
511 | totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1]; | |
512 | totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1]; | |
513 | ||
514 | Int_t chipsHitByPrimaries[2]; | |
515 | chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2]; | |
516 | chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2]; | |
517 | ||
518 | Int_t clustersLayer[2]; | |
519 | clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3]; | |
520 | clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3]; | |
521 | ||
522 | fChipsLayer1->Fill(totalNumberOfFO[0]); | |
523 | fChipsLayer2->Fill(totalNumberOfFO[1]); | |
524 | ||
525 | fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]); | |
526 | ||
527 | fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]); | |
528 | fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]); | |
529 | ||
530 | fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]); | |
531 | fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]); | |
532 | } | |
533 | } | |
534 | ||
b3b260d1 | 535 | TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut, Int_t upperCut) |
9608df41 | 536 | { |
f6093269 | 537 | // |
538 | // returns the trigger efficiency as function of multiplicity with a given cut | |
539 | // | |
9608df41 | 540 | |
541 | //cut and multiply with x-section | |
f6093269 | 542 | TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001); |
9608df41 | 543 | //allEvents->Sumw2(); |
544 | ||
545 | //cut and multiply with x-section | |
b3b260d1 | 546 | TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, upperCut); |
9608df41 | 547 | //proj->Sumw2(); |
548 | ||
f6093269 | 549 | //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy(); |
550 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); | |
9608df41 | 551 | |
552 | // make probability distribution out of it | |
553 | // TODO binomial errors do not work??? weird... | |
554 | proj->Divide(proj, allEvents, 1, 1, "B"); | |
555 | ||
f6093269 | 556 | return proj; |
557 | } | |
558 | ||
559 | TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut) | |
560 | { | |
561 | // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed | |
562 | ||
563 | TH1* proj = GetTriggerEfficiency(multVsLayer, cut); | |
564 | ||
565 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); | |
9608df41 | 566 | |
567 | for (Int_t i=1; i<=proj->GetNbinsX(); i++) | |
568 | { | |
569 | if (i <= xSection->GetNbinsX()) | |
570 | { | |
571 | Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i); | |
572 | Double_t error = 0; | |
573 | ||
574 | if (value != 0) | |
575 | error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i)); | |
576 | ||
577 | proj->SetBinContent(i, value); | |
578 | proj->SetBinError(i, error); | |
579 | } | |
580 | else | |
581 | { | |
582 | proj->SetBinContent(i, 0); | |
583 | proj->SetBinError(i, 0); | |
584 | } | |
585 | } | |
586 | ||
f6093269 | 587 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); |
9608df41 | 588 | |
589 | return proj; | |
590 | } | |
591 | ||
f6093269 | 592 | void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL) |
9608df41 | 593 | { |
a9aed06c | 594 | // creates graphs |
595 | ||
f6093269 | 596 | TGraph* effGraph = new TGraph; |
597 | effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title)); | |
598 | TGraph* ratioGraph = new TGraph; | |
599 | ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title)); | |
600 | TGraph* totalGraph = new TGraph; | |
601 | totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title)); | |
602 | ||
603 | for (Int_t cut = 0; cut <= 300; cut+=50) | |
604 | { | |
605 | TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone"); | |
606 | ||
607 | //proj->Rebin(3); | |
608 | //proj->Scale(1.0 / 3); | |
609 | ||
610 | new TCanvas; proj->DrawCopy(); | |
9608df41 | 611 | |
f6093269 | 612 | Int_t limitBin = proj->GetNbinsX()+1; |
613 | while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5) | |
614 | limitBin--; | |
615 | ||
616 | Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin); | |
617 | ||
618 | effGraph->SetPoint(effGraph->GetN(), cut, limit); | |
619 | ||
620 | proj = GetXSectionCut(xSection, fMvsL, cut); | |
621 | ||
622 | Double_t ratio = 0; | |
623 | Double_t total = 0; | |
624 | if (proj->Integral(1, 1001) > 0) | |
625 | { | |
626 | ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001); | |
627 | total = proj->Integral(proj->FindBin(limit), 1001); | |
628 | } | |
629 | ||
630 | ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio); | |
631 | totalGraph->SetPoint(totalGraph->GetN(), cut, total); | |
632 | ||
633 | Printf("Cut at %d --> trigger eff. is > 0.5 for mult. >= %.2f. That is the case for %f of the triggered, %e of all events", cut, limit, ratio, total); | |
634 | } | |
635 | ||
636 | TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800); | |
637 | canvas->Divide(2, 2); | |
638 | ||
639 | canvas->cd(1); | |
640 | effGraph->Draw("A*"); | |
641 | ||
642 | for (Int_t i=8; i<=10; ++i) | |
643 | { | |
644 | TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i); | |
645 | line->Draw(); | |
646 | } | |
647 | ||
648 | canvas->cd(2); ratioGraph->Draw("A*"); | |
649 | canvas->cd(3); gPad->SetLogy(); totalGraph->Draw("A*"); | |
650 | ||
651 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
652 | } | |
653 | ||
654 | void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit) | |
655 | { | |
9608df41 | 656 | // relative x-section, once we have a collision |
a9aed06c | 657 | |
9608df41 | 658 | xSection->Scale(1.0 / xSection->Integral()); |
659 | ||
660 | TGraph* ratioGraph = new TGraph; | |
661 | ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit)); | |
662 | TGraph* totalGraph = new TGraph; | |
663 | totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit)); | |
664 | ||
665 | Double_t max = 0; | |
666 | Int_t bestCut = -1; | |
667 | Double_t bestRatio = -1; | |
668 | Double_t bestTotal = -1; | |
669 | Int_t fullCut = -1; | |
670 | Double_t fullRatio = -1; | |
671 | Double_t fullTotal = -1; | |
672 | ||
673 | fMvsL->Sumw2(); | |
674 | ||
675 | for (Int_t cut = 50; cut <= 300; cut+=2) | |
676 | { | |
677 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
678 | ||
679 | Double_t ratio = 0; | |
680 | Double_t total = 0; | |
681 | if (proj->Integral(1, 1000) > 0) | |
682 | { | |
683 | ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000); | |
684 | total = proj->Integral(limit, 1000); | |
685 | } | |
686 | ||
687 | max = TMath::Max(max, total); | |
688 | ||
689 | //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio); | |
690 | ||
691 | if (total < max * 0.9 && bestCut == -1) | |
692 | { | |
693 | bestCut = cut; | |
694 | bestRatio = ratio; | |
695 | bestTotal = total; | |
696 | } | |
697 | ||
698 | if (ratio == 1 && fullCut == -1) | |
699 | { | |
700 | fullCut = cut; | |
701 | fullRatio = ratio; | |
702 | fullTotal = total; | |
703 | } | |
704 | ||
705 | ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio); | |
706 | totalGraph->SetPoint(totalGraph->GetN(), cut, total); | |
707 | } | |
708 | ||
709 | if (bestCut != -1) | |
710 | printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio); | |
711 | if (fullCut != -1) | |
712 | printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio); | |
713 | ||
714 | TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400); | |
715 | ratioGraph->Draw("A*"); | |
716 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
717 | ||
718 | canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400); | |
719 | totalGraph->Draw("A*"); | |
720 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
721 | } | |
722 | ||
723 | void AliHighMultiplicitySelector::JPRPlots() | |
724 | { | |
a9aed06c | 725 | // more plots |
726 | ||
9608df41 | 727 | /* |
728 | ||
729 | gSystem->Load("libPWG0base"); | |
730 | .L AliHighMultiplicitySelector.cxx+ | |
731 | x = new AliHighMultiplicitySelector(); | |
732 | x->ReadHistograms("highmult_hijing100k.root"); | |
733 | x->JPRPlots(); | |
734 | ||
735 | */ | |
736 | ||
737 | // get x-sections | |
738 | TFile* file = TFile::Open("crosssectionEx.root"); | |
739 | if (!file) | |
740 | return; | |
741 | ||
742 | TH1* xSections[2]; | |
743 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
744 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
745 | ||
746 | for (Int_t i=0; i<2; ++i) | |
747 | { | |
748 | if (!xSections[i]) | |
749 | continue; | |
750 | ||
751 | TH1* xSection = xSections[i]; | |
752 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
753 | //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean | |
754 | //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean | |
755 | Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean | |
756 | ||
f6093269 | 757 | // limit is N times the mean |
758 | Int_t limit = (Int_t) (xSection->GetMean() * 10); | |
759 | ||
760 | // 10^28 lum --> 1.2 kHz | |
761 | // 10^31 lum --> 1200 kHz | |
762 | Float_t rate = 1200e3; | |
763 | ||
764 | // time in s | |
765 | Float_t lengthRun = 1e6; | |
766 | ||
9608df41 | 767 | xSection->SetStats(kFALSE); |
f6093269 | 768 | xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2"); |
9608df41 | 769 | xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5)); |
770 | xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 400 : 350); | |
f6093269 | 771 | //xSection->GetYaxis()->SetTitle("relative cross section"); |
9608df41 | 772 | xSection->GetYaxis()->SetTitleOffset(1.2); |
773 | ||
9608df41 | 774 | // relative x-section, once we have a collision |
775 | xSection->Scale(1.0 / xSection->Integral()); | |
776 | ||
777 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
778 | ||
779 | Double_t ratio = 0; | |
780 | Double_t total = 0; | |
781 | if (proj->Integral(1, 1000) > 0) | |
782 | { | |
783 | ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000); | |
784 | total = proj->Integral(limit, 1000); | |
785 | } | |
786 | ||
787 | printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio); | |
788 | ||
789 | TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600); | |
790 | canvas->SetLogy(); | |
791 | xSection->DrawCopy(); | |
792 | proj->SetLineColor(2); | |
793 | proj->SetStats(kFALSE); | |
794 | proj->DrawCopy("SAME"); | |
795 | ||
796 | TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3); | |
797 | legend->SetFillColor(0); | |
798 | legend->AddEntry(xSection, "no trigger"); | |
799 | legend->AddEntry(proj, Form("FO trigger > %d chips", cut)); | |
800 | legend->Draw(); | |
801 | ||
802 | TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2); | |
803 | line->SetLineWidth(2); | |
804 | line->Draw(); | |
805 | ||
806 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
807 | ||
808 | TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600); | |
f6093269 | 809 | //canvas2->SetTopMargin(0.05); |
810 | //canvas2->SetRightMargin(0.05); | |
9608df41 | 811 | canvas2->SetLogy(); |
f6093269 | 812 | xSection->DrawCopy("HIST"); |
9608df41 | 813 | |
f6093269 | 814 | TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3); |
9608df41 | 815 | legend2->SetFillColor(0); |
816 | legend2->AddEntry(xSection, "no trigger"); | |
817 | ||
818 | TH1* proj2 = (TH1*) proj->Clone("random"); | |
819 | proj2->Reset(); | |
f6093269 | 820 | // MB lengthRun s 100 Hz |
821 | Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000)); | |
9608df41 | 822 | proj2->FillRandom(proj, nTrigger); |
823 | ||
824 | TH1* proj3 = (TH1*) proj2->Clone("random_clone"); | |
825 | proj3->Divide(proj); | |
826 | proj3->Fit("pol0", "0", ""); | |
827 | proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0)); | |
828 | ||
f6093269 | 829 | /* |
9608df41 | 830 | proj2->DrawCopy("SAME"); |
831 | legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio))); | |
f6093269 | 832 | */ |
9608df41 | 833 | |
834 | proj2 = (TH1*) proj->Clone("random2"); | |
835 | proj2->Reset(); | |
f6093269 | 836 | // 10^31 lum --> 1200 kHz; lengthRun s |
837 | nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun); | |
9608df41 | 838 | proj2->FillRandom(proj, nTrigger); |
839 | ||
840 | proj3 = (TH1*) proj2->Clone("random_clone2"); | |
841 | proj3->Divide(proj); | |
842 | proj3->Fit("pol0", "0", ""); | |
843 | proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0)); | |
844 | ||
845 | proj2->SetLineColor(4); | |
f6093269 | 846 | proj2->SetMarkerStyle(7); |
847 | proj2->SetMarkerColor(4); | |
848 | proj2->DrawCopy("SAME P"); | |
849 | //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio))); | |
850 | legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut)); | |
9608df41 | 851 | |
852 | legend2->Draw(); | |
853 | line->Draw(); | |
854 | ||
855 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
f6093269 | 856 | canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); |
9608df41 | 857 | } |
858 | } | |
859 | ||
d1594f8a | 860 | void AliHighMultiplicitySelector::Ntrigger(Bool_t relative) |
f91bc12d | 861 | { |
862 | // | |
863 | // produces a spectrum created with N triggers | |
864 | // number of triggers and thresholds for the moment fixed | |
865 | // | |
866 | ||
867 | /* | |
868 | ||
dd367a14 | 869 | gSystem->Load("libANALYSIS"); |
f91bc12d | 870 | gSystem->Load("libPWG0base"); |
871 | .L AliHighMultiplicitySelector.cxx+g | |
872 | x = new AliHighMultiplicitySelector(); | |
873 | x->ReadHistograms("highmult_hijing100k.root"); | |
874 | x->Ntrigger(); | |
875 | ||
d1594f8a | 876 | gSystem->Load("libPWG0base"); |
877 | .L AliHighMultiplicitySelector.cxx+g | |
878 | x = new AliHighMultiplicitySelector(); | |
879 | x->ReadHistograms("highmult_hijing100k.root"); | |
880 | x->Ntrigger(kFALSE); | |
881 | ||
f91bc12d | 882 | */ |
883 | ||
884 | // get x-sections | |
885 | TFile* file = TFile::Open("crosssectionEx.root"); | |
886 | if (!file) | |
887 | return; | |
888 | ||
889 | TH1* xSections[2]; | |
890 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
891 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
892 | ||
51f6de65 | 893 | // 10^28 lum --> 1.4 kHz |
894 | // 10^31 lum --> 1400 kHz | |
895 | //Float_t rate = 1400e3; | |
896 | Float_t rate = 1.4e3; | |
f91bc12d | 897 | |
898 | // time in s | |
51f6de65 | 899 | Float_t lengthRun = 1e5; |
f91bc12d | 900 | |
901 | Int_t colors[] = { 2, 3, 4, 6, 7, 8 }; | |
902 | Int_t markers[] = { 7, 2, 4, 5, 6, 27 }; | |
903 | ||
904 | // put to 2 for second layer | |
905 | for (Int_t i=0; i<1; ++i) | |
906 | { | |
907 | if (!xSections[i]) | |
908 | continue; | |
909 | ||
910 | TH1* xSection = xSections[i]; | |
911 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
912 | ||
a9aed06c | 913 | //Int_t nCuts = 6; |
914 | //Int_t cuts[] = { 0, 164, 178, 190, 204, 216 }; | |
f91bc12d | 915 | //Int_t nCuts = 4; |
916 | //Int_t cuts[] = { 0, 164, 190, 216 }; | |
d1594f8a | 917 | |
918 | //Int_t nCuts = 4; | |
919 | //Int_t cuts[] = { 0, 114, 145, 165 }; | |
920 | //Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 }; | |
921 | ||
922 | Int_t nCuts = 3; | |
51f6de65 | 923 | Int_t cuts[] = { 0, 114, 148 }; |
924 | ||
925 | //Int_t nCuts = 3; | |
926 | //Int_t cuts[] = { 0, 126, 162 }; | |
927 | //Float_t ratePerTrigger[] = { 60, 20.0, 20.0 }; | |
a9aed06c | 928 | |
f91bc12d | 929 | // desired trigger rate in Hz |
51f6de65 | 930 | Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 }; |
f91bc12d | 931 | |
932 | xSection->SetStats(kFALSE); | |
933 | xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2"); | |
934 | xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5)); | |
935 | xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350); | |
936 | //xSection->GetYaxis()->SetTitle("relative cross section"); | |
937 | xSection->GetYaxis()->SetTitleOffset(1.2); | |
938 | ||
939 | // relative x-section, once we have a collision | |
940 | xSection->Scale(1.0 / xSection->Integral()); | |
941 | ||
942 | TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600); | |
943 | canvas2->SetTopMargin(0.05); | |
944 | canvas2->SetRightMargin(0.05); | |
945 | canvas2->SetLogy(); | |
d1594f8a | 946 | |
947 | if (relative) | |
948 | xSection->DrawCopy("HIST"); | |
f91bc12d | 949 | |
950 | TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4); | |
951 | legend2->SetFillColor(0); | |
952 | legend2->AddEntry(xSection, "cross-section"); | |
953 | ||
954 | for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut) | |
955 | { | |
956 | Int_t cut = cuts[currentCut]; | |
957 | ||
d1594f8a | 958 | TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); |
959 | ||
f91bc12d | 960 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); |
961 | ||
d1594f8a | 962 | Float_t triggerLimit = 0; |
963 | for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++) | |
964 | if (triggerEff->GetBinContent(bin) < 0.5) | |
965 | triggerLimit = triggerEff->GetXaxis()->GetBinCenter(bin); | |
966 | ||
967 | Printf("Efficiency limit (50%%) is at multiplicity %f", triggerLimit); | |
968 | ||
f91bc12d | 969 | Double_t total = 0; |
970 | if (proj->Integral(1, 1000) > 0) | |
971 | total = proj->Integral(1, 1000); | |
972 | ||
973 | printf("Cut at %d: rel. x-section = %e\n", cut, total); | |
974 | ||
975 | TH1* proj2 = (TH1*) proj->Clone("random2"); | |
976 | proj2->Reset(); | |
977 | ||
978 | // calculate downscale factor | |
979 | Float_t normalRate = rate * proj->Integral(1, 1000); | |
980 | Float_t downScale = normalRate / ratePerTrigger[currentCut]; | |
981 | if (downScale < 1) | |
982 | downScale = 1; | |
983 | Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun); | |
dd367a14 | 984 | nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000; |
f91bc12d | 985 | |
986 | Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger); | |
51f6de65 | 987 | if (nTrigger == 0) |
988 | continue; | |
989 | ||
f91bc12d | 990 | proj2->FillRandom(proj, nTrigger); |
991 | ||
992 | Int_t removed = 0; | |
993 | for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin) | |
994 | if (proj2->GetBinContent(bin) < 5) | |
995 | { | |
996 | removed += (Int_t) proj2->GetBinContent(bin); | |
997 | proj2->SetBinContent(bin, 0); | |
998 | } | |
999 | ||
1000 | Printf("Removed %d events", removed); | |
1001 | ||
d1594f8a | 1002 | if (relative) |
1003 | proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000)); | |
f91bc12d | 1004 | |
1005 | proj2->SetLineColor(colors[currentCut]); | |
1006 | proj2->SetMarkerStyle(markers[currentCut]); | |
1007 | proj2->SetMarkerColor(colors[currentCut]); | |
d1594f8a | 1008 | |
1009 | if (relative || currentCut > 0) { | |
1010 | proj2->DrawCopy("SAME P"); | |
1011 | } else | |
1012 | proj2->DrawCopy(" P"); | |
1013 | ||
745d6088 | 1014 | TString eventStr; |
1015 | if (nTrigger > 1e6) | |
1016 | { | |
1017 | eventStr.Form("%lld M", nTrigger / 1000 / 1000); | |
1018 | } | |
1019 | else if (nTrigger > 1e3) | |
1020 | { | |
1021 | eventStr.Form("%lld K", nTrigger / 1000); | |
1022 | } | |
1023 | else | |
1024 | eventStr.Form("%lld", nTrigger); | |
1025 | ||
1026 | TString triggerStr; | |
1027 | if (cut == 0) | |
1028 | { | |
1029 | triggerStr = "minimum bias"; | |
1030 | } | |
1031 | else | |
1032 | triggerStr.Form("FO > %d chips", cut); | |
d1594f8a | 1033 | |
745d6088 | 1034 | legend2->AddEntry(proj2, Form("%s evts, %s", eventStr.Data(), triggerStr.Data())); |
1035 | ||
1036 | if (triggerLimit > 1) | |
1037 | { | |
1038 | TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum()); | |
1039 | line->SetLineColor(colors[currentCut]); | |
1040 | line->Draw(); | |
1041 | } | |
f91bc12d | 1042 | } |
1043 | ||
1044 | legend2->Draw(); | |
1045 | ||
1046 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
1047 | canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); | |
1048 | } | |
1049 | } | |
1050 | ||
d1594f8a | 1051 | void AliHighMultiplicitySelector::Contamination() |
1052 | { | |
1053 | // | |
d1594f8a | 1054 | // |
1055 | ||
1056 | /* | |
1057 | ||
1058 | gSystem->Load("libANALYSIS"); | |
1059 | gSystem->Load("libPWG0base"); | |
1060 | .L AliHighMultiplicitySelector.cxx+g | |
1061 | x = new AliHighMultiplicitySelector(); | |
1062 | x->ReadHistograms("highmult_hijing100k.root"); | |
1063 | x->Contamination(); | |
1064 | ||
1065 | */ | |
1066 | ||
1067 | // get x-sections | |
1068 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1069 | if (!file) | |
1070 | return; | |
1071 | ||
1072 | TH1* xSections[2]; | |
1073 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1074 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1075 | ||
1076 | // rate = L * sigma | |
1077 | // sigma ~ 80 mb (Pythia 14 TeV) | |
1078 | // 10^28 lum --> 8e2 Hz | |
1079 | // 10^31 lum --> 8e5 Hz | |
1080 | Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 }; | |
1081 | ||
1082 | Int_t nCuts = 4; | |
1083 | Int_t cuts[] = { 104, 134, 154, 170 }; | |
1084 | ||
1085 | // put to 2 for second layer | |
1086 | for (Int_t i=0; i<1; ++i) | |
1087 | { | |
1088 | if (!xSections[i]) | |
1089 | continue; | |
1090 | ||
1091 | // relative x-section, once we have a collision | |
1092 | xSections[i]->Scale(1.0 / xSections[i]->Integral()); | |
1093 | ||
1094 | Int_t max = xSections[i]->GetNbinsX(); | |
1095 | max = 500; | |
1096 | ||
1097 | Float_t* xSection = new Float_t[max]; | |
1098 | for (Int_t mult = 0; mult < max; mult++) | |
1099 | xSection[mult] = xSections[i]->GetBinContent(mult+1); | |
1100 | ||
1101 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
1102 | ||
1103 | TGraph* graph = new TGraph; | |
1104 | ||
1105 | for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut) | |
1106 | { | |
1107 | Int_t cut = cuts[currentCut]; | |
1108 | Double_t rate = rates[currentCut]; | |
1109 | //Double_t rate = rates[3]; | |
1110 | ||
1111 | // coll. in 100 ns window | |
69b09e3b | 1112 | Double_t windowSize = 100e-9; |
1113 | //Double_t windowSize = 25e-9; | |
d1594f8a | 1114 | Double_t collPerWindow = windowSize * rate; |
1115 | Printf("coll/window = %f", collPerWindow); | |
1116 | Double_t windowsPerSecond = 1.0 / windowSize; | |
1117 | ||
1118 | TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); | |
1119 | Float_t* triggerEff = new Float_t[max]; | |
1120 | for (Int_t mult = 0; mult < max; mult++) | |
1121 | triggerEff[mult] = triggerEffHist->GetBinContent(mult+1); | |
1122 | ||
1123 | Double_t triggerRate = 0; | |
1124 | for (Int_t mult = 0; mult < max; mult++) | |
1125 | triggerRate += xSection[mult] * triggerEff[mult]; | |
1126 | ||
1127 | triggerRate *= TMath::Poisson(1, collPerWindow) * windowsPerSecond; | |
1128 | ||
1129 | Printf("Rate for 1 collision is %f Hz", triggerRate); | |
1130 | ||
1131 | Double_t triggerRate2 = 0; | |
1132 | for (Int_t mult = 0; mult < max; mult++) | |
1133 | for (Int_t mult2 = mult; mult2 < max; mult2++) | |
1134 | if (mult+mult2 < max) | |
1135 | triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2]; | |
1136 | ||
1137 | triggerRate2 *= TMath::Poisson(2, collPerWindow) * windowsPerSecond; | |
d1594f8a | 1138 | |
1139 | Printf("Rate for 2 collisions is %f Hz --> %.1f%%", triggerRate2, triggerRate2 / triggerRate * 100); | |
1140 | ||
1141 | Double_t triggerRate3 = 0; | |
144ff489 | 1142 | |
d1594f8a | 1143 | for (Int_t mult = 0; mult < max; mult++) |
1144 | for (Int_t mult2 = mult; mult2 < max-mult; mult2++) | |
1145 | for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++) | |
1146 | //if (mult+mult2+mult3 < max) | |
1147 | triggerRate3 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3]; | |
1148 | ||
1149 | triggerRate3 *= TMath::Poisson(3, collPerWindow) * windowsPerSecond; | |
1150 | //triggerRate3 *= collPerWindow * collPerWindow * rate; | |
1151 | ||
1152 | Printf("Rate for 3 collisions is %f Hz --> %.1f%%", triggerRate3, triggerRate3 / triggerRate * 100); | |
1153 | ||
1154 | Float_t totalContamination = (triggerRate2 + triggerRate3) / (triggerRate + triggerRate2 + triggerRate3); | |
1155 | ||
1156 | Printf("Total contamination is %.1f%%", totalContamination * 100); | |
1157 | ||
1158 | graph->SetPoint(graph->GetN(), cut, totalContamination); | |
1159 | ||
1160 | continue; | |
1161 | ||
1162 | Double_t triggerRate4 = 0; | |
1163 | for (Int_t mult = 0; mult < max; mult++) | |
1164 | for (Int_t mult2 = mult; mult2 < max-mult; mult2++) | |
1165 | for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++) | |
1166 | for (Int_t mult4 = 0; mult4 < max-mult-mult2-mult3; mult4++) | |
1167 | //if (mult+mult2+mult3+mult4 < max) | |
1168 | triggerRate4 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * xSection[mult4] * triggerEff[mult+mult2+mult3+mult4]; | |
1169 | ||
1170 | //triggerRate4 *= collPerWindow * collPerWindow * collPerWindow * rate; | |
1171 | triggerRate4 *= TMath::Poisson(4, collPerWindow) * windowsPerSecond; | |
1172 | ||
1173 | Printf("Rate for 4 collisions is %f Hz --> %.1f%%", triggerRate4, triggerRate4 / triggerRate * 100); | |
1174 | ||
1175 | // general code for n collisions follows, however much slower... | |
1176 | /* | |
1177 | const Int_t maxdepth = 4; | |
1178 | for (Int_t depth = 1; depth <= maxdepth; depth++) { | |
1179 | Double_t triggerRate = 0; | |
1180 | ||
1181 | Int_t m[maxdepth]; | |
1182 | for (Int_t d=0; d<maxdepth; d++) | |
1183 | m[d] = 0; | |
1184 | ||
1185 | while (m[0] < max) { | |
1186 | Double_t value = 1; | |
1187 | Int_t sum = 0; | |
1188 | for (Int_t d=0; d<depth; d++) { | |
1189 | value *= xSection[m[d]]; | |
1190 | sum += m[d]; | |
1191 | } | |
1192 | ||
1193 | if (sum < max) { | |
1194 | value *= triggerEff[sum]; | |
1195 | triggerRate += value; | |
1196 | } | |
1197 | ||
1198 | Int_t increase = depth-1; | |
1199 | ++m[increase]; | |
1200 | while (m[increase] == max && increase > 0) { | |
1201 | m[increase] = 0; | |
1202 | --increase; | |
1203 | ++m[increase]; | |
1204 | } | |
1205 | } | |
1206 | ||
1207 | triggerRate *= rate * TMath::Power(collPerWindow, depth - 1); | |
1208 | ||
1209 | Printf("Rate for %d collisions is %f Hz", depth, triggerRate); | |
1210 | }*/ | |
1211 | } | |
1212 | ||
1213 | new TCanvas; graph->Draw("AP*"); | |
1214 | } | |
1215 | } | |
1216 | ||
144ff489 | 1217 | void AliHighMultiplicitySelector::Contamination2() |
1218 | { | |
1219 | // | |
1220 | // produces a spectrum created with N triggers | |
1221 | // number of triggers and thresholds for the moment fixed | |
1222 | // | |
1223 | ||
1224 | /* | |
1225 | ||
1226 | gSystem->Load("libANALYSIS"); | |
1227 | gSystem->Load("libPWG0base"); | |
1228 | .L AliHighMultiplicitySelector.cxx+g | |
1229 | x = new AliHighMultiplicitySelector(); | |
1230 | x->ReadHistograms("highmult_hijing100k.root"); | |
1231 | x->Contamination2(); | |
1232 | ||
1233 | */ | |
1234 | ||
1235 | // get x-sections | |
1236 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1237 | if (!file) | |
1238 | return; | |
1239 | ||
1240 | TH1* xSections[2]; | |
1241 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1242 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1243 | ||
1244 | Int_t nCuts = 4; | |
1245 | Int_t cuts[] = { 104, 134, 154, 170 }; | |
1246 | ||
1247 | new TCanvas; | |
1248 | ||
1249 | Int_t colors[] = { 2, 3, 4, 6, 7, 8 }; | |
1250 | Int_t markers[] = { 7, 2, 4, 5, 6, 27 }; | |
1251 | ||
1252 | // put to 2 for second layer | |
1253 | for (Int_t i=0; i<1; ++i) | |
1254 | { | |
1255 | if (!xSections[i]) | |
1256 | continue; | |
1257 | ||
1258 | // relative x-section, once we have a collision | |
1259 | xSections[i]->Scale(1.0 / xSections[i]->Integral()); | |
1260 | ||
1261 | Int_t max = xSections[i]->GetNbinsX(); | |
1262 | max = 500; | |
1263 | ||
1264 | Float_t* xSection = new Float_t[max]; | |
1265 | for (Int_t mult = 0; mult < max; mult++) | |
1266 | xSection[mult] = xSections[i]->GetBinContent(mult+1); | |
1267 | ||
1268 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
1269 | ||
1270 | for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut) | |
1271 | { | |
1272 | TGraph* graph = new TGraph; | |
1273 | graph->SetMarkerColor(colors[currentCut]); | |
1274 | graph->SetMarkerStyle(markers[currentCut]); | |
1275 | ||
1276 | Int_t cut = cuts[currentCut]; | |
1277 | ||
1278 | TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); | |
1279 | Float_t* triggerEff = new Float_t[max]; | |
1280 | for (Int_t mult = 0; mult < max; mult++) | |
1281 | triggerEff[mult] = triggerEffHist->GetBinContent(mult+1); | |
1282 | ||
1283 | Double_t triggerRate = 0; | |
1284 | for (Int_t mult = 0; mult < max; mult++) | |
1285 | triggerRate += xSection[mult] * triggerEff[mult]; | |
1286 | ||
1287 | Printf("Raw value for 1 collision is %e", triggerRate); | |
1288 | ||
1289 | Double_t triggerRate2 = 0; | |
1290 | for (Int_t mult = 0; mult < max; mult++) | |
1291 | for (Int_t mult2 = mult; mult2 < max; mult2++) | |
1292 | if (mult+mult2 < max) | |
1293 | triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2]; | |
1294 | ||
1295 | Printf("Raw value for 2 collisions is %e", triggerRate2); | |
1296 | ||
69b09e3b | 1297 | for (Double_t doubleRate = 0; doubleRate <= 0.3; doubleRate += 0.005) |
144ff489 | 1298 | { |
1299 | Float_t totalContamination = (triggerRate2 * doubleRate) / (triggerRate + triggerRate2 * doubleRate); | |
1300 | ||
1301 | //Printf("Total contamination is %.1f%%", totalContamination * 100); | |
1302 | ||
1303 | graph->SetPoint(graph->GetN(), doubleRate, totalContamination); | |
1304 | } | |
1305 | ||
1306 | graph->Draw((currentCut == 0) ? "A*" : "* SAME"); | |
1307 | graph->GetXaxis()->SetRangeUser(0, 1); | |
1308 | } | |
1309 | } | |
1310 | } | |
1311 | ||
69b09e3b | 1312 | void AliHighMultiplicitySelector::Contamination3() |
1313 | { | |
1314 | // | |
b3b260d1 | 1315 | // draws the contamination as function of treshold depending on a number a set of input MC and rate parameters |
69b09e3b | 1316 | // |
1317 | ||
1318 | /* | |
1319 | ||
1320 | gSystem->Load("libANALYSIS"); | |
1321 | gSystem->Load("libPWG0base"); | |
1322 | .L AliHighMultiplicitySelector.cxx+g | |
1323 | x = new AliHighMultiplicitySelector(); | |
1324 | x->ReadHistograms("highmult_hijing100k.root"); | |
1325 | x->Contamination3(); | |
1326 | ||
1327 | */ | |
1328 | ||
b3b260d1 | 1329 | // output file |
1330 | TFile* output = TFile::Open("contamination3.root", "RECREATE"); | |
1331 | ||
69b09e3b | 1332 | // get x-sections |
b3b260d1 | 1333 | TFile* file = TFile::Open("crosssectionEx_10TeV.root"); |
69b09e3b | 1334 | if (!file) |
1335 | return; | |
b3b260d1 | 1336 | |
1337 | TCanvas* c = new TCanvas; | |
1338 | c->SetGridx(); | |
1339 | c->SetGridy(); | |
1340 | ||
1341 | TLegend* legend = new TLegend(0.7, 0.2, 1, 0.5); | |
1342 | legend->SetNColumns(2); | |
1343 | ||
1344 | TH2* dummy = new TH2F("dummy", ";Layer 1 Threshold;Contamination", 100, 95, 255, 100, 0, 1); | |
1345 | dummy->SetStats(kFALSE); | |
1346 | dummy->Draw(); | |
1347 | ||
1348 | for (Int_t mc = 0; mc < 6; mc++) | |
69b09e3b | 1349 | { |
b3b260d1 | 1350 | TH1* xSections[2]; |
1351 | TString str; | |
1352 | str.Form("xSection2Ex_%d_%d", mc/3, mc%3); | |
1353 | Printf("%s", str.Data()); | |
1354 | file->cd(); | |
1355 | xSections[0] = dynamic_cast<TH1*> (gFile->Get(str)); | |
1356 | //xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1357 | ||
1358 | // prob for a collision in a bunch crossing | |
1359 | Int_t nRates = 1; | |
1360 | //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2}; | |
1361 | Float_t rates[] = {0.0636}; | |
1362 | ||
1363 | // bunch crossing rate in Hz | |
1364 | Float_t bunchCrossingRate = 24. * 11245.5; | |
1365 | ||
1366 | Int_t colors[] = { 2, 3, 4, 6, 7, 8 }; | |
1367 | Int_t markers[] = { 7, 2, 4, 5, 6, 27 }; | |
1368 | ||
1369 | // put to 2 for second layer | |
1370 | for (Int_t i=0; i<1; ++i) | |
69b09e3b | 1371 | { |
b3b260d1 | 1372 | if (!xSections[i]) |
1373 | continue; | |
1374 | ||
1375 | // relative x-section, once we have a collision | |
1376 | xSections[i]->Scale(1.0 / xSections[i]->Integral()); | |
1377 | ||
1378 | Int_t max = xSections[i]->GetNbinsX(); | |
1379 | max = 500; | |
1380 | ||
1381 | Float_t* xSection = new Float_t[max]; | |
1382 | for (Int_t mult = 0; mult < max; mult++) | |
1383 | xSection[mult] = xSections[i]->GetBinContent(mult+1); | |
1384 | ||
1385 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
1386 | ||
1387 | for (Int_t currentRate = 0; currentRate<nRates; ++currentRate) | |
69b09e3b | 1388 | { |
b3b260d1 | 1389 | TGraph* graph = new TGraph; |
1390 | graph->SetMarkerColor(colors[currentRate]); | |
1391 | graph->SetMarkerStyle(markers[currentRate]); | |
1392 | ||
1393 | TGraph* graph2 = new TGraph; | |
1394 | ||
1395 | Float_t rate = rates[currentRate]; | |
1396 | ||
1397 | Double_t singleRate = TMath::Poisson(1, rate); | |
1398 | Double_t doubleRate = TMath::Poisson(2, rate); | |
1399 | Double_t tripleRate = TMath::Poisson(3, rate); | |
1400 | ||
1401 | Printf("single = %f, double = %f, triple = %f", singleRate, doubleRate, tripleRate); | |
1402 | ||
1403 | for (Int_t cut = 100; cut <= 251; cut += 10) | |
1404 | { | |
1405 | Printf("Cut at %d", cut); | |
1406 | ||
1407 | TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); | |
1408 | Float_t* triggerEff = new Float_t[max]; | |
1409 | for (Int_t mult = 0; mult < max; mult++) | |
1410 | triggerEff[mult] = triggerEffHist->GetBinContent(mult+1); | |
1411 | ||
1412 | Double_t triggerRate = 0; | |
1413 | for (Int_t mult = 0; mult < max; mult++) | |
1414 | triggerRate += xSection[mult] * triggerEff[mult]; | |
1415 | ||
1416 | //Printf(" Raw value for 1 collision is %e; Rate: %.1f Hz", triggerRate, triggerRate * singleRate * bunchCrossingRate); | |
1417 | ||
1418 | Double_t triggerRate2 = 0; | |
1419 | for (Int_t mult = 0; mult < max; mult++) | |
1420 | for (Int_t mult2 = mult; mult2 < max; mult2++) | |
1421 | if (mult+mult2 < max) | |
1422 | triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2]; | |
1423 | ||
1424 | //Printf(" Raw value for 2 collisions is %e; Rate: %.1f Hz", triggerRate2, triggerRate2 * doubleRate * bunchCrossingRate); | |
1425 | ||
1426 | Double_t triggerRate3 = 0; | |
1427 | for (Int_t mult = 0; mult < max; mult++) | |
1428 | for (Int_t mult2 = 0; mult2 < max; mult2++) | |
1429 | for (Int_t mult3 = 0; mult3 < max; mult3++) | |
1430 | if (mult+mult2+mult3 < max) | |
1431 | triggerRate3 += xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3]; | |
1432 | ||
1433 | //Printf(" Raw value for 3 collisions is %e; Rate: %.1f Hz", triggerRate3, triggerRate3 * tripleRate * bunchCrossingRate); | |
1434 | ||
1435 | Printf(" Rates: %.1f Hz; %.1f Hz; %.1f Hz", triggerRate * singleRate * bunchCrossingRate, triggerRate2 * doubleRate * bunchCrossingRate, triggerRate3 * tripleRate * bunchCrossingRate); | |
1436 | ||
1437 | Float_t totalTrigger = (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate); | |
1438 | ||
1439 | Printf(" Total trigger rate: %.1f Hz", totalTrigger * bunchCrossingRate); | |
1440 | ||
1441 | //if (totalTrigger * bunchCrossingRate > 200) | |
1442 | // continue; | |
1443 | ||
1444 | Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / totalTrigger; | |
1445 | //if (totalContamination > 0.99) | |
1446 | // break; | |
1447 | ||
1448 | Printf(" Total contamination is %.1f%%", totalContamination * 100); | |
1449 | ||
1450 | graph->SetPoint(graph->GetN(), cut, totalContamination); | |
1451 | graph2->SetPoint(graph->GetN(), cut, totalTrigger * bunchCrossingRate); | |
1452 | } | |
1453 | ||
1454 | graph->SetMarkerStyle(mc+20); | |
1455 | graph->SetMarkerColor(currentRate+1); | |
1456 | graph->Draw("P SAME"); | |
1457 | graph->GetXaxis()->SetTitle("Layer 1 threshold"); | |
1458 | graph->GetYaxis()->SetTitle("Contamination"); | |
1459 | graph->GetYaxis()->SetRangeUser(0, 1); | |
1460 | ||
1461 | if (currentRate == 0) | |
1462 | { | |
1463 | const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" }; | |
1464 | legend->AddEntry(graph, legendLabel[mc], "P"); | |
1465 | } | |
69b09e3b | 1466 | |
b3b260d1 | 1467 | output->cd(); |
1468 | graph->Write(Form("%s_%d_cont", str.Data(), currentRate)); | |
1469 | graph2->Write(Form("%s_%d_rate", str.Data(), currentRate)); | |
1470 | } | |
1471 | } | |
1472 | } | |
1473 | ||
1474 | output->Close(); | |
1475 | ||
1476 | legend->Draw(); | |
1477 | } | |
69b09e3b | 1478 | |
b3b260d1 | 1479 | void AliHighMultiplicitySelector::Contamination_Reach() |
1480 | { | |
1481 | // plot the multiplicity reach based on the output from Contamination3() | |
1482 | // for each rate case and each MC, a certain number of events is required starting counting from the highest multiplicity | |
1483 | // note that the reach of different MC cannot be compared with each other | |
69b09e3b | 1484 | |
b3b260d1 | 1485 | /* |
69b09e3b | 1486 | |
b3b260d1 | 1487 | gSystem->Load("libANALYSIS"); |
1488 | gSystem->Load("libPWG0base"); | |
1489 | .L AliHighMultiplicitySelector.cxx+g | |
1490 | x = new AliHighMultiplicitySelector(); | |
1491 | x->ReadHistograms("highmult_hijing100k.root"); | |
1492 | x->Contamination_Reach(); | |
69b09e3b | 1493 | |
b3b260d1 | 1494 | */ |
69b09e3b | 1495 | |
b3b260d1 | 1496 | TCanvas* c = new TCanvas("c", "c", 800, 600); |
1497 | c->Divide(2, 3); | |
1498 | ||
1499 | // prob for a collision in a bunch crossing | |
1500 | Int_t nRates = 1; | |
1501 | //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2}; | |
1502 | Float_t rates[] = {0.0636}; | |
69b09e3b | 1503 | |
b3b260d1 | 1504 | // bunch crossing rate in Hz |
1505 | Float_t bunchCrossingRate = 24. * 11245.5; | |
69b09e3b | 1506 | |
b3b260d1 | 1507 | TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;multiplicity reach", 100, 0, 0.3, 100, 50, 350); |
1508 | //TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;fractional cross-section", 100, 0, 0.3, 1000, 1e-6, 0.1); | |
1509 | dummy->SetStats(kFALSE); | |
69b09e3b | 1510 | |
b3b260d1 | 1511 | const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" }; |
69b09e3b | 1512 | |
b3b260d1 | 1513 | TFile* mcFile = TFile::Open("crosssectionEx_10TeV.root"); |
1514 | TFile* contFile = TFile::Open("contamination3.root"); | |
69b09e3b | 1515 | |
b3b260d1 | 1516 | // for comparison: how many MB events can one take at the same time |
1517 | Int_t mbEvents = 2e6 * 500; | |
69b09e3b | 1518 | |
b3b260d1 | 1519 | for (Int_t mc=0; mc<6; mc++) |
1520 | { | |
1521 | mcFile->cd(); | |
1522 | TH1* mcHist = (TH1*) gFile->Get(Form("xSection2Ex_%d_%d", mc/3, mc%3)); | |
1523 | mcHist->Scale(1.0 / mcHist->Integral()); | |
1524 | ||
1525 | c->cd(mc+1);//->SetLogy(); | |
1526 | c->SetGridx(); | |
1527 | c->SetGridy(); | |
1528 | dummy->Draw(); | |
1529 | ||
1530 | Int_t color = 0; | |
1531 | for (Int_t requiredEvents = 300; requiredEvents <= 3000000; requiredEvents *= 10) | |
1532 | { | |
1533 | TGraph* reach = new TGraph; | |
1534 | ||
1535 | color++; | |
1536 | if (color == 5) | |
1537 | color++; | |
1538 | ||
1539 | Float_t requiredRate = (Float_t) requiredEvents / 1e6; | |
1540 | Printf("Required rate is %f", requiredRate); | |
1541 | ||
1542 | // find reach without trigger | |
1543 | Int_t mbReach = 1000; | |
1544 | while (mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) < (Float_t) requiredEvents / mbEvents && mbReach > 1) | |
1545 | mbReach--; | |
1546 | Printf("MB reach is %d with %f events", mbReach, mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) * mbEvents); | |
1547 | ||
1548 | for (Int_t rate=0; rate<nRates; rate++) | |
1549 | { | |
1550 | contFile->cd(); | |
1551 | TGraph* cont = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_cont", mc/3, mc%3, rate)); | |
1552 | TGraph* rateh = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_rate", mc/3, mc%3, rate)); | |
1553 | ||
1554 | Double_t singleRate = TMath::Poisson(1, rates[rate]); | |
1555 | Double_t totalCollRate = singleRate * bunchCrossingRate; | |
1556 | Printf("collisions/bc: %f; coll. rate: %f", singleRate, totalCollRate); | |
1557 | ||
1558 | // find 200 Hz limit | |
1559 | Int_t low = 100; | |
1560 | while (rateh->Eval(low) > 200) | |
1561 | low++; | |
1562 | ||
1563 | // find contamination limit | |
1564 | Int_t high = 100; | |
1565 | while (cont->Eval(high) < 0.9 && high < 250) | |
1566 | high++; | |
1567 | ||
1568 | Printf("MC %d Rate %d: Acceptable threshold region is %d to %d", mc, rate, low, high); | |
1569 | // find reachable multiplicity; include contamination in rate calculation | |
1570 | // TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); | |
1571 | ||
1572 | // trigger efficiency as function of multiplicity in range mult <= ... <= high | |
1573 | //new TCanvas; fMvsL1->Draw("COLZ"); | |
1574 | TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL1, low, high); | |
1575 | ||
1576 | //new TCanvas; triggerEffHist->DrawCopy(); | |
1577 | triggerEffHist->Multiply(mcHist); | |
1578 | ||
1579 | Float_t fractionXSection = triggerEffHist->Integral(); | |
1580 | Printf("The fraction of the cross-section is %f", fractionXSection); | |
1581 | ||
1582 | //new TCanvas; triggerEffHist->DrawCopy(); | |
1583 | triggerEffHist->Scale(totalCollRate); | |
1584 | //new TCanvas; triggerEffHist->DrawCopy(); gPad->SetLogy(); | |
1585 | ||
1586 | Float_t achievedRate = 0; | |
1587 | Int_t mult = 1000; | |
1588 | while (1) | |
1589 | { | |
1590 | achievedRate = triggerEffHist->Integral(triggerEffHist->FindBin(mult), triggerEffHist->GetNbinsX()); | |
1591 | ||
1592 | if (achievedRate >= requiredRate) | |
1593 | break; | |
1594 | ||
1595 | if (mult == 1) | |
1596 | break; | |
1597 | ||
1598 | mult--; | |
1599 | } | |
1600 | ||
1601 | Printf("Achieved rate %f above multiplicity %d", achievedRate, mult); | |
1602 | ||
1603 | if (achievedRate < requiredRate) | |
1604 | { | |
1605 | Printf("Achieved rate too low"); | |
1606 | continue; | |
1607 | } | |
1608 | ||
1609 | reach->SetPoint(reach->GetN(), rates[rate], mult); | |
1610 | //reach->SetPoint(reach->GetN(), rates[rate], fractionXSection); | |
1611 | ||
1612 | ||
1613 | //return; | |
69b09e3b | 1614 | } |
b3b260d1 | 1615 | |
1616 | reach->SetMarkerColor(color); | |
1617 | reach->Draw("SAME*"); | |
1618 | ||
1619 | TLine* line = new TLine(0, mbReach, 0.3, mbReach); | |
1620 | line->SetLineColor(color); | |
1621 | //line->Draw(); | |
1622 | } | |
1623 | ||
1624 | TText* text = new TText; | |
1625 | text->DrawText(0.2, 325, legendLabel[mc]); | |
1626 | //return; | |
69b09e3b | 1627 | } |
1628 | } | |
1629 | ||
9608df41 | 1630 | void AliHighMultiplicitySelector::DrawHistograms() |
1631 | { | |
a9aed06c | 1632 | // draws the histograms |
1633 | ||
9608df41 | 1634 | /* |
1635 | ||
1636 | gSystem->Load("libPWG0base"); | |
1637 | .L AliHighMultiplicitySelector.cxx+ | |
1638 | x = new AliHighMultiplicitySelector(); | |
1639 | x->ReadHistograms("highmult_pythia.root"); | |
1640 | x->DrawHistograms(); | |
1641 | ||
1642 | ||
1643 | gSystem->Load("libPWG0base"); | |
1644 | .L AliHighMultiplicitySelector.cxx+ | |
1645 | x = new AliHighMultiplicitySelector(); | |
1646 | x->ReadHistograms("highmult_hijing.root"); | |
1647 | x->DrawHistograms(); | |
1648 | ||
1649 | gSystem->Load("libPWG0base"); | |
1650 | .L AliHighMultiplicitySelector.cxx+ | |
1651 | x = new AliHighMultiplicitySelector(); | |
1652 | x->ReadHistograms("highmult_central.root"); | |
1653 | x->DrawHistograms(); | |
1654 | ||
745d6088 | 1655 | gSystem->Load("libANALYSIS"); |
9608df41 | 1656 | gSystem->Load("libPWG0base"); |
1657 | .L AliHighMultiplicitySelector.cxx+ | |
1658 | x = new AliHighMultiplicitySelector(); | |
1659 | x->ReadHistograms("highmult_hijing100k.root"); | |
1660 | x->DrawHistograms(); | |
1661 | ||
1662 | */ | |
1663 | ||
1664 | /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400); | |
1665 | ||
1666 | fChipsLayer2->SetLineColor(2); | |
1667 | fChipsLayer2->SetStats(kFALSE); | |
1668 | fChipsLayer1->SetStats(kFALSE); | |
1669 | fChipsLayer2->SetTitle(""); | |
1670 | fChipsLayer2->DrawCopy(); | |
1671 | fChipsLayer1->DrawCopy("SAME"); | |
1672 | canvas->SaveAs("chips.gif"); | |
1673 | ||
1674 | canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400); | |
1675 | fL1vsL2->SetStats(kFALSE); | |
1676 | fL1vsL2->DrawCopy("COLZ"); | |
1677 | gPad->SetLogz(); | |
1678 | canvas->SaveAs("L1vsL2.gif");*/ | |
1679 | ||
f6093269 | 1680 | TCanvas *canvas = new TCanvas("L1", "L1", 800, 600); |
1681 | canvas->SetTopMargin(0.05); | |
1682 | canvas->SetRightMargin(0.12); | |
9608df41 | 1683 | fMvsL1->SetStats(kFALSE); |
1684 | fMvsL1->DrawCopy("COLZ"); | |
1685 | gPad->SetLogz(); | |
1686 | ||
1687 | canvas->SaveAs("L1NoCurve.gif"); | |
f6093269 | 1688 | canvas->SaveAs("L1NoCurve.eps"); |
9608df41 | 1689 | |
745d6088 | 1690 | TLine* line = new TLine(fMvsL1->GetXaxis()->GetXmin(), 150, fMvsL1->GetXaxis()->GetXmax(), 150); |
1691 | line->SetLineWidth(2); | |
1692 | line->SetLineColor(kRed); | |
1693 | line->Draw(); | |
1694 | ||
1695 | canvas->SaveAs("L1NoCurveCut.gif"); | |
1696 | canvas->SaveAs("L1NoCurveCut.eps"); | |
1697 | ||
1698 | return; | |
1699 | ||
9608df41 | 1700 | // draw corresponding theoretical curve |
1701 | TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000); | |
1702 | func->SetParameter(0, 400-5*2); | |
1703 | func->DrawCopy("SAME"); | |
1704 | ||
1705 | canvas->SaveAs("L1.gif"); | |
1706 | ||
1707 | canvas = new TCanvas("L2", "L2", 600, 400); | |
1708 | //fMvsL2->GetYaxis()->SetRangeUser(0, 150); | |
1709 | fMvsL2->SetStats(kFALSE); | |
1710 | fMvsL2->DrawCopy("COLZ"); | |
1711 | gPad->SetLogz(); | |
1712 | func->SetParameter(0, 800-5*4); | |
1713 | func->DrawCopy("SAME"); | |
1714 | canvas->SaveAs("L2.gif"); | |
1715 | ||
f6093269 | 1716 | // get x-sections |
1717 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1718 | if (file) | |
1719 | { | |
1720 | TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1721 | TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1722 | ||
1723 | MakeGraphs2("Layer1", xSection2, fMvsL1); | |
1724 | return; | |
1725 | ||
1726 | // 5 times the mean | |
1727 | //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2); | |
1728 | //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2)); | |
1729 | ||
1730 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8)); | |
1731 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8)); | |
1732 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9)); | |
1733 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9)); | |
1734 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10)); | |
1735 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10)); | |
1736 | ||
1737 | file->Close(); | |
1738 | } | |
1739 | ||
9608df41 | 1740 | // make spread hists |
1741 | TGraph* spread = new TGraph; | |
1742 | spread->SetTitle("Spread L1;true multiplicity;RMS"); | |
1743 | ||
1744 | for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i) | |
1745 | { | |
1746 | TH1* proj = fMvsL1->ProjectionY("proj", i, i); | |
1747 | spread->SetPoint(spread->GetN(), i, proj->GetRMS()); | |
1748 | } | |
1749 | ||
1750 | canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400); | |
1751 | spread->Draw("A*"); | |
1752 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1753 | ||
1754 | TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150); | |
1755 | log->SetParLimits(0, 0, 10); | |
1756 | log->SetParLimits(1, 1e-5, 10); | |
1757 | ||
1758 | spread->Fit(log, "", "", 1, 150); | |
1759 | log->DrawCopy("SAME"); | |
1760 | ||
1761 | TGraph* spread2 = new TGraph; | |
1762 | spread2->SetTitle("Spread L2;true multiplicity;RMS"); | |
1763 | ||
1764 | for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i) | |
1765 | { | |
1766 | TH1* proj = fMvsL2->ProjectionY("proj", i, i); | |
1767 | spread2->SetPoint(spread2->GetN(), i, proj->GetRMS()); | |
1768 | } | |
1769 | ||
1770 | canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400); | |
1771 | spread2->Draw("A*"); | |
1772 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1773 | ||
1774 | spread2->Fit(log, "", "", 1, 150); | |
1775 | log->DrawCopy("SAME"); | |
1776 | ||
9608df41 | 1777 | canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400); |
1778 | fClvsL1->SetStats(kFALSE); | |
1779 | fClvsL1->DrawCopy("COLZ"); | |
1780 | gPad->SetLogz(); | |
1781 | ||
1782 | func->SetParameter(0, 400-5*2); | |
1783 | func->DrawCopy("SAME"); | |
1784 | ||
1785 | canvas->SaveAs("Clusters_L1.gif"); | |
1786 | ||
1787 | canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400); | |
1788 | fClvsL2->SetStats(kFALSE); | |
1789 | fClvsL2->DrawCopy("COLZ"); | |
1790 | gPad->SetLogz(); | |
1791 | func->SetParameter(0, 800-5*4); | |
1792 | func->DrawCopy("SAME"); | |
1793 | canvas->SaveAs("Clusters_L2.gif"); | |
1794 | ||
1795 | canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400); | |
1796 | //fChipsFired->GetYaxis()->SetRangeUser(0, 150); | |
1797 | fChipsFired->SetStats(kFALSE); | |
1798 | fChipsFired->DrawCopy("COLZ"); | |
1799 | canvas->SaveAs("ChipsFired.gif"); | |
1800 | ||
1801 | /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1); | |
1802 | TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2); | |
1803 | ||
1804 | for (Int_t treshold = 0; treshold < 800; treshold++) | |
1805 | { | |
1806 | if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0) | |
1807 | { | |
1808 | TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram()); | |
1809 | if (mult) | |
1810 | tresholdHistL1->Fill(treshold, mult->GetMean()); | |
1811 | } | |
1812 | if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0) | |
1813 | { | |
1814 | TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram()); | |
1815 | if (mult) | |
1816 | tresholdHistL2->Fill(treshold, mult->GetMean()); | |
1817 | } | |
1818 | } | |
1819 | ||
1820 | canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400); | |
1821 | tresholdHistL1->Draw(); | |
1822 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1823 | ||
1824 | canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400); | |
1825 | tresholdHistL2->Draw(); | |
1826 | canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/ | |
1827 | ||
1828 | fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff"); | |
1829 | fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff"); | |
1830 | ||
1831 | canvas = new TCanvas("Efficiency", "Efficiency", 600, 400); | |
1832 | fPrimaryL1->GetHistogram()->SetStats(kFALSE); | |
1833 | fPrimaryL1->GetHistogram()->Draw(); | |
1834 | fPrimaryL2->GetHistogram()->SetLineColor(2); | |
1835 | fPrimaryL2->GetHistogram()->SetStats(kFALSE); | |
1836 | fPrimaryL2->GetHistogram()->Draw("SAME"); | |
1837 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1838 | ||
1839 | canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400); | |
1840 | fClusterZL1->Rebin(2); | |
1841 | fClusterZL1->Draw(); | |
1842 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1843 | ||
1844 | canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400); | |
1845 | fClusterZL2->Draw(); | |
1846 | fClusterZL2->Rebin(2); | |
1847 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1848 | } | |
a9aed06c | 1849 | |
1850 | TGraph* AliHighMultiplicitySelector::IntFractRate() | |
1851 | { | |
1852 | // A plot which shows the fractional rate above any threshold | |
1853 | // as function of threshold (i.e. the integral of dSigma/dN as function of | |
1854 | // N, normalised to 1 for N=0) | |
1855 | ||
1856 | /* | |
d1594f8a | 1857 | gSystem->Load("libANALYSIS"); |
a9aed06c | 1858 | gSystem->Load("libPWG0base"); |
1859 | .L AliHighMultiplicitySelector.cxx+ | |
1860 | x = new AliHighMultiplicitySelector(); | |
1861 | x->ReadHistograms("highmult_hijing100k.root"); | |
1862 | x->IntFractRate(); | |
1863 | */ | |
1864 | ||
1865 | // get x-sections | |
1866 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1867 | if (!file) | |
1868 | return 0; | |
1869 | ||
1870 | TH1* xSection; | |
1871 | xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1872 | ||
1873 | TGraph* result = new TGraph; | |
1874 | ||
1875 | for (Int_t threshold = 0; threshold < 300; threshold += 2) | |
1876 | { | |
1877 | TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold); | |
1878 | ||
1879 | //new TCanvas; proj->DrawCopy(); | |
1880 | ||
1881 | Double_t integral = proj->Integral(); | |
1882 | ||
d1594f8a | 1883 | Printf("Cut at %d, integral is %e", threshold, integral); |
a9aed06c | 1884 | |
1885 | result->SetPoint(result->GetN(), threshold, integral); | |
1886 | } | |
1887 | ||
1888 | TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400); | |
1889 | gPad->SetLogy(); | |
1890 | ||
1891 | result->Draw("A*"); | |
1892 | result->GetXaxis()->SetTitle("threshold"); | |
1893 | result->GetYaxis()->SetTitle("integrated fractional rate above threshold"); | |
1894 | ||
1895 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1896 | ||
1897 | return result; | |
1898 | } | |
d1594f8a | 1899 | |
1900 | void AliHighMultiplicitySelector::MBComparison() | |
1901 | { | |
1902 | // | |
1903 | // finds the threshold from which onwards the number of found events above N times the mean | |
1904 | // is higher using a high mult. trigger than just triggering with MB | |
1905 | // | |
1906 | ||
1907 | /* | |
1908 | ||
1909 | gSystem->Load("libANALYSIS"); | |
1910 | gSystem->Load("libPWG0base"); | |
1911 | .L AliHighMultiplicitySelector.cxx+g | |
1912 | x = new AliHighMultiplicitySelector(); | |
1913 | x->ReadHistograms("highmult_hijing100k.root"); | |
1914 | x->MBComparison(); | |
1915 | ||
1916 | */ | |
1917 | ||
1918 | // get x-sections | |
1919 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1920 | if (!file) | |
1921 | return; | |
1922 | ||
1923 | TH1* xSections[2]; | |
1924 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1925 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1926 | ||
1927 | // rate = L * sigma | |
1928 | // sigma ~ 80 mb (Pythia 14 TeV) | |
1929 | // 10^28 lum --> 8e2 Hz | |
1930 | // 10^31 lum --> 8e5 Hz | |
1931 | Int_t nRates = 4; | |
1932 | Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 }; | |
1933 | ||
1934 | // threshold in number of fired chips for corresponding rate | |
1935 | //Int_t cuts[] = { 104, 134, 154, 170 }; // values for 20 Hz | |
1936 | Int_t cuts[] = { 82, 124, 147, 164 }; // values for 50 Hz | |
1937 | ||
1938 | // bandwidth, fractions (for MB, high mult.) | |
1939 | Float_t bandwidth = 1e3; | |
1940 | Float_t fractionMB = 0.5; | |
1941 | Float_t fractionHM = 0.05; | |
1942 | ||
1943 | // different limits to define "interesting events" | |
1944 | Int_t nLimits = 9; | |
1945 | Int_t limits[] = { 0, 1, 2, 4, 6, 7, 8, 9, 10 }; | |
1946 | ||
1947 | // put to 2 for second layer | |
1948 | for (Int_t i=0; i<1; ++i) | |
1949 | { | |
1950 | if (!xSections[i]) | |
1951 | continue; | |
1952 | ||
1953 | TH1* xSection = xSections[i]; | |
1954 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
1955 | ||
1956 | // relative x-section, once we have a collision | |
1957 | xSection->Scale(1.0 / xSection->Integral()); | |
1958 | ||
1959 | xSection->SetStats(kFALSE); | |
1960 | xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2"); | |
1961 | xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5)); | |
1962 | xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350); | |
1963 | xSection->GetYaxis()->SetTitleOffset(1.2); | |
1964 | ||
1965 | TCanvas* canvas = new TCanvas("MBComparison", "MBComparison", 1000, 800); | |
1966 | canvas->Divide(3, 3); | |
1967 | ||
1968 | for (Int_t currentLimit = 0; currentLimit<nLimits; currentLimit++) | |
1969 | { | |
1970 | // limit is N times the mean | |
1971 | Int_t limit = (Int_t) (xSection->GetMean() * limits[currentLimit]); | |
1972 | if (limit < 1) | |
1973 | limit = 1; | |
1974 | ||
1975 | TGraph* graphMB = new TGraph; | |
1976 | graphMB->SetTitle(Form("Events with %d times above <n> (i.e. n >= %d)", limits[currentLimit], limit)); | |
1977 | graphMB->SetMarkerStyle(20); | |
1978 | ||
1979 | TGraph* graphBoth = new TGraph; | |
1980 | graphBoth->SetMarkerStyle(21); | |
1981 | ||
1982 | Float_t min = bandwidth; | |
1983 | Float_t max = 0; | |
1984 | ||
1985 | for (Int_t current = 0; current<nRates; ++current) | |
1986 | { | |
1987 | Float_t rate = rates[current]; | |
1988 | Int_t cut = cuts[current]; | |
1989 | ||
1990 | TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff"); | |
1991 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
1992 | ||
1993 | Float_t downScaleMB1 = rate / bandwidth; | |
1994 | if (downScaleMB1 < 1) | |
1995 | downScaleMB1 = 1; | |
1996 | ||
1997 | Float_t downScaleMB2 = rate / (bandwidth * fractionMB); | |
1998 | if (downScaleMB2 < 1) | |
1999 | downScaleMB2 = 1; | |
2000 | ||
2001 | Float_t downScaleHM = rate * proj->Integral(1, xSection->GetNbinsX()) / (bandwidth * fractionHM); | |
2002 | if (downScaleHM < 1) | |
2003 | downScaleHM = 1; | |
2004 | ||
2005 | Float_t rateMB1 = rate / downScaleMB1 * xSection->Integral(limit, xSection->GetNbinsX()); | |
2006 | Float_t rateMB2 = rate / downScaleMB2 * xSection->Integral(limit, xSection->GetNbinsX()); | |
2007 | Float_t rateHM = rate / downScaleHM * proj->Integral(limit, xSection->GetNbinsX()); | |
2008 | Float_t combinedRate = rateMB2 + rateHM; | |
2009 | ||
2010 | graphMB->SetPoint(graphMB->GetN(), rate, rateMB1); | |
2011 | graphBoth->SetPoint(graphBoth->GetN(), rate, combinedRate); | |
2012 | ||
2013 | min = TMath::Min(min, TMath::Min(rateMB1, combinedRate)); | |
2014 | max = TMath::Max(min, TMath::Max(rateMB1, combinedRate)); | |
2015 | ||
2016 | Printf("The rates for events with %d times above <n> (i.e. n >= %d) at a rate of %.2e Hz is:", limits[currentLimit], limit, rate); | |
2017 | Printf(" %.2e Hz in MB-only mode", rateMB1); | |
2018 | Printf(" %.2e Hz = %.2e Hz + %.2e Hz in MB + high mult. mode", combinedRate, rateMB2, rateHM); | |
2019 | ||
2020 | Printf(" The downscale factors are: %.2f %.2f %.2f", downScaleMB1, downScaleMB2, downScaleHM); | |
2021 | ||
2022 | Int_t triggerLimit = 0; | |
2023 | for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++) | |
2024 | if (triggerEff->GetBinContent(bin) < 0.5) | |
2025 | triggerLimit = (Int_t) triggerEff->GetXaxis()->GetBinCenter(bin); | |
2026 | ||
2027 | Printf(" Efficiency limit (50%%) is at multiplicity %d", triggerLimit); | |
144ff489 | 2028 | Float_t fractionGood = proj->Integral(triggerLimit, proj->GetNbinsX()) / proj->Integral(); |
2029 | Printf(" %.2f %% of the events are above the trigger limit", fractionGood * 100); | |
d1594f8a | 2030 | |
2031 | if (triggerLimit > limit) | |
2032 | Printf(" WARNING: interesting events also counted inside the trigger limit"); | |
2033 | ||
2034 | Printf(""); | |
2035 | } | |
2036 | ||
2037 | canvas->cd(currentLimit+1)->SetLogx(); | |
2038 | canvas->cd(currentLimit+1)->SetLogy(); | |
2039 | ||
2040 | graphMB->Draw("AP"); | |
2041 | graphBoth->Draw("P SAME"); | |
2042 | ||
2043 | graphMB->GetYaxis()->SetRangeUser(0.5 * min, 2 * max); | |
2044 | graphMB->GetXaxis()->SetTitle("Raw rate in Hz"); | |
2045 | graphMB->GetYaxis()->SetTitle("Event rate in Hz"); | |
2046 | } | |
2047 | } | |
2048 | } |