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