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