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