]>
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> | |
19 | ||
20 | #include <AliLog.h> | |
21 | #include <AliESD.h> | |
22 | #include <AliRunLoader.h> | |
23 | #include <AliStack.h> | |
24 | ||
25 | #include <AliITSgeom.h> | |
26 | #include <AliITSLoader.h> | |
27 | #include <AliITSdigitSPD.h> | |
28 | #include <AliITSRecPoint.h> | |
29 | ||
30 | #include "AliPWG0Helper.h" | |
31 | ||
32 | // | |
33 | // | |
34 | ||
35 | ClassImp(AliHighMultiplicitySelector) | |
36 | ||
37 | AliHighMultiplicitySelector::AliHighMultiplicitySelector() : | |
38 | AliSelectorRL(), | |
39 | fChipsLayer1(0), | |
40 | fChipsLayer2(0), | |
41 | fL1vsL2(0), | |
42 | fMvsL1(0), | |
43 | fMvsL2(0), | |
44 | fChipsFired(0), | |
45 | fPrimaryL1(0), | |
46 | fPrimaryL2(0), | |
47 | fClusterZL1(0), | |
48 | fClusterZL2(0), | |
49 | fClvsL1(0), | |
50 | fClvsL2(0), | |
f91bc12d | 51 | centralRegion(kFALSE) |
9608df41 | 52 | { |
53 | // | |
54 | // Constructor. Initialization of pointers | |
55 | // | |
56 | } | |
57 | ||
58 | AliHighMultiplicitySelector::~AliHighMultiplicitySelector() | |
59 | { | |
60 | // | |
61 | // Destructor | |
62 | // | |
63 | ||
64 | // histograms are in the output list and deleted when the output | |
65 | // list is deleted by the TSelector dtor | |
66 | } | |
67 | ||
68 | void AliHighMultiplicitySelector::SlaveBegin(TTree *tree) | |
69 | { | |
70 | AliSelectorRL::SlaveBegin(tree); | |
71 | ||
72 | fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5); | |
73 | ||
74 | fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters"); | |
75 | fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters"); | |
76 | ||
77 | fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20); | |
78 | fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20); | |
79 | } | |
80 | ||
81 | void AliHighMultiplicitySelector::Init(TTree* tree) | |
82 | { | |
83 | // read the user objects | |
84 | ||
85 | AliSelectorRL::Init(tree); | |
86 | ||
87 | // enable only the needed branches | |
88 | if (tree) | |
89 | { | |
90 | tree->SetBranchStatus("*", 0); | |
91 | tree->SetBranchStatus("fTriggerMask", 1); | |
92 | ||
93 | /*if (fTree->GetCurrentFile()) | |
94 | { | |
95 | TString fileName(fTree->GetCurrentFile()->GetName()); | |
96 | fileName.ReplaceAll("AliESDs", "geometry"); | |
97 | ||
98 | // load geometry | |
99 | TGeoManager::Import(fileName); | |
100 | }*/ | |
101 | } | |
102 | } | |
103 | ||
104 | Bool_t AliHighMultiplicitySelector::Process(Long64_t entry) | |
105 | { | |
106 | // | |
107 | // processing | |
108 | // | |
109 | ||
110 | if (AliSelectorRL::Process(entry) == kFALSE) | |
111 | return kFALSE; | |
112 | ||
113 | // Check prerequisites | |
114 | if (!fESD) | |
115 | { | |
116 | AliDebug(AliLog::kError, "ESD branch not available"); | |
117 | return kFALSE; | |
118 | } | |
119 | ||
120 | Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1); | |
121 | ||
122 | if (!eventTriggered) | |
123 | { | |
124 | AliDebug(AliLog::kDebug, "Event not triggered. Skipping."); | |
125 | return kTRUE; | |
126 | } | |
127 | ||
128 | AliStack* stack = GetStack(); | |
129 | if (!stack) | |
130 | { | |
131 | AliDebug(AliLog::kError, "Stack not available"); | |
132 | return kFALSE; | |
133 | } | |
134 | ||
135 | Int_t nPrim = stack->GetNprimary(); | |
136 | Int_t multiplicity21 = 0; | |
137 | Int_t multiplicity16 = 0; | |
138 | ||
139 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
140 | { | |
141 | AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc)); | |
142 | ||
143 | TParticle* particle = stack->Particle(iMc); | |
144 | ||
145 | if (!particle) | |
146 | { | |
147 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc)); | |
148 | continue; | |
149 | } | |
150 | ||
151 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
152 | continue; | |
153 | ||
154 | if (centralRegion) | |
155 | { | |
156 | if (TMath::Abs(particle->Eta()) < 1.05) | |
157 | multiplicity21++; | |
158 | if (TMath::Abs(particle->Eta()) < 0.8) | |
159 | multiplicity16++; | |
160 | } | |
161 | else | |
162 | { | |
163 | if (TMath::Abs(particle->Eta()) < 2.1) | |
164 | multiplicity21++; | |
165 | if (TMath::Abs(particle->Eta()) < 1.6) | |
166 | multiplicity16++; | |
167 | } | |
168 | }// end of mc particle | |
169 | ||
170 | AliRunLoader* runLoader = GetRunLoader(); | |
171 | if (!runLoader) | |
172 | { | |
173 | AliDebug(AliLog::kError, "runloader not available"); | |
174 | return kFALSE; | |
175 | } | |
176 | ||
177 | // TDirectory::TContext restores the current directory is restored when the scope ends. | |
178 | // This helps around ROOT bug #26025 and is good behaviour anyway | |
179 | TDirectory::TContext context(0); | |
180 | AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" ); | |
181 | loader->LoadDigits("READ"); | |
182 | TTree* treeD = loader->TreeD(); | |
183 | if (!treeD) | |
184 | { | |
185 | AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS"); | |
186 | return kFALSE; | |
187 | } | |
188 | ||
189 | treeD->SetBranchStatus("*", 0); | |
190 | treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1); | |
191 | treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1); | |
192 | ||
193 | TClonesArray* digits = 0; | |
194 | treeD->SetBranchAddress("ITSDigitsSPD", &digits); | |
195 | if (digits); | |
196 | digits->Clear(); | |
197 | ||
198 | // each value for both layers | |
199 | Int_t totalNumberOfFO[2]; | |
200 | Int_t chipsHitByPrimaries[2]; | |
201 | //Int_t chipsHitBySecondaries[2]; | |
202 | ||
203 | for (Int_t i=0; i<2; ++i) | |
204 | { | |
205 | totalNumberOfFO[i] = 0; | |
206 | chipsHitByPrimaries[i] = 0; | |
207 | //chipsHitBySecondaries[i] = 0; | |
208 | } | |
209 | ||
210 | Int_t startSPD = 0; //geom->GetStartSPD(); | |
211 | Int_t lastSPD = 239; //geom->GetLastSPD(); | |
212 | ||
213 | //printf("%d %d\n", startSPD, lastSPD); | |
214 | // for (Int_t l=0; l<240; ++l) { AliITSgeomTGeo::GetModuleId(l, i, j, k); printf("%d --> %d\n", l, i); } | |
215 | ||
216 | // loop over modules (ladders) | |
217 | for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++) | |
218 | { | |
219 | if (centralRegion) | |
220 | { | |
221 | if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3) | |
222 | continue; | |
223 | } | |
224 | ||
225 | Int_t currentLayer = 0; | |
226 | if (moduleIndex >= 80) | |
227 | currentLayer = 1; | |
228 | ||
229 | treeD->GetEvent(moduleIndex); | |
230 | ||
231 | // get number of digits in this module | |
232 | Int_t ndigitsInModule = digits->GetEntriesFast(); | |
233 | ||
234 | // get number of digits in each chip of the module | |
235 | Int_t ndigitsInChip[5]; | |
236 | Bool_t hitByPrimary[5]; | |
237 | for( Int_t iChip=0; iChip<5; iChip++) | |
238 | { | |
239 | ndigitsInChip[iChip]=0; | |
240 | hitByPrimary[iChip] = kFALSE; | |
241 | } | |
242 | ||
243 | // loop over digits in this module | |
244 | for (Int_t iDig=0; iDig<ndigitsInModule; iDig++) | |
245 | { | |
246 | AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig); | |
247 | Int_t column = dp->GetCoord1(); | |
248 | Int_t isChip = column / 32; | |
249 | ||
250 | //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip); | |
251 | ||
252 | fChipsFired->Fill(moduleIndex, isChip); | |
253 | ||
254 | ndigitsInChip[isChip]++; | |
255 | ||
256 | Bool_t debug = kFALSE; | |
257 | ||
258 | // find out which particle caused this chip to fire | |
259 | // if we find at least one primary we consider this chip being fired by a primary | |
260 | for (Int_t track=0; track<10; ++track) | |
261 | { | |
262 | Int_t label = dp->GetTrack(track); | |
263 | ||
264 | if (label < 0) | |
265 | continue; | |
266 | ||
267 | if (debug) | |
268 | printf("track %d contains label %d\n", track, label); | |
269 | ||
270 | TParticle* particle = stack->Particle(label); | |
271 | ||
272 | if (!particle) | |
273 | { | |
274 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label)); | |
275 | continue; | |
276 | } | |
277 | ||
278 | if (debug) | |
279 | { | |
280 | particle->Print(); | |
281 | printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother()); | |
282 | } | |
283 | ||
284 | // TODO delta electrons should be traced back to their mother. this is e.g. solved in AliITSClusterFinderV2::CheckLabels2 | |
285 | while (particle->P() < 0.02 && particle->GetStatusCode() == 0 && particle->GetFirstMother() >= 0) | |
286 | { | |
287 | label = particle->GetFirstMother(); | |
288 | particle = stack->Particle(label); | |
289 | ||
290 | if (!particle) | |
291 | break; | |
292 | ||
293 | if (debug) | |
294 | { | |
295 | printf("-->\n"); | |
296 | printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother()); | |
297 | particle->Print(); | |
298 | } | |
299 | } | |
300 | ||
301 | if (!particle) | |
302 | { | |
303 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label)); | |
304 | continue; | |
305 | } | |
306 | ||
307 | if (label > nPrim) | |
308 | continue; | |
309 | ||
310 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
311 | continue; | |
312 | ||
313 | if (debug) | |
314 | printf("This was a primary (or delta-electron of a primary)!\n"); | |
315 | ||
316 | hitByPrimary[isChip] = kTRUE; | |
317 | } | |
318 | } | |
319 | ||
320 | // get number of FOs in the module | |
321 | for (Int_t ifChip=0; ifChip<5; ifChip++) | |
322 | if( ndigitsInChip[ifChip] >= 1 ) | |
323 | { | |
324 | totalNumberOfFO[currentLayer]++; | |
325 | if (hitByPrimary[ifChip]) | |
326 | { | |
327 | chipsHitByPrimaries[currentLayer]++; | |
328 | } | |
329 | //else | |
330 | // chipsHitBySecondaries[currentLayer]++; | |
331 | } | |
332 | } | |
333 | ||
334 | //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2); | |
335 | ||
336 | // now find clusters | |
337 | Int_t clustersLayer[2]; | |
338 | clustersLayer[0] = 0; | |
339 | clustersLayer[1] = 0; | |
340 | ||
341 | loader->LoadRecPoints("READ"); | |
342 | TTree* treeR = loader->TreeR(); | |
343 | if (!treeR) | |
344 | { | |
345 | AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS"); | |
346 | return kFALSE; | |
347 | } | |
348 | ||
349 | // TODO branches! | |
350 | //treeR->SetBranchStatus("*", 0); | |
351 | ||
352 | TClonesArray* itsClusters = 0; | |
353 | treeR->SetBranchAddress("ITSRecPoints", &itsClusters); | |
354 | ||
355 | Int_t nTreeEntries = treeR->GetEntries(); | |
356 | for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry) | |
357 | { | |
358 | if (!treeR->GetEvent(iEntry)) | |
359 | continue; | |
360 | ||
361 | Int_t nClusters = itsClusters->GetEntriesFast(); | |
362 | ||
363 | while(nClusters--) | |
364 | { | |
365 | AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters); | |
366 | ||
367 | if (cluster->GetLayer() == 0) | |
368 | { | |
369 | clustersLayer[0]++; | |
370 | fClusterZL1->Fill(cluster->GetZ()); | |
371 | } | |
372 | else if (cluster->GetLayer() == 1) | |
373 | { | |
374 | clustersLayer[1]++; | |
375 | fClusterZL2->Fill(cluster->GetZ()); | |
376 | } | |
377 | } | |
378 | } | |
379 | ||
380 | fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]); | |
381 | fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]); | |
382 | ||
383 | return kTRUE; | |
384 | } | |
385 | ||
386 | Bool_t AliHighMultiplicitySelector::Notify() | |
387 | { | |
388 | AliRunLoader* runLoader = GetRunLoader(); | |
389 | ||
390 | if (runLoader) | |
391 | { | |
392 | AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" ); | |
393 | if (loader) | |
394 | { | |
395 | loader->UnloadDigits(); | |
396 | loader->UnloadRecPoints(); | |
397 | } | |
398 | } | |
399 | ||
400 | return AliSelectorRL::Notify(); | |
401 | } | |
402 | ||
403 | void AliHighMultiplicitySelector::SlaveTerminate() | |
404 | { | |
405 | // The SlaveTerminate() function is called after all entries or objects | |
406 | // have been processed. When running with PROOF SlaveTerminate() is called | |
407 | // on each slave server. | |
408 | ||
409 | AliSelectorRL::SlaveTerminate(); | |
410 | ||
411 | // Add the histograms to the output on each slave server | |
412 | if (!fOutput) | |
413 | { | |
414 | AliDebug(AliLog::kError, "ERROR: Output list not initialized."); | |
415 | return; | |
416 | } | |
417 | ||
418 | fOutput->Add(fChipsFired); | |
419 | fOutput->Add(fPrimaryL1); | |
420 | fOutput->Add(fPrimaryL2); | |
421 | fOutput->Add(fClusterZL1); | |
422 | fOutput->Add(fClusterZL2); | |
423 | } | |
424 | ||
425 | void AliHighMultiplicitySelector::Terminate() | |
426 | { | |
427 | // The Terminate() function is the last function to be called during | |
428 | // a query. It always runs on the client, it can be used to present | |
429 | // the results graphically or save the results to file. | |
430 | ||
431 | AliSelectorRL::Terminate(); | |
432 | ||
433 | fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired")); | |
434 | fPrimaryL1 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL1")); | |
435 | fPrimaryL2 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL2")); | |
436 | fClusterZL1 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL1")); | |
437 | fClusterZL2 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL2")); | |
438 | ||
439 | if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2) | |
440 | { | |
441 | AliError("Histograms not available"); | |
442 | return; | |
443 | } | |
444 | ||
445 | WriteHistograms(); | |
446 | } | |
447 | ||
448 | void AliHighMultiplicitySelector::WriteHistograms(const char* filename) | |
449 | { | |
450 | TFile* file = TFile::Open(filename, "RECREATE"); | |
451 | ||
452 | fChipsFired->Write(); | |
453 | fPrimaryL1->Write(); | |
454 | fPrimaryL2->Write(); | |
455 | fClusterZL1->Write(); | |
456 | fClusterZL2->Write(); | |
457 | ||
458 | file->Close(); | |
459 | } | |
460 | ||
461 | void AliHighMultiplicitySelector::ReadHistograms(const char* filename) | |
462 | { | |
463 | TFile* file = TFile::Open(filename); | |
464 | ||
465 | if (!file) | |
466 | return; | |
467 | ||
468 | fPrimaryL1 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL1")); | |
469 | fPrimaryL2 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL2")); | |
470 | fChipsFired = dynamic_cast<TH2F*> (file->Get("fChipsFired")); | |
471 | fClusterZL1 = dynamic_cast<TH1F*> (file->Get("fClusterZL1")); | |
472 | fClusterZL2 = dynamic_cast<TH1F*> (file->Get("fClusterZL2")); | |
473 | ||
474 | #define MULT 1001, -0.5, 1000.5 | |
475 | #define BINNING_LAYER1 401, -0.5, 400.5 | |
476 | #define BINNING_LAYER2 801, -0.5, 800.5 | |
477 | ||
478 | fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1); | |
479 | fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2); | |
480 | ||
481 | fL1vsL2 = new TH2F("fL1vsL2", ";Fired Chips Layer 1;Fired Chips Layer 2", BINNING_LAYER1, BINNING_LAYER2); | |
482 | fMvsL1 = new TH2F("fMvsL1", ";true multiplicity;fired chips layer1", MULT, BINNING_LAYER1); | |
483 | fMvsL2 = new TH2F("fMvsL2", ";true multiplicity;fired chips layer2", MULT, BINNING_LAYER2); | |
484 | ||
485 | fClvsL1 = new TH2F("fClvsL1", ";clusters layer1;fired chips layer1", MULT, BINNING_LAYER1); | |
486 | fClvsL2 = new TH2F("fClvsL2", ";clusters layer2;fired chips layer2", MULT, BINNING_LAYER2); | |
487 | ||
488 | for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++) | |
489 | { | |
490 | fPrimaryL1->GetEvent(i); | |
491 | fPrimaryL2->GetEvent(i); | |
492 | ||
493 | Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0]; | |
494 | Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0]; | |
495 | ||
496 | Int_t totalNumberOfFO[2]; | |
497 | totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1]; | |
498 | totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1]; | |
499 | ||
500 | Int_t chipsHitByPrimaries[2]; | |
501 | chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2]; | |
502 | chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2]; | |
503 | ||
504 | Int_t clustersLayer[2]; | |
505 | clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3]; | |
506 | clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3]; | |
507 | ||
508 | fChipsLayer1->Fill(totalNumberOfFO[0]); | |
509 | fChipsLayer2->Fill(totalNumberOfFO[1]); | |
510 | ||
511 | fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]); | |
512 | ||
513 | fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]); | |
514 | fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]); | |
515 | ||
516 | fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]); | |
517 | fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]); | |
518 | } | |
519 | } | |
520 | ||
f6093269 | 521 | TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut) |
9608df41 | 522 | { |
f6093269 | 523 | // |
524 | // returns the trigger efficiency as function of multiplicity with a given cut | |
525 | // | |
9608df41 | 526 | |
527 | //cut and multiply with x-section | |
f6093269 | 528 | TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001); |
9608df41 | 529 | //allEvents->Sumw2(); |
530 | ||
531 | //cut and multiply with x-section | |
f6093269 | 532 | TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001); |
9608df41 | 533 | //proj->Sumw2(); |
534 | ||
f6093269 | 535 | //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy(); |
536 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); | |
9608df41 | 537 | |
538 | // make probability distribution out of it | |
539 | // TODO binomial errors do not work??? weird... | |
540 | proj->Divide(proj, allEvents, 1, 1, "B"); | |
541 | ||
f6093269 | 542 | return proj; |
543 | } | |
544 | ||
545 | TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut) | |
546 | { | |
547 | // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed | |
548 | ||
549 | TH1* proj = GetTriggerEfficiency(multVsLayer, cut); | |
550 | ||
551 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); | |
9608df41 | 552 | |
553 | for (Int_t i=1; i<=proj->GetNbinsX(); i++) | |
554 | { | |
555 | if (i <= xSection->GetNbinsX()) | |
556 | { | |
557 | Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i); | |
558 | Double_t error = 0; | |
559 | ||
560 | if (value != 0) | |
561 | error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i)); | |
562 | ||
563 | proj->SetBinContent(i, value); | |
564 | proj->SetBinError(i, error); | |
565 | } | |
566 | else | |
567 | { | |
568 | proj->SetBinContent(i, 0); | |
569 | proj->SetBinError(i, 0); | |
570 | } | |
571 | } | |
572 | ||
f6093269 | 573 | //new TCanvas; proj->DrawCopy(); gPad->SetLogy(); |
9608df41 | 574 | |
575 | return proj; | |
576 | } | |
577 | ||
f6093269 | 578 | void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL) |
9608df41 | 579 | { |
f6093269 | 580 | TGraph* effGraph = new TGraph; |
581 | effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title)); | |
582 | TGraph* ratioGraph = new TGraph; | |
583 | ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title)); | |
584 | TGraph* totalGraph = new TGraph; | |
585 | totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title)); | |
586 | ||
587 | for (Int_t cut = 0; cut <= 300; cut+=50) | |
588 | { | |
589 | TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone"); | |
590 | ||
591 | //proj->Rebin(3); | |
592 | //proj->Scale(1.0 / 3); | |
593 | ||
594 | new TCanvas; proj->DrawCopy(); | |
9608df41 | 595 | |
f6093269 | 596 | Int_t limitBin = proj->GetNbinsX()+1; |
597 | while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5) | |
598 | limitBin--; | |
599 | ||
600 | Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin); | |
601 | ||
602 | effGraph->SetPoint(effGraph->GetN(), cut, limit); | |
603 | ||
604 | proj = GetXSectionCut(xSection, fMvsL, cut); | |
605 | ||
606 | Double_t ratio = 0; | |
607 | Double_t total = 0; | |
608 | if (proj->Integral(1, 1001) > 0) | |
609 | { | |
610 | ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001); | |
611 | total = proj->Integral(proj->FindBin(limit), 1001); | |
612 | } | |
613 | ||
614 | ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio); | |
615 | totalGraph->SetPoint(totalGraph->GetN(), cut, total); | |
616 | ||
617 | 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); | |
618 | } | |
619 | ||
620 | TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800); | |
621 | canvas->Divide(2, 2); | |
622 | ||
623 | canvas->cd(1); | |
624 | effGraph->Draw("A*"); | |
625 | ||
626 | for (Int_t i=8; i<=10; ++i) | |
627 | { | |
628 | TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i); | |
629 | line->Draw(); | |
630 | } | |
631 | ||
632 | canvas->cd(2); ratioGraph->Draw("A*"); | |
633 | canvas->cd(3); gPad->SetLogy(); totalGraph->Draw("A*"); | |
634 | ||
635 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
636 | } | |
637 | ||
638 | void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit) | |
639 | { | |
9608df41 | 640 | // relative x-section, once we have a collision |
641 | xSection->Scale(1.0 / xSection->Integral()); | |
642 | ||
643 | TGraph* ratioGraph = new TGraph; | |
644 | ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit)); | |
645 | TGraph* totalGraph = new TGraph; | |
646 | totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit)); | |
647 | ||
648 | Double_t max = 0; | |
649 | Int_t bestCut = -1; | |
650 | Double_t bestRatio = -1; | |
651 | Double_t bestTotal = -1; | |
652 | Int_t fullCut = -1; | |
653 | Double_t fullRatio = -1; | |
654 | Double_t fullTotal = -1; | |
655 | ||
656 | fMvsL->Sumw2(); | |
657 | ||
658 | for (Int_t cut = 50; cut <= 300; cut+=2) | |
659 | { | |
660 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
661 | ||
662 | Double_t ratio = 0; | |
663 | Double_t total = 0; | |
664 | if (proj->Integral(1, 1000) > 0) | |
665 | { | |
666 | ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000); | |
667 | total = proj->Integral(limit, 1000); | |
668 | } | |
669 | ||
670 | max = TMath::Max(max, total); | |
671 | ||
672 | //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio); | |
673 | ||
674 | if (total < max * 0.9 && bestCut == -1) | |
675 | { | |
676 | bestCut = cut; | |
677 | bestRatio = ratio; | |
678 | bestTotal = total; | |
679 | } | |
680 | ||
681 | if (ratio == 1 && fullCut == -1) | |
682 | { | |
683 | fullCut = cut; | |
684 | fullRatio = ratio; | |
685 | fullTotal = total; | |
686 | } | |
687 | ||
688 | ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio); | |
689 | totalGraph->SetPoint(totalGraph->GetN(), cut, total); | |
690 | } | |
691 | ||
692 | if (bestCut != -1) | |
693 | printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio); | |
694 | if (fullCut != -1) | |
695 | printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio); | |
696 | ||
697 | TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400); | |
698 | ratioGraph->Draw("A*"); | |
699 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
700 | ||
701 | canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400); | |
702 | totalGraph->Draw("A*"); | |
703 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
704 | } | |
705 | ||
706 | void AliHighMultiplicitySelector::JPRPlots() | |
707 | { | |
708 | /* | |
709 | ||
710 | gSystem->Load("libPWG0base"); | |
711 | .L AliHighMultiplicitySelector.cxx+ | |
712 | x = new AliHighMultiplicitySelector(); | |
713 | x->ReadHistograms("highmult_hijing100k.root"); | |
714 | x->JPRPlots(); | |
715 | ||
716 | */ | |
717 | ||
718 | // get x-sections | |
719 | TFile* file = TFile::Open("crosssectionEx.root"); | |
720 | if (!file) | |
721 | return; | |
722 | ||
723 | TH1* xSections[2]; | |
724 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
725 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
726 | ||
727 | for (Int_t i=0; i<2; ++i) | |
728 | { | |
729 | if (!xSections[i]) | |
730 | continue; | |
731 | ||
732 | TH1* xSection = xSections[i]; | |
733 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
734 | //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean | |
735 | //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean | |
736 | Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean | |
737 | ||
f6093269 | 738 | // limit is N times the mean |
739 | Int_t limit = (Int_t) (xSection->GetMean() * 10); | |
740 | ||
741 | // 10^28 lum --> 1.2 kHz | |
742 | // 10^31 lum --> 1200 kHz | |
743 | Float_t rate = 1200e3; | |
744 | ||
745 | // time in s | |
746 | Float_t lengthRun = 1e6; | |
747 | ||
9608df41 | 748 | xSection->SetStats(kFALSE); |
f6093269 | 749 | xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2"); |
9608df41 | 750 | xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5)); |
751 | xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 400 : 350); | |
f6093269 | 752 | //xSection->GetYaxis()->SetTitle("relative cross section"); |
9608df41 | 753 | xSection->GetYaxis()->SetTitleOffset(1.2); |
754 | ||
9608df41 | 755 | // relative x-section, once we have a collision |
756 | xSection->Scale(1.0 / xSection->Integral()); | |
757 | ||
758 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
759 | ||
760 | Double_t ratio = 0; | |
761 | Double_t total = 0; | |
762 | if (proj->Integral(1, 1000) > 0) | |
763 | { | |
764 | ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000); | |
765 | total = proj->Integral(limit, 1000); | |
766 | } | |
767 | ||
768 | printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio); | |
769 | ||
770 | TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600); | |
771 | canvas->SetLogy(); | |
772 | xSection->DrawCopy(); | |
773 | proj->SetLineColor(2); | |
774 | proj->SetStats(kFALSE); | |
775 | proj->DrawCopy("SAME"); | |
776 | ||
777 | TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3); | |
778 | legend->SetFillColor(0); | |
779 | legend->AddEntry(xSection, "no trigger"); | |
780 | legend->AddEntry(proj, Form("FO trigger > %d chips", cut)); | |
781 | legend->Draw(); | |
782 | ||
783 | TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2); | |
784 | line->SetLineWidth(2); | |
785 | line->Draw(); | |
786 | ||
787 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
788 | ||
789 | TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600); | |
f6093269 | 790 | //canvas2->SetTopMargin(0.05); |
791 | //canvas2->SetRightMargin(0.05); | |
9608df41 | 792 | canvas2->SetLogy(); |
f6093269 | 793 | xSection->DrawCopy("HIST"); |
9608df41 | 794 | |
f6093269 | 795 | TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3); |
9608df41 | 796 | legend2->SetFillColor(0); |
797 | legend2->AddEntry(xSection, "no trigger"); | |
798 | ||
799 | TH1* proj2 = (TH1*) proj->Clone("random"); | |
800 | proj2->Reset(); | |
f6093269 | 801 | // MB lengthRun s 100 Hz |
802 | Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000)); | |
9608df41 | 803 | proj2->FillRandom(proj, nTrigger); |
804 | ||
805 | TH1* proj3 = (TH1*) proj2->Clone("random_clone"); | |
806 | proj3->Divide(proj); | |
807 | proj3->Fit("pol0", "0", ""); | |
808 | proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0)); | |
809 | ||
f6093269 | 810 | /* |
9608df41 | 811 | proj2->DrawCopy("SAME"); |
812 | legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio))); | |
f6093269 | 813 | */ |
9608df41 | 814 | |
815 | proj2 = (TH1*) proj->Clone("random2"); | |
816 | proj2->Reset(); | |
f6093269 | 817 | // 10^31 lum --> 1200 kHz; lengthRun s |
818 | nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun); | |
9608df41 | 819 | proj2->FillRandom(proj, nTrigger); |
820 | ||
821 | proj3 = (TH1*) proj2->Clone("random_clone2"); | |
822 | proj3->Divide(proj); | |
823 | proj3->Fit("pol0", "0", ""); | |
824 | proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0)); | |
825 | ||
826 | proj2->SetLineColor(4); | |
f6093269 | 827 | proj2->SetMarkerStyle(7); |
828 | proj2->SetMarkerColor(4); | |
829 | proj2->DrawCopy("SAME P"); | |
830 | //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio))); | |
831 | legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut)); | |
9608df41 | 832 | |
833 | legend2->Draw(); | |
834 | line->Draw(); | |
835 | ||
836 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
f6093269 | 837 | canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); |
9608df41 | 838 | } |
839 | } | |
840 | ||
f91bc12d | 841 | void AliHighMultiplicitySelector::Ntrigger() |
842 | { | |
843 | // | |
844 | // produces a spectrum created with N triggers | |
845 | // number of triggers and thresholds for the moment fixed | |
846 | // | |
847 | ||
848 | /* | |
849 | ||
850 | gSystem->Load("libPWG0base"); | |
851 | .L AliHighMultiplicitySelector.cxx+g | |
852 | x = new AliHighMultiplicitySelector(); | |
853 | x->ReadHistograms("highmult_hijing100k.root"); | |
854 | x->Ntrigger(); | |
855 | ||
856 | */ | |
857 | ||
858 | // get x-sections | |
859 | TFile* file = TFile::Open("crosssectionEx.root"); | |
860 | if (!file) | |
861 | return; | |
862 | ||
863 | TH1* xSections[2]; | |
864 | xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
865 | xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
866 | ||
867 | // 10^28 lum --> 1.2 kHz | |
868 | // 10^31 lum --> 1200 kHz | |
869 | //Float_t rate = 1200e3; | |
870 | Float_t rate = 1200e3; | |
871 | ||
872 | // time in s | |
873 | Float_t lengthRun = 1e6; | |
874 | ||
875 | Int_t colors[] = { 2, 3, 4, 6, 7, 8 }; | |
876 | Int_t markers[] = { 7, 2, 4, 5, 6, 27 }; | |
877 | ||
878 | // put to 2 for second layer | |
879 | for (Int_t i=0; i<1; ++i) | |
880 | { | |
881 | if (!xSections[i]) | |
882 | continue; | |
883 | ||
884 | TH1* xSection = xSections[i]; | |
885 | TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2; | |
886 | ||
887 | Int_t nCuts = 6; | |
888 | Int_t cuts[] = { 0, 164, 178, 190, 204, 216 }; | |
889 | //Int_t nCuts = 4; | |
890 | //Int_t cuts[] = { 0, 164, 190, 216 }; | |
891 | // desired trigger rate in Hz | |
892 | Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 }; | |
893 | ||
894 | xSection->SetStats(kFALSE); | |
895 | xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2"); | |
896 | xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5)); | |
897 | xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350); | |
898 | //xSection->GetYaxis()->SetTitle("relative cross section"); | |
899 | xSection->GetYaxis()->SetTitleOffset(1.2); | |
900 | ||
901 | // relative x-section, once we have a collision | |
902 | xSection->Scale(1.0 / xSection->Integral()); | |
903 | ||
904 | TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600); | |
905 | canvas2->SetTopMargin(0.05); | |
906 | canvas2->SetRightMargin(0.05); | |
907 | canvas2->SetLogy(); | |
908 | xSection->DrawCopy("HIST"); | |
909 | ||
910 | TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4); | |
911 | legend2->SetFillColor(0); | |
912 | legend2->AddEntry(xSection, "cross-section"); | |
913 | ||
914 | for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut) | |
915 | { | |
916 | Int_t cut = cuts[currentCut]; | |
917 | ||
918 | TH1* proj = GetXSectionCut(xSection, fMvsL, cut); | |
919 | ||
920 | Double_t total = 0; | |
921 | if (proj->Integral(1, 1000) > 0) | |
922 | total = proj->Integral(1, 1000); | |
923 | ||
924 | printf("Cut at %d: rel. x-section = %e\n", cut, total); | |
925 | ||
926 | TH1* proj2 = (TH1*) proj->Clone("random2"); | |
927 | proj2->Reset(); | |
928 | ||
929 | // calculate downscale factor | |
930 | Float_t normalRate = rate * proj->Integral(1, 1000); | |
931 | Float_t downScale = normalRate / ratePerTrigger[currentCut]; | |
932 | if (downScale < 1) | |
933 | downScale = 1; | |
934 | Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun); | |
935 | ||
936 | Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger); | |
937 | proj2->FillRandom(proj, nTrigger); | |
938 | ||
939 | Int_t removed = 0; | |
940 | for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin) | |
941 | if (proj2->GetBinContent(bin) < 5) | |
942 | { | |
943 | removed += (Int_t) proj2->GetBinContent(bin); | |
944 | proj2->SetBinContent(bin, 0); | |
945 | } | |
946 | ||
947 | Printf("Removed %d events", removed); | |
948 | ||
949 | proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000)); | |
950 | ||
951 | proj2->SetLineColor(colors[currentCut]); | |
952 | proj2->SetMarkerStyle(markers[currentCut]); | |
953 | proj2->SetMarkerColor(colors[currentCut]); | |
954 | proj2->DrawCopy("SAME P"); | |
955 | legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut)); | |
956 | } | |
957 | ||
958 | legend2->Draw(); | |
959 | ||
960 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
961 | canvas2->SaveAs(Form("%s.eps", canvas2->GetName())); | |
962 | } | |
963 | } | |
964 | ||
9608df41 | 965 | void AliHighMultiplicitySelector::DrawHistograms() |
966 | { | |
967 | /* | |
968 | ||
969 | gSystem->Load("libPWG0base"); | |
970 | .L AliHighMultiplicitySelector.cxx+ | |
971 | x = new AliHighMultiplicitySelector(); | |
972 | x->ReadHistograms("highmult_pythia.root"); | |
973 | x->DrawHistograms(); | |
974 | ||
975 | ||
976 | gSystem->Load("libPWG0base"); | |
977 | .L AliHighMultiplicitySelector.cxx+ | |
978 | x = new AliHighMultiplicitySelector(); | |
979 | x->ReadHistograms("highmult_hijing.root"); | |
980 | x->DrawHistograms(); | |
981 | ||
982 | gSystem->Load("libPWG0base"); | |
983 | .L AliHighMultiplicitySelector.cxx+ | |
984 | x = new AliHighMultiplicitySelector(); | |
985 | x->ReadHistograms("highmult_central.root"); | |
986 | x->DrawHistograms(); | |
987 | ||
988 | gSystem->Load("libPWG0base"); | |
989 | .L AliHighMultiplicitySelector.cxx+ | |
990 | x = new AliHighMultiplicitySelector(); | |
991 | x->ReadHistograms("highmult_hijing100k.root"); | |
992 | x->DrawHistograms(); | |
993 | ||
994 | */ | |
995 | ||
996 | /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400); | |
997 | ||
998 | fChipsLayer2->SetLineColor(2); | |
999 | fChipsLayer2->SetStats(kFALSE); | |
1000 | fChipsLayer1->SetStats(kFALSE); | |
1001 | fChipsLayer2->SetTitle(""); | |
1002 | fChipsLayer2->DrawCopy(); | |
1003 | fChipsLayer1->DrawCopy("SAME"); | |
1004 | canvas->SaveAs("chips.gif"); | |
1005 | ||
1006 | canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400); | |
1007 | fL1vsL2->SetStats(kFALSE); | |
1008 | fL1vsL2->DrawCopy("COLZ"); | |
1009 | gPad->SetLogz(); | |
1010 | canvas->SaveAs("L1vsL2.gif");*/ | |
1011 | ||
f6093269 | 1012 | TCanvas *canvas = new TCanvas("L1", "L1", 800, 600); |
1013 | canvas->SetTopMargin(0.05); | |
1014 | canvas->SetRightMargin(0.12); | |
9608df41 | 1015 | fMvsL1->SetStats(kFALSE); |
1016 | fMvsL1->DrawCopy("COLZ"); | |
1017 | gPad->SetLogz(); | |
1018 | ||
1019 | canvas->SaveAs("L1NoCurve.gif"); | |
f6093269 | 1020 | canvas->SaveAs("L1NoCurve.eps"); |
9608df41 | 1021 | |
1022 | // draw corresponding theoretical curve | |
1023 | TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000); | |
1024 | func->SetParameter(0, 400-5*2); | |
1025 | func->DrawCopy("SAME"); | |
1026 | ||
1027 | canvas->SaveAs("L1.gif"); | |
1028 | ||
1029 | canvas = new TCanvas("L2", "L2", 600, 400); | |
1030 | //fMvsL2->GetYaxis()->SetRangeUser(0, 150); | |
1031 | fMvsL2->SetStats(kFALSE); | |
1032 | fMvsL2->DrawCopy("COLZ"); | |
1033 | gPad->SetLogz(); | |
1034 | func->SetParameter(0, 800-5*4); | |
1035 | func->DrawCopy("SAME"); | |
1036 | canvas->SaveAs("L2.gif"); | |
1037 | ||
f6093269 | 1038 | // get x-sections |
1039 | TFile* file = TFile::Open("crosssectionEx.root"); | |
1040 | if (file) | |
1041 | { | |
1042 | TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex")); | |
1043 | TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex")); | |
1044 | ||
1045 | MakeGraphs2("Layer1", xSection2, fMvsL1); | |
1046 | return; | |
1047 | ||
1048 | // 5 times the mean | |
1049 | //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2); | |
1050 | //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2)); | |
1051 | ||
1052 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8)); | |
1053 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8)); | |
1054 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9)); | |
1055 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9)); | |
1056 | MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10)); | |
1057 | MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10)); | |
1058 | ||
1059 | file->Close(); | |
1060 | } | |
1061 | ||
9608df41 | 1062 | // make spread hists |
1063 | TGraph* spread = new TGraph; | |
1064 | spread->SetTitle("Spread L1;true multiplicity;RMS"); | |
1065 | ||
1066 | for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i) | |
1067 | { | |
1068 | TH1* proj = fMvsL1->ProjectionY("proj", i, i); | |
1069 | spread->SetPoint(spread->GetN(), i, proj->GetRMS()); | |
1070 | } | |
1071 | ||
1072 | canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400); | |
1073 | spread->Draw("A*"); | |
1074 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1075 | ||
1076 | TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150); | |
1077 | log->SetParLimits(0, 0, 10); | |
1078 | log->SetParLimits(1, 1e-5, 10); | |
1079 | ||
1080 | spread->Fit(log, "", "", 1, 150); | |
1081 | log->DrawCopy("SAME"); | |
1082 | ||
1083 | TGraph* spread2 = new TGraph; | |
1084 | spread2->SetTitle("Spread L2;true multiplicity;RMS"); | |
1085 | ||
1086 | for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i) | |
1087 | { | |
1088 | TH1* proj = fMvsL2->ProjectionY("proj", i, i); | |
1089 | spread2->SetPoint(spread2->GetN(), i, proj->GetRMS()); | |
1090 | } | |
1091 | ||
1092 | canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400); | |
1093 | spread2->Draw("A*"); | |
1094 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1095 | ||
1096 | spread2->Fit(log, "", "", 1, 150); | |
1097 | log->DrawCopy("SAME"); | |
1098 | ||
9608df41 | 1099 | canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400); |
1100 | fClvsL1->SetStats(kFALSE); | |
1101 | fClvsL1->DrawCopy("COLZ"); | |
1102 | gPad->SetLogz(); | |
1103 | ||
1104 | func->SetParameter(0, 400-5*2); | |
1105 | func->DrawCopy("SAME"); | |
1106 | ||
1107 | canvas->SaveAs("Clusters_L1.gif"); | |
1108 | ||
1109 | canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400); | |
1110 | fClvsL2->SetStats(kFALSE); | |
1111 | fClvsL2->DrawCopy("COLZ"); | |
1112 | gPad->SetLogz(); | |
1113 | func->SetParameter(0, 800-5*4); | |
1114 | func->DrawCopy("SAME"); | |
1115 | canvas->SaveAs("Clusters_L2.gif"); | |
1116 | ||
1117 | canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400); | |
1118 | //fChipsFired->GetYaxis()->SetRangeUser(0, 150); | |
1119 | fChipsFired->SetStats(kFALSE); | |
1120 | fChipsFired->DrawCopy("COLZ"); | |
1121 | canvas->SaveAs("ChipsFired.gif"); | |
1122 | ||
1123 | /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1); | |
1124 | TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2); | |
1125 | ||
1126 | for (Int_t treshold = 0; treshold < 800; treshold++) | |
1127 | { | |
1128 | if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0) | |
1129 | { | |
1130 | TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram()); | |
1131 | if (mult) | |
1132 | tresholdHistL1->Fill(treshold, mult->GetMean()); | |
1133 | } | |
1134 | if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0) | |
1135 | { | |
1136 | TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram()); | |
1137 | if (mult) | |
1138 | tresholdHistL2->Fill(treshold, mult->GetMean()); | |
1139 | } | |
1140 | } | |
1141 | ||
1142 | canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400); | |
1143 | tresholdHistL1->Draw(); | |
1144 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1145 | ||
1146 | canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400); | |
1147 | tresholdHistL2->Draw(); | |
1148 | canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/ | |
1149 | ||
1150 | fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff"); | |
1151 | fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff"); | |
1152 | ||
1153 | canvas = new TCanvas("Efficiency", "Efficiency", 600, 400); | |
1154 | fPrimaryL1->GetHistogram()->SetStats(kFALSE); | |
1155 | fPrimaryL1->GetHistogram()->Draw(); | |
1156 | fPrimaryL2->GetHistogram()->SetLineColor(2); | |
1157 | fPrimaryL2->GetHistogram()->SetStats(kFALSE); | |
1158 | fPrimaryL2->GetHistogram()->Draw("SAME"); | |
1159 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1160 | ||
1161 | canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400); | |
1162 | fClusterZL1->Rebin(2); | |
1163 | fClusterZL1->Draw(); | |
1164 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1165 | ||
1166 | canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400); | |
1167 | fClusterZL2->Draw(); | |
1168 | fClusterZL2->Rebin(2); | |
1169 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
1170 | } |