PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx

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