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