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(kFALSE)
  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::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
 522 {
 523   //
 524   // returns the trigger efficiency as function of multiplicity with a given cut
 525   //
 526
 527   //cut and multiply with x-section
 528   TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
 529   //allEvents->Sumw2();
 530
 531   //cut and multiply with x-section
 532   TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
 533   //proj->Sumw2();
 534
 535   //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
 536   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
 537
 538   // make probability distribution out of it
 539   // TODO binomial errors do not work??? weird...
 540   proj->Divide(proj, allEvents, 1, 1, "B");
 541
 542   return proj;
 543 }
 544
 545 TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut)
 546 {
 547   // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed
 548
 549   TH1* proj = GetTriggerEfficiency(multVsLayer, cut);
 550
 551   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
 552
 553   for (Int_t i=1; i<=proj->GetNbinsX(); i++)
 554   {
 555     if (i <= xSection->GetNbinsX())
 556     {
 557       Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i);
 558       Double_t error = 0;
 559
 560       if (value != 0)
 561         error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i));
 562
 563       proj->SetBinContent(i, value);
 564       proj->SetBinError(i, error);
 565     }
 566     else
 567     {
 568       proj->SetBinContent(i, 0);
 569       proj->SetBinError(i, 0);
 570     }
 571   }
 572
 573   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
 574
 575   return proj;
 576 }
 577
 578 void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
 579 {
 580   TGraph* effGraph = new TGraph;
 581   effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title));
 582   TGraph* ratioGraph = new TGraph;
 583   ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title));
 584   TGraph* totalGraph = new TGraph;
 585   totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title));
 586
 587   for (Int_t cut = 0; cut <= 300; cut+=50)
 588   {
 589     TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone");
 590
 591     //proj->Rebin(3);
 592     //proj->Scale(1.0 / 3);
 593
 594     new TCanvas; proj->DrawCopy();
 595
 596     Int_t limitBin = proj->GetNbinsX()+1;
 597     while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5)
 598       limitBin--;
 599
 600     Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin);
 601
 602     effGraph->SetPoint(effGraph->GetN(), cut, limit);
 603
 604     proj = GetXSectionCut(xSection, fMvsL, cut);
 605
 606     Double_t ratio = 0;
 607     Double_t total = 0;
 608     if (proj->Integral(1, 1001) > 0)
 609     {
 610       ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001);
 611       total = proj->Integral(proj->FindBin(limit), 1001);
 612     }
 613
 614     ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
 615     totalGraph->SetPoint(totalGraph->GetN(), cut, total);
 616
 617     Printf("Cut at %d --> trigger eff. is > 0.5 for mult. >= %.2f. That is the case for %f of the triggered, %e of all events", cut, limit, ratio, total);
 618   }
 619
 620   TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800);
 621   canvas->Divide(2, 2);
 622
 623   canvas->cd(1);
 624   effGraph->Draw("A*");
 625
 626   for (Int_t i=8; i<=10; ++i)
 627   {
 628     TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i);
 629     line->Draw();
 630   }
 631
 632   canvas->cd(2);  ratioGraph->Draw("A*");
 633   canvas->cd(3);  gPad->SetLogy(); totalGraph->Draw("A*");
 634
 635   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
 636 }
 637
 638 void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
 639 {
 640   // relative x-section, once we have a collision
 641   xSection->Scale(1.0 / xSection->Integral());
 642
 643   TGraph* ratioGraph = new TGraph;
 644   ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit));
 645   TGraph* totalGraph = new TGraph;
 646   totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit));
 647
 648   Double_t max = 0;
 649   Int_t bestCut = -1;
 650   Double_t bestRatio = -1;
 651   Double_t bestTotal = -1;
 652   Int_t fullCut = -1;
 653   Double_t fullRatio = -1;
 654   Double_t fullTotal = -1;
 655
 656   fMvsL->Sumw2();
 657
 658   for (Int_t cut = 50; cut <= 300; cut+=2)
 659   {
 660     TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
 661
 662     Double_t ratio = 0;
 663     Double_t total = 0;
 664     if (proj->Integral(1, 1000) > 0)
 665     {
 666       ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
 667       total = proj->Integral(limit, 1000);
 668     }
 669
 670     max = TMath::Max(max, total);
 671
 672     //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
 673
 674     if (total < max * 0.9 && bestCut == -1)
 675     {
 676       bestCut = cut;
 677       bestRatio = ratio;
 678       bestTotal = total;
 679     }
 680
 681     if (ratio == 1 && fullCut == -1)
 682     {
 683       fullCut = cut;
 684       fullRatio = ratio;
 685       fullTotal = total;
 686     }
 687
 688     ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
 689     totalGraph->SetPoint(totalGraph->GetN(), cut, total);
 690   }
 691
 692   if (bestCut != -1)
 693     printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio);
 694   if (fullCut != -1)
 695     printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio);
 696
 697   TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400);
 698   ratioGraph->Draw("A*");
 699   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
 700
 701   canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400);
 702   totalGraph->Draw("A*");
 703   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
 704 }
 705
 706 void AliHighMultiplicitySelector::JPRPlots()
 707 {
 708   /*
 709
 710   gSystem->Load("libPWG0base");
 711   .L AliHighMultiplicitySelector.cxx+
 712   x = new AliHighMultiplicitySelector();
 713   x->ReadHistograms("highmult_hijing100k.root");
 714   x->JPRPlots();
 715
 716   */
 717
 718   // get x-sections
 719   TFile* file = TFile::Open("crosssectionEx.root");
 720   if (!file)
 721     return;
 722
 723   TH1* xSections[2];
 724   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
 725   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
 726
 727   for (Int_t i=0; i<2; ++i)
 728   {
 729     if (!xSections[i])
 730       continue;
 731
 732     TH1* xSection = xSections[i];
 733     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
 734     //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean
 735     //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean
 736     Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean
 737
 738     // limit is N times the mean
 739     Int_t limit = (Int_t) (xSection->GetMean() * 10);
 740
 741     // 10^28 lum --> 1.2 kHz
 742     // 10^31 lum --> 1200 kHz
 743     Float_t rate = 1200e3;
 744
 745     // time in s
 746     Float_t lengthRun = 1e6;
 747
 748     xSection->SetStats(kFALSE);
 749     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
 750     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
 751     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 400 : 350);
 752     //xSection->GetYaxis()->SetTitle("relative cross section");
 753     xSection->GetYaxis()->SetTitleOffset(1.2);
 754
 755     // relative x-section, once we have a collision
 756     xSection->Scale(1.0 / xSection->Integral());
 757
 758     TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
 759
 760     Double_t ratio = 0;
 761     Double_t total = 0;
 762     if (proj->Integral(1, 1000) > 0)
 763     {
 764       ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
 765       total = proj->Integral(limit, 1000);
 766     }
 767
 768     printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
 769
 770     TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600);
 771     canvas->SetLogy();
 772     xSection->DrawCopy();
 773     proj->SetLineColor(2);
 774     proj->SetStats(kFALSE);
 775     proj->DrawCopy("SAME");
 776
 777     TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3);
 778     legend->SetFillColor(0);
 779     legend->AddEntry(xSection, "no trigger");
 780     legend->AddEntry(proj, Form("FO trigger > %d chips", cut));
 781     legend->Draw();
 782
 783     TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2);
 784     line->SetLineWidth(2);
 785     line->Draw();
 786
 787     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
 788
 789     TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600);
 790     //canvas2->SetTopMargin(0.05);
 791     //canvas2->SetRightMargin(0.05);
 792     canvas2->SetLogy();
 793     xSection->DrawCopy("HIST");
 794
 795     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
 796     legend2->SetFillColor(0);
 797     legend2->AddEntry(xSection, "no trigger");
 798
 799     TH1* proj2 = (TH1*) proj->Clone("random");
 800     proj2->Reset();
 801     // MB lengthRun s 100 Hz
 802     Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
 803     proj2->FillRandom(proj, nTrigger);
 804
 805     TH1* proj3 = (TH1*) proj2->Clone("random_clone");
 806     proj3->Divide(proj);
 807     proj3->Fit("pol0", "0", "");
 808     proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
 809
 810     /*
 811     proj2->DrawCopy("SAME");
 812     legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
 813     */
 814
 815     proj2 = (TH1*) proj->Clone("random2");
 816     proj2->Reset();
 817     // 10^31 lum --> 1200 kHz; lengthRun s
 818     nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
 819     proj2->FillRandom(proj, nTrigger);
 820
 821     proj3 = (TH1*) proj2->Clone("random_clone2");
 822     proj3->Divide(proj);
 823     proj3->Fit("pol0", "0", "");
 824     proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
 825
 826     proj2->SetLineColor(4);
 827     proj2->SetMarkerStyle(7);
 828     proj2->SetMarkerColor(4);
 829     proj2->DrawCopy("SAME P");
 830     //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
 831     legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut));
 832
 833     legend2->Draw();
 834     line->Draw();
 835
 836     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
 837     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
 838   }
 839 }
 840
 841 void AliHighMultiplicitySelector::Ntrigger()
 842 {
 843   //
 844   // produces a spectrum created with N triggers
 845   // number of triggers and thresholds for the moment fixed
 846   //
 847
 848   /*
 849
 850   gSystem->Load("libPWG0base");
 851   .L AliHighMultiplicitySelector.cxx+g
 852   x = new AliHighMultiplicitySelector();
 853   x->ReadHistograms("highmult_hijing100k.root");
 854   x->Ntrigger();
 855
 856   */
 857
 858   // get x-sections
 859   TFile* file = TFile::Open("crosssectionEx.root");
 860   if (!file)
 861     return;
 862
 863   TH1* xSections[2];
 864   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
 865   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
 866
 867   // 10^28 lum --> 1.2 kHz
 868   // 10^31 lum --> 1200 kHz
 869   //Float_t rate = 1200e3;
 870   Float_t rate = 1200e3;
 871
 872   // time in s
 873   Float_t lengthRun = 1e6;
 874
 875   Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
 876   Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
 877
 878   // put to 2 for second layer
 879   for (Int_t i=0; i<1; ++i)
 880   {
 881     if (!xSections[i])
 882       continue;
 883
 884     TH1* xSection = xSections[i];
 885     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
 886
 887     Int_t nCuts = 6;
 888     Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
 889     //Int_t nCuts = 4;
 890     //Int_t cuts[] = { 0, 164, 190, 216 };
 891     // desired trigger rate in Hz
 892     Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
 893
 894     xSection->SetStats(kFALSE);
 895     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
 896     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
 897     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
 898     //xSection->GetYaxis()->SetTitle("relative cross section");
 899     xSection->GetYaxis()->SetTitleOffset(1.2);
 900
 901     // relative x-section, once we have a collision
 902     xSection->Scale(1.0 / xSection->Integral());
 903
 904     TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
 905     canvas2->SetTopMargin(0.05);
 906     canvas2->SetRightMargin(0.05);
 907     canvas2->SetLogy();
 908     xSection->DrawCopy("HIST");
 909
 910     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
 911     legend2->SetFillColor(0);
 912     legend2->AddEntry(xSection, "cross-section");
 913
 914     for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
 915     {
 916       Int_t cut = cuts[currentCut];
 917
 918       TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
 919
 920       Double_t total = 0;
 921       if (proj->Integral(1, 1000) > 0)
 922         total = proj->Integral(1, 1000);
 923
 924       printf("Cut at %d: rel. x-section = %e\n", cut, total);
 925
 926       TH1* proj2 = (TH1*) proj->Clone("random2");
 927       proj2->Reset();
 928
 929       // calculate downscale factor
 930       Float_t normalRate = rate * proj->Integral(1, 1000);
 931       Float_t downScale = normalRate / ratePerTrigger[currentCut];
 932       if (downScale < 1)
 933         downScale = 1;
 934       Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
 935
 936       Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
 937       proj2->FillRandom(proj, nTrigger);
 938
 939       Int_t removed = 0;
 940       for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
 941         if (proj2->GetBinContent(bin) < 5)
 942         {
 943           removed += (Int_t) proj2->GetBinContent(bin);
 944           proj2->SetBinContent(bin, 0);
 945         }
 946
 947       Printf("Removed %d events", removed);
 948
 949       proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
 950
 951       proj2->SetLineColor(colors[currentCut]);
 952       proj2->SetMarkerStyle(markers[currentCut]);
 953       proj2->SetMarkerColor(colors[currentCut]);
 954       proj2->DrawCopy("SAME P");
 955       legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
 956     }
 957
 958     legend2->Draw();
 959
 960     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
 961     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
 962   }
 963 }
 964
 965 void AliHighMultiplicitySelector::DrawHistograms()
 966 {
 967   /*
 968
 969   gSystem->Load("libPWG0base");
 970   .L AliHighMultiplicitySelector.cxx+
 971   x = new AliHighMultiplicitySelector();
 972   x->ReadHistograms("highmult_pythia.root");
 973   x->DrawHistograms();
 974
 975
 976   gSystem->Load("libPWG0base");
 977   .L AliHighMultiplicitySelector.cxx+
 978   x = new AliHighMultiplicitySelector();
 979   x->ReadHistograms("highmult_hijing.root");
 980   x->DrawHistograms();
 981
 982   gSystem->Load("libPWG0base");
 983   .L AliHighMultiplicitySelector.cxx+
 984   x = new AliHighMultiplicitySelector();
 985   x->ReadHistograms("highmult_central.root");
 986   x->DrawHistograms();
 987
 988   gSystem->Load("libPWG0base");
 989   .L AliHighMultiplicitySelector.cxx+
 990   x = new AliHighMultiplicitySelector();
 991   x->ReadHistograms("highmult_hijing100k.root");
 992   x->DrawHistograms();
 993
 994   */
 995
 996   /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
 997
 998   fChipsLayer2->SetLineColor(2);
 999   fChipsLayer2->SetStats(kFALSE);
1000   fChipsLayer1->SetStats(kFALSE);
1001   fChipsLayer2->SetTitle("");
1002   fChipsLayer2->DrawCopy();
1003   fChipsLayer1->DrawCopy("SAME");
1004   canvas->SaveAs("chips.gif");
1005
1006   canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1007   fL1vsL2->SetStats(kFALSE);
1008   fL1vsL2->DrawCopy("COLZ");
1009   gPad->SetLogz();
1010   canvas->SaveAs("L1vsL2.gif");*/
1011
1012   TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1013   canvas->SetTopMargin(0.05);
1014   canvas->SetRightMargin(0.12);
1015   fMvsL1->SetStats(kFALSE);
1016   fMvsL1->DrawCopy("COLZ");
1017   gPad->SetLogz();
1018
1019   canvas->SaveAs("L1NoCurve.gif");
1020   canvas->SaveAs("L1NoCurve.eps");
1021
1022   // draw corresponding theoretical curve
1023   TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1024   func->SetParameter(0, 400-5*2);
1025   func->DrawCopy("SAME");
1026
1027   canvas->SaveAs("L1.gif");
1028
1029   canvas = new TCanvas("L2", "L2", 600, 400);
1030   //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1031   fMvsL2->SetStats(kFALSE);
1032   fMvsL2->DrawCopy("COLZ");
1033   gPad->SetLogz();
1034   func->SetParameter(0, 800-5*4);
1035   func->DrawCopy("SAME");
1036   canvas->SaveAs("L2.gif");
1037
1038   // get x-sections
1039   TFile* file = TFile::Open("crosssectionEx.root");
1040   if (file)
1041   {
1042     TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1043     TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1044
1045     MakeGraphs2("Layer1", xSection2, fMvsL1);
1046     return;
1047
1048     // 5 times the mean
1049     //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1050     //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1051
1052     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1053     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1054     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1055     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1056     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1057     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1058
1059     file->Close();
1060   }
1061
1062   // make spread hists
1063   TGraph* spread = new TGraph;
1064   spread->SetTitle("Spread L1;true multiplicity;RMS");
1065
1066   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1067   {
1068     TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1069     spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1070   }
1071
1072   canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1073   spread->Draw("A*");
1074   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1075
1076   TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1077   log->SetParLimits(0, 0, 10);
1078   log->SetParLimits(1, 1e-5, 10);
1079
1080   spread->Fit(log, "", "", 1, 150);
1081   log->DrawCopy("SAME");
1082
1083   TGraph* spread2 = new TGraph;
1084   spread2->SetTitle("Spread L2;true multiplicity;RMS");
1085
1086   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1087   {
1088     TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1089     spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1090   }
1091
1092   canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1093   spread2->Draw("A*");
1094   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1095
1096   spread2->Fit(log, "", "", 1, 150);
1097   log->DrawCopy("SAME");
1098
1099   canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1100   fClvsL1->SetStats(kFALSE);
1101   fClvsL1->DrawCopy("COLZ");
1102   gPad->SetLogz();
1103
1104   func->SetParameter(0, 400-5*2);
1105   func->DrawCopy("SAME");
1106
1107   canvas->SaveAs("Clusters_L1.gif");
1108
1109   canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1110   fClvsL2->SetStats(kFALSE);
1111   fClvsL2->DrawCopy("COLZ");
1112   gPad->SetLogz();
1113   func->SetParameter(0, 800-5*4);
1114   func->DrawCopy("SAME");
1115   canvas->SaveAs("Clusters_L2.gif");
1116
1117   canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1118   //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1119   fChipsFired->SetStats(kFALSE);
1120   fChipsFired->DrawCopy("COLZ");
1121   canvas->SaveAs("ChipsFired.gif");
1122
1123   /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1124   TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1125
1126   for (Int_t treshold = 0; treshold < 800; treshold++)
1127   {
1128     if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1129     {
1130       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1131       if (mult)
1132         tresholdHistL1->Fill(treshold, mult->GetMean());
1133     }
1134     if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1135     {
1136       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1137       if (mult)
1138         tresholdHistL2->Fill(treshold, mult->GetMean());
1139     }
1140   }
1141
1142   canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1143   tresholdHistL1->Draw();
1144   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1145
1146   canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1147   tresholdHistL2->Draw();
1148   canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1149
1150   fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1151   fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1152
1153   canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1154   fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1155   fPrimaryL1->GetHistogram()->Draw();
1156   fPrimaryL2->GetHistogram()->SetLineColor(2);
1157   fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1158   fPrimaryL2->GetHistogram()->Draw("SAME");
1159   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1160
1161   canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1162   fClusterZL1->Rebin(2);
1163   fClusterZL1->Draw();
1164   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1165
1166   canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1167   fClusterZL2->Draw();
1168   fClusterZL2->Rebin(2);
1169   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1170 }