PHOS/AliPHOSIhepAnalyze.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Algorythm class to analyze PHOSv1 events:
  18 // Construct histograms and displays them.
  19 // IHEP CPV/PHOS reconstruction algorithm used.
  20 // Use the macro EditorBar.C for best access to the functionnalities
  21 //*--
  22 //*-- Author: B. Polichtchouk (IHEP)
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 // --- ROOT system ---
  26
  27 #include "TFile.h"
  28 #include "TH1.h"
  29 #include "TPad.h"
  30 #include "TH2.h"
  31 #include "TParticle.h"
  32 #include "TClonesArray.h"
  33 #include "TTree.h"
  34 #include "TMath.h"
  35 #include "TCanvas.h"
  36 #include "TStyle.h"
  37
  38 // --- Standard library ---
  39
  40 // --- AliRoot header files ---
  41
  42 #include "AliRunLoader.h"
  43 #include "AliHeader.h"
  44
  45 // --- PHOS header files ---
  46 #include "AliPHOSIhepAnalyze.h"
  47 #include "AliPHOSDigit.h"
  48 #include "AliPHOSRecParticle.h"
  49 #include "AliPHOSLoader.h"
  50 #include "AliPHOSHit.h"
  51 #include "AliPHOSImpact.h"
  52 #include "AliPHOSvImpacts.h"
  53 #include "AliPHOSCpvRecPoint.h"
  54 #include "AliRun.h"
  55 #include "AliPHOSGeometry.h"
  56 #include "AliPHOSEvalRecPoint.h"
  57
  58 ClassImp(AliPHOSIhepAnalyze)
  59
  60 //____________________________________________________________________________
  61
  62 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze()
  63  {
  64    fRunLoader = 0x0;
  65  }
  66
  67 //____________________________________________________________________________
  68
  69 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) : fFileName(name) {
  70   // Constructor: open a header file
  71   fRunLoader = AliRunLoader::Open(fFileName);
  72   if (fRunLoader == 0x0)
  73    {
  74      Fatal("AliPHOSIhepAnalyze","Can not load event from file %s",name);
  75    }
  76  }
  77
  78 //____________________________________________________________________________
  79 void AliPHOSIhepAnalyze::AnalyzeCPV1(Int_t Nevents)
  80 {
  81   //
  82   // Analyzes CPV characteristics: resolutions, cluster multiplicity,
  83   // cluster lengths along Z and Phi.
  84   // Author: Yuri Kharlov
  85   // 9 October 2000
  86   // Modified by Boris Polichtchouk, 3.07.2001
  87   //
  88
  89   // Book histograms
  90
  91   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
  92   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
  93   TH1F *hChi2 = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 10.);
  94   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
  95   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
  96   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
  97
  98   TH1F *hEg   = new TH1F("hEg"  ,"Energies of impacts",30,0.,6.);
  99   TH1F *hEr   = new TH1F("hEr"  ,"Amplitudes of rec. points",50,0.,20.);
 100
 101   TList * fCpvImpacts ;
 102   TBranch * branchCPVimpacts;
 103
 104
 105
 106   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 107   if ( please == 0 )
 108    {
 109      Error("AnalyzeCPV1","Could not obtain the Loader object !");
 110      return ;
 111    }
 112
 113   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 114
 115   Info("AnalyzeCPV1", "Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths") ;
 116   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 117
 118     Int_t nTotalGen = 0;
 119     Int_t nChargedGen = 0;
 120
 121     Int_t ntracks = gAlice->GetEvent(ievent);
 122     Info("AnalyzeCPV1", ">>>>>>>Event %d .<<<<<<<", ievent) ;
 123
 124   /******************************************************************/
 125       TTree* treeH = please->TreeH();
 126       if (treeH == 0x0)
 127        {
 128          Error("AnalyzeCPV1","Can not get TreeH");
 129          return;
 130        }
 131 /******************************************************************/
 132
 133     // Get branch of CPV impacts
 134     if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) {
 135       Info("AnalyzeCPV1", "Couldn't find branch PHOSCpvImpacts. Exit.") ;
 136       return;
 137     }
 138
 139     // Create and fill arrays of hits for each CPV module
 140
 141     Int_t nOfModules = fGeom->GetNModules();
 142     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 143     Int_t iModule = 0;
 144     for (iModule=0; iModule < nOfModules; iModule++)
 145       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 146
 147     TClonesArray    *impacts;
 148     AliPHOSImpact   *impact;
 149     TLorentzVector   p;
 150
 151     // First go through all primary tracks and fill the arrays
 152     // of hits per each CPV module
 153
 154     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 155       branchCPVimpacts ->SetAddress(&fCpvImpacts);
 156       branchCPVimpacts ->GetEntry(itrack,0);
 157
 158       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 159         impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 160         // Do loop over impacts in the module
 161         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 162           impact=(AliPHOSImpact*)impacts->At(iImpact);
 163           hEg->Fill(impact->GetMomentum().E());
 164           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 165           if(IsCharged(impact->GetPid())) nChargedGen++;
 166           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 167         }
 168       }
 169       fCpvImpacts->Clear();
 170     }
 171     for (iModule=0; iModule < nOfModules; iModule++) {
 172       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 173       printf("CPV module %d has %d impacts\n",iModule,nsum);
 174       nTotalGen += nsum;
 175     }
 176
 177     // Then go through reconstructed points and for each find
 178     // the closeset hit
 179     // The distance from the rec.point to the closest hit
 180     // gives the coordinate resolution of the CPV
 181
 182     fRunLoader->GetEvent(ievent);
 183     TIter nextRP(please->CpvRecPoints()) ;
 184     AliPHOSCpvRecPoint *cpvRecPoint ;
 185     Float_t xgen, ygen, zgen;
 186     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
 187
 188       Float_t chi2dof = ((AliPHOSEvalRecPoint*)cpvRecPoint)->Chi2Dof();
 189       hChi2->Fill(chi2dof);
 190       hEr->Fill(cpvRecPoint->GetEnergy());
 191
 192       TVector3  locpos;
 193       cpvRecPoint->GetLocalPosition(locpos);
 194       Int_t phosModule = cpvRecPoint->GetPHOSMod();
 195       Int_t rpMult     = cpvRecPoint->GetMultiplicity();
 196       Int_t rpMultX, rpMultZ;
 197       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
 198       Float_t xrec  = locpos.X();
 199       Float_t zrec  = locpos.Z();
 200       Float_t dxmin = 1.e+10;
 201       Float_t dzmin = 1.e+10;
 202       Float_t r2min = 1.e+10;
 203       Float_t r2;
 204
 205       Int_t nCPVhits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
 206       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
 207       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
 208       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 209         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
 210         xgen   = impact->X();
 211         zgen   = impact->Z();
 212         ygen   = impact->Y();
 213
 214         //Transform to the local ref.frame
 215         const AliPHOSGeometry* geom = please->PHOSGeometry();
 216         Float_t phig = geom->GetPHOSAngle(phosModule);
 217         Float_t phi = TMath::Pi()/180*phig;
 218         Float_t distanceIPtoCPV = geom->GetIPtoOuterCoverDistance() -
 219                           (geom->GetFTPosition(1)+
 220                           geom->GetFTPosition(2)+
 221                           geom->GetCPVTextoliteThickness()
 222                           )/2;
 223         Float_t xoL,yoL,zoL ;
 224 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
 225 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoCPV;
 226         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
 227         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoCPV;
 228         zoL = zgen;
 229
 230         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
 231         if ( r2 < r2min ) {
 232           r2min = r2;
 233           dxmin = xoL - xrec;
 234           dzmin = zoL - zrec;
 235           locImpX = xoL;
 236           locImpZ = zoL;
 237           gImpX = xgen;
 238           gImpZ = zgen;
 239           gImpY = ygen;
 240         }
 241       }
 242       Info("AnalyzeCPV1", "Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY);
 243       Info("AnalyzeCPV1", "Impact local (X,Z) = %f %f", locImpX, locImpZ);
 244       Info("AnalyzeCPV1", "Reconstructed (X,Z) = %f %f", xrec, zrec);
 245       Info("AnalyzeCPV1", "dxmin %f dzmin %f", dxmin, dzmin);
 246       hDx  ->Fill(dxmin);
 247       hDz  ->Fill(dzmin);
 248 //        hDr  ->Fill(TMath::Sqrt(r2min));
 249       hNrp ->Fill(rpMult);
 250       hNrpX->Fill(rpMultX);
 251       hNrpZ->Fill(rpMultZ);
 252     }
 253     delete [] hitsPerModule;
 254
 255     Info("AnalyzeCPV1", "++++ Event %d : total %d impacts, %d charged impacts and %d  rec. points.",
 256           ievent, nTotalGen, nChargedGen, please->CpvRecPoints()->GetEntries()) ;
 257   }
 258   // Save histograms
 259
 260   Text_t outputname[80] ;
 261   sprintf(outputname,"%s.analyzed",GetFileName().Data());
 262   TFile output(outputname,"RECREATE");
 263   output.cd();
 264
 265   // Plot histograms
 266
 267   TCanvas *cpvCanvas = new TCanvas("Cpv1","CPV analysis-I",20,20,800,600);
 268   gStyle->SetOptStat(111111);
 269   gStyle->SetOptFit(1);
 270   gStyle->SetOptDate(1);
 271   cpvCanvas->Divide(3,3);
 272
 273   cpvCanvas->cd(1);
 274   gPad->SetFillColor(10);
 275   hNrp->SetFillColor(16);
 276   hNrp->Draw();
 277
 278   cpvCanvas->cd(2);
 279   gPad->SetFillColor(10);
 280   hNrpX->SetFillColor(16);
 281   hNrpX->Draw();
 282
 283   cpvCanvas->cd(3);
 284   gPad->SetFillColor(10);
 285   hNrpZ->SetFillColor(16);
 286   hNrpZ->Draw();
 287
 288   cpvCanvas->cd(4);
 289   gPad->SetFillColor(10);
 290   hDx->SetFillColor(16);
 291   hDx->Fit("gaus");
 292   hDx->Draw();
 293
 294   cpvCanvas->cd(5);
 295   gPad->SetFillColor(10);
 296   hDz->SetFillColor(16);
 297   hDz->Fit("gaus");
 298   hDz->Draw();
 299
 300   cpvCanvas->cd(6);
 301   gPad->SetFillColor(10);
 302   hChi2->SetFillColor(16);
 303   hChi2->Draw();
 304
 305   cpvCanvas->cd(7);
 306   gPad->SetFillColor(10);
 307   hEg->SetFillColor(16);
 308   hEg->Draw();
 309
 310   cpvCanvas->cd(8);
 311   gPad->SetFillColor(10);
 312   hEr->SetFillColor(16);
 313   hEr->Draw();
 314
 315   cpvCanvas->Write(0,kOverwrite);
 316
 317 }
 318
 319
 320 void AliPHOSIhepAnalyze::AnalyzeEMC1(Int_t Nevents)
 321 {
 322   //
 323   // Analyzes Emc characteristics: resolutions, cluster multiplicity,
 324   // cluster lengths along Z and Phi.
 325   // Author: Boris Polichtchouk, 24.08.2001
 326   //
 327
 328   // Book histograms
 329
 330   TH1F *hDx   = new TH1F("hDx"  ,"EMC x-resolution@reconstruction",100,-5. , 5.);
 331   TH1F *hDz   = new TH1F("hDz"  ,"EMC z-resolution@reconstruction",100,-5. , 5.);
 332   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 3.);
 333   TH1S *hNrp  = new TH1S("hNrp" ,"EMC rec.point multiplicity",      21,-0.5,20.5);
 334   TH1S *hNrpX = new TH1S("hNrpX","EMC rec.point Phi-length"  ,      21,-0.5,20.5);
 335   TH1S *hNrpZ = new TH1S("hNrpZ","EMC rec.point Z-length"    ,      21,-0.5,20.5);
 336
 337   TList * fEmcImpacts ;
 338   TBranch * branchEMCimpacts;
 339
 340   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 341   if ( please == 0 )
 342    {
 343      Error("AnalyzeEMC1","Could not obtain the Loader object !");
 344      return ;
 345    }
 346
 347   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 348
 349   Info("AnalyzeCPV1", "Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths");
 350   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 351
 352     Int_t nTotalGen = 0;
 353
 354     Int_t ntracks = gAlice->GetEvent(ievent);
 355
 356     Info("AnalyzeCPV1", " >>>>>>>Event %d .<<<<<<<", ievent) ;
 357
 358     TTree* treeH = please->TreeH();
 359     if (treeH == 0x0)
 360      {
 361        Error("AnalyzeEMC1","Can not get TreeH");
 362        return;
 363      }
 364
 365
 366     // Get branch of EMC impacts
 367     if (! (branchEMCimpacts =treeH->GetBranch("PHOSEmcImpacts")) ) {
 368       Info("AnalyzeCPV1", " Couldn't find branch PHOSEmcImpacts. Exit.");
 369       return;
 370     }
 371
 372     // Create and fill arrays of hits for each EMC module
 373
 374     Int_t nOfModules = fGeom->GetNModules();
 375     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 376     Int_t iModule = 0;
 377     for (iModule=0; iModule < nOfModules; iModule++)
 378       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 379
 380     TClonesArray    *impacts;
 381     AliPHOSImpact   *impact;
 382     TLorentzVector   p;
 383
 384     // First go through all primary tracks and fill the arrays
 385     // of hits per each EMC module
 386
 387     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 388       branchEMCimpacts ->SetAddress(&fEmcImpacts);
 389       branchEMCimpacts ->GetEntry(itrack,0);
 390
 391       for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 392         impacts = (TClonesArray *)fEmcImpacts->At(iModule);
 393         // Do loop over impacts in the module
 394         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 395           impact=(AliPHOSImpact*)impacts->At(iImpact);
 396           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 397           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 398         }
 399       }
 400       fEmcImpacts->Clear();
 401     }
 402     for (iModule=0; iModule < nOfModules; iModule++) {
 403       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 404       printf("EMC module %d has %d hits\n",iModule,nsum);
 405       nTotalGen += nsum;
 406     }
 407
 408     // Then go through reconstructed points and for each find
 409     // the closeset hit
 410     // The distance from the rec.point to the closest hit
 411     // gives the coordinate resolution of the EMC
 412
 413     fRunLoader->GetEvent(ievent);
 414     TIter nextRP(please->EmcRecPoints()) ;
 415     AliPHOSEmcRecPoint *emcRecPoint ;
 416     Float_t xgen, ygen, zgen;
 417     while( ( emcRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
 418
 419       Float_t chi2dof = ((AliPHOSEvalRecPoint*)emcRecPoint)->Chi2Dof();
 420       hChi2->Fill(chi2dof);
 421
 422       TVector3  locpos;
 423       emcRecPoint->GetLocalPosition(locpos);
 424       Int_t phosModule = emcRecPoint->GetPHOSMod();
 425       Int_t rpMult     = emcRecPoint->GetMultiplicity();
 426       Int_t rpMultX, rpMultZ;
 427       ((AliPHOSEvalRecPoint*)emcRecPoint)->GetClusterLengths(rpMultX,rpMultZ);
 428       Float_t xrec  = locpos.X();
 429       Float_t zrec  = locpos.Z();
 430       Float_t dxmin = 1.e+10;
 431       Float_t dzmin = 1.e+10;
 432       Float_t r2min = 1.e+10;
 433       Float_t r2;
 434
 435       Int_t nEMChits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
 436       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
 437       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
 438       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
 439         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
 440         xgen   = impact->X();
 441         zgen   = impact->Z();
 442         ygen   = impact->Y();
 443
 444
 445         //Transform to the local ref.frame
 446         const AliPHOSGeometry* geom = please->PHOSGeometry();
 447         Float_t phig = geom->GetPHOSAngle(phosModule);
 448         Float_t phi = TMath::Pi()/180*phig;
 449         Float_t distanceIPtoEMC = geom->GetIPtoCrystalSurface();
 450         Float_t xoL,yoL,zoL ;
 451 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
 452 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoEMC;
 453         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
 454         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoEMC;
 455         zoL = zgen;
 456
 457         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
 458         if ( r2 < r2min ) {
 459           r2min = r2;
 460           dxmin = xoL - xrec;
 461           dzmin = zoL - zrec;
 462           locImpX = xoL;
 463           locImpZ = zoL;
 464           gImpX = xgen;
 465           gImpZ = zgen;
 466           gImpY = ygen;
 467         }
 468       }
 469       Info("AnalyzeCPV1", " Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY);
 470       Info("AnalyzeCPV1", " Impact local (X,Z) = %f %f", locImpX, locImpZ);
 471       Info("AnalyzeCPV1", " Reconstructed (X,Z) = %f %f", xrec, zrec);
 472       Info("AnalyzeCPV1", " dxmin %f dzmin %f", dxmin, dzmin) ;
 473       hDx  ->Fill(dxmin);
 474       hDz  ->Fill(dzmin);
 475 //        hDr  ->Fill(TMath::Sqrt(r2min));
 476       hNrp ->Fill(rpMult);
 477       hNrpX->Fill(rpMultX);
 478       hNrpZ->Fill(rpMultZ);
 479     }
 480     delete [] hitsPerModule;
 481
 482     Info("AnalyzeCPV1", "++++ Event %d : total  %d impacts,  %d Emc rec. points.",
 483          ievent, nTotalGen, please->EmcRecPoints()->GetEntriesFast()) ;
 484
 485   }
 486   // Save histograms
 487
 488   Text_t outputname[80] ;
 489   sprintf(outputname,"%s.analyzed",GetFileName().Data());
 490   TFile output(outputname,"update");
 491   output.cd();
 492
 493   // Plot histograms
 494
 495   TCanvas *emcCanvas = new TCanvas("Emc1","EMC analysis-I",20,20,800,400);
 496   gStyle->SetOptStat(111111);
 497   gStyle->SetOptFit(1);
 498   gStyle->SetOptDate(1);
 499   emcCanvas->Divide(3,2);
 500
 501   emcCanvas->cd(1);
 502   gPad->SetFillColor(10);
 503   hNrp->SetFillColor(16);
 504   hNrp->Draw();
 505
 506   emcCanvas->cd(2);
 507   gPad->SetFillColor(10);
 508   hNrpX->SetFillColor(16);
 509   hNrpX->Draw();
 510
 511   emcCanvas->cd(3);
 512   gPad->SetFillColor(10);
 513   hNrpZ->SetFillColor(16);
 514   hNrpZ->Draw();
 515
 516   emcCanvas->cd(4);
 517   gPad->SetFillColor(10);
 518   hDx->SetFillColor(16);
 519   hDx->Fit("gaus");
 520   hDx->Draw();
 521
 522   emcCanvas->cd(5);
 523   gPad->SetFillColor(10);
 524   hDz->SetFillColor(16);
 525   hDz->Fit("gaus");
 526   hDz->Draw();
 527
 528   emcCanvas->cd(6);
 529   gPad->SetFillColor(10);
 530   hChi2->SetFillColor(16);
 531   hChi2->Draw();
 532
 533   emcCanvas->Write(0,kOverwrite);
 534 }
 535
 536 //____________________________________________________________________________
 537 void AliPHOSIhepAnalyze::AnalyzeCPV2(Int_t Nevents)
 538 {
 539   // CPV analysis - part II.
 540   // Ratio of combinatoric distances between generated
 541   // and reconstructed hits.
 542   // Author: Boris Polichtchouk (polishchuk@mx.ihep.su)
 543   // 24 March 2001
 544
 545
 546   TH1F* hDrijCPVr = new TH1F("Drij_cpv_r","Distance between reconstructed hits in CPV",140,0,50);
 547   TH1F* hDrijCPVg = new TH1F("Drij_cpv_g","Distance between generated hits in CPV",140,0,50);
 548   TH1F* hDrijCPVratio = new TH1F("Drij_cpv_ratio","R_{ij}^{rec}/R_{ij}^{gen} in CPV",140,0,50);
 549
 550 //    TH1F* hT0 = new TH1F("hT0","Type of entering particle",20000,-10000,10000);
 551
 552   hDrijCPVr->Sumw2();
 553   hDrijCPVg->Sumw2();
 554   hDrijCPVratio->Sumw2(); //correct treatment of errors
 555
 556   TList * fCpvImpacts = new TList();
 557   TBranch * branchCPVimpacts;
 558
 559   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 560   if ( please == 0 )
 561    {
 562      Error("AnalyzeCPV2","Could not obtain the Loader object !");
 563      return ;
 564    }
 565   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 566   fRunLoader->LoadHeader();
 567
 568   for (Int_t nev=0; nev<Nevents; nev++)
 569     {
 570       printf("\n=================== Event %10d ===================\n",nev);
 571       fRunLoader->GetEvent(nev);
 572       Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
 573
 574       Int_t nRecCPV = 0; // Reconstructed points in event
 575       Int_t nGenCPV = 0; // Impacts in event
 576
 577       // Get branch of CPV impacts
 578       TTree* treeH = please->TreeH();
 579       if (treeH == 0x0)
 580        {
 581         Error("AnalyzeCPV2","Can not get TreeH");
 582         return;
 583        }
 584
 585       if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) )  return;
 586
 587       // Create and fill arrays of hits for each CPV module
 588       Int_t nOfModules = fGeom->GetNModules();
 589       TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 590       Int_t iModule = 0;
 591       for (iModule=0; iModule < nOfModules; iModule++)
 592         hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 593
 594       TClonesArray    *impacts;
 595       AliPHOSImpact   *impact;
 596
 597       for (Int_t itrack=0; itrack<ntracks; itrack++) {
 598         branchCPVimpacts ->SetAddress(&fCpvImpacts);
 599         Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
 600         branchCPVimpacts ->GetEntry(itrack,0);
 601
 602         for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 603           impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 604           // Do loop over impacts in the module
 605           for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 606             impact=(AliPHOSImpact*)impacts->At(iImpact);
 607             TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 608             if(IsCharged(impact->GetPid()))
 609               new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 610           }
 611         }
 612         fCpvImpacts->Clear();
 613       }
 614
 615       for (iModule=0; iModule < nOfModules; iModule++) {
 616         Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 617         printf("CPV module %d has %d hits\n",iModule,nsum);
 618
 619         AliPHOSImpact* genHit1;
 620         AliPHOSImpact* genHit2;
 621         Int_t irp1,irp2;
 622         for(irp1=0; irp1< nsum; irp1++)
 623           {
 624             genHit1 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp1));
 625             for(irp2 = irp1+1; irp2<nsum; irp2++)
 626               {
 627                 genHit2 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp2));
 628                 Float_t dx = genHit1->X() - genHit2->X();
 629                 Float_t dz = genHit1->Z() - genHit2->Z();
 630                 Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
 631                 hDrijCPVg->Fill(dr);
 632 //                      Info("AnalyzeCPV1", "(dx dz dr): %f %f", dx, dz);
 633               }
 634           }
 635       }
 636
 637
 638   //--------- Combinatoric distance between rec. hits in CPV
 639
 640       TObjArray* cpvRecPoints = please->CpvRecPoints();
 641       nRecCPV =  cpvRecPoints->GetEntriesFast();
 642
 643       if(nRecCPV)
 644         {
 645           AliPHOSCpvRecPoint* recHit1;
 646           AliPHOSCpvRecPoint* recHit2;
 647           TIter nextCPVrec1(cpvRecPoints);
 648           while(TObject* obj1 = nextCPVrec1() )
 649             {
 650               TIter nextCPVrec2(cpvRecPoints);
 651               while (TObject* obj2 = nextCPVrec2())
 652                 {
 653                   if(!obj2->IsEqual(obj1))
 654                     {
 655                       recHit1 = (AliPHOSCpvRecPoint*)obj1;
 656                       recHit2 = (AliPHOSCpvRecPoint*)obj2;
 657                       TVector3 locpos1;
 658                       TVector3 locpos2;
 659                       recHit1->GetLocalPosition(locpos1);
 660                       recHit2->GetLocalPosition(locpos2);
 661                       Float_t dx = locpos1.X() - locpos2.X();
 662                       Float_t dz = locpos1.Z() - locpos2.Z();
 663                       Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
 664                       if(recHit1->GetPHOSMod() == recHit2->GetPHOSMod())
 665                         hDrijCPVr->Fill(dr);
 666                     }
 667                 }
 668             }
 669         }
 670
 671       Info("AnalyzeCPV1", " Event %d . Total of %d hits, %d rec.points.",
 672            nev, nGenCPV, nRecCPV) ;
 673
 674       delete [] hitsPerModule;
 675
 676     } // End of loop over events.
 677
 678
 679 //    hDrijCPVg->Draw();
 680 //    hDrijCPVr->Draw();
 681   hDrijCPVratio->Divide(hDrijCPVr,hDrijCPVg);
 682   hDrijCPVratio->Draw();
 683
 684 //    hT0->Draw();
 685
 686 }
 687
 688
 689 void AliPHOSIhepAnalyze::CpvSingle(Int_t nevents)
 690 {
 691   // Distributions of coordinates and cluster lengths of rec. points
 692   // in the case of single initial particle.
 693
 694   TH1F* hZr = new TH1F("Zrec","Reconstructed Z distribution",150,-5,5);
 695   TH1F* hXr = new TH1F("Xrec","Reconstructed X distribution",150,-14,-2);
 696   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",100, 0. , 20.);
 697
 698   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",21,-0.5,20.5);
 699   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,21,-0.5,20.5);
 700   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,21,-0.5,20.5);
 701
 702   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 703   if ( gime == 0 )
 704    {
 705      Error("CpvSingle","Could not obtain the Loader object !");
 706      return ;
 707    }
 708
 709   for(Int_t ievent=0; ievent<nevents; ievent++)
 710     {
 711       fRunLoader->GetEvent(ievent);
 712       if(gime->CpvRecPoints()->GetEntriesFast()>1) continue;
 713
 714       AliPHOSCpvRecPoint* pt = (AliPHOSCpvRecPoint*)(gime->CpvRecPoints())->At(0);
 715       if(pt) {
 716         TVector3 lpos;
 717         pt->GetLocalPosition(lpos);
 718         hXr->Fill(lpos.X());
 719         hZr->Fill(lpos.Z());
 720
 721         Int_t rpMult = pt->GetMultiplicity();
 722         hNrp->Fill(rpMult);
 723         Int_t rpMultX, rpMultZ;
 724         pt->GetClusterLengths(rpMultX,rpMultZ);
 725         hNrpX->Fill(rpMultX);
 726         hNrpZ->Fill(rpMultZ);
 727         hChi2->Fill(((AliPHOSEvalRecPoint*)pt)->Chi2Dof());
 728         Info("AnalyzeCPV1", "+++++ Event %d . (Mult,MultX,MultZ) = %d %d %d +++++",
 729              ievent, rpMult, rpMultX, rpMultZ) ;
 730
 731       }
 732
 733     }
 734
 735   Text_t outputname[80] ;
 736   sprintf(outputname,"%s.analyzed.single",GetFileName().Data());
 737   TFile output(outputname,"RECREATE");
 738   output.cd();
 739
 740   // Plot histograms
 741   TCanvas *cpvCanvas = new TCanvas("SingleParticle","Single particle events",20,20,800,400);
 742   gStyle->SetOptStat(111111);
 743   gStyle->SetOptFit(1);
 744   gStyle->SetOptDate(1);
 745   cpvCanvas->Divide(3,2);
 746
 747   cpvCanvas->cd(1);
 748   gPad->SetFillColor(10);
 749   hXr->SetFillColor(16);
 750   hXr->Draw();
 751
 752   cpvCanvas->cd(2);
 753   gPad->SetFillColor(10);
 754   hZr->SetFillColor(16);
 755   hZr->Draw();
 756
 757   cpvCanvas->cd(3);
 758   gPad->SetFillColor(10);
 759   hChi2->SetFillColor(16);
 760   hChi2->Draw();
 761
 762   cpvCanvas->cd(4);
 763   gPad->SetFillColor(10);
 764   hNrp->SetFillColor(16);
 765   hNrp->Draw();
 766
 767   cpvCanvas->cd(5);
 768   gPad->SetFillColor(10);
 769   hNrpX->SetFillColor(16);
 770   hNrpX->Draw();
 771
 772   cpvCanvas->cd(6);
 773   gPad->SetFillColor(10);
 774   hNrpZ->SetFillColor(16);
 775   hNrpZ->Draw();
 776
 777   cpvCanvas->Write(0,kOverwrite);
 778
 779 }
 780
 781 void AliPHOSIhepAnalyze::HitsCPV(Int_t nev)
 782 {
 783   // Cumulative list of charged CPV impacts in event nev.
 784
 785   TList * fCpvImpacts ;
 786   TBranch * branchCPVimpacts;
 787
 788   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 789   if ( please == 0 )
 790    {
 791      Error("HitsCPV","Could not obtain the Loader object !");
 792      return ;
 793    }
 794   const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 795
 796
 797   printf("\n=================== Event %10d ===================\n",nev);
 798   fRunLoader->GetEvent(nev);
 799   Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
 800
 801 //    Int_t nRecCPV = 0; // Reconstructed points in event // 01.10.2001
 802 //    Int_t nGenCPV = 0; // Impacts in event
 803
 804   // Get branch of CPV impacts
 805    TTree* treeH = please->TreeH();
 806    if (treeH == 0x0)
 807     {
 808       Error("CPVSingle","Can not get TreeH");
 809       return;
 810     }
 811
 812   if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) )  return;
 813
 814   // Create and fill arrays of hits for each CPV module
 815   Int_t nOfModules = fGeom->GetNModules();
 816   TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 817   Int_t iModule = 0;
 818   for (iModule=0; iModule < nOfModules; iModule++)
 819     hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 820
 821   TClonesArray    *impacts;
 822   AliPHOSImpact   *impact;
 823
 824   for (Int_t itrack=0; itrack<ntracks; itrack++) {
 825     branchCPVimpacts ->SetAddress(&fCpvImpacts);
 826     Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
 827     branchCPVimpacts ->GetEntry(itrack,0);
 828
 829     for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 830       impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 831       // Do loop over impacts in the module
 832       for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 833         impact=(AliPHOSImpact*)impacts->At(iImpact);
 834         TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 835         new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 836       }
 837     }
 838     fCpvImpacts->Clear();
 839   }
 840
 841   for (iModule=0; iModule < nOfModules; iModule++) {
 842     Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 843     printf("CPV module %d has %d hits\n",iModule,nsum);
 844   }
 845
 846 //    TList * fCpvImpacts ;
 847 //    TBranch * branchCPVimpacts;
 848 //    AliPHOSImpact* impact;
 849 //    TClonesArray* impacts;
 850
 851 //    AliPHOSLoader * please = AliPHOSLoader::GetInstance(GetFileName().Data(),"PHOS");
 852 //    const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 853
 854 //    Int_t ntracks = gAlice->GetEvent(ievent);
 855 //    Int_t nOfModules = fGeom->GetNModules();
 856 //    Info("AnalyzeCPV1", " Tracks: "<<ntracks<<"  Modules: "<<nOfModules);
 857
 858 //    if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) )  return;
 859
 860 //    for (Int_t itrack=0; itrack<ntracks; itrack++) {
 861 //      branchCPVimpacts ->SetAddress(&fCpvImpacts);
 862 //      Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
 863 //      branchCPVimpacts ->GetEntry(itrack,0);
 864 //      Info("AnalyzeCPV1", " branchCPVimpacts ->GetEntry(itrack,0) OK.");
 865
 866 //      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 867 //        impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 868 //        Info("AnalyzeCPV1", " fCpvImpacts->At(iModule) OK.");
 869 //        // Do loop over impacts in the module
 870 //        for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 871 //      impact=(AliPHOSImpact*)impacts->At(iImpact);
 872 //      impact->Print();
 873 //      if(IsCharged(impact->GetPid()))
 874 //        {
 875 //          Info("AnalyzeCPV1", " Add charged hit..";
 876 //          new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
 877 //          Info("AnalyzeCPV1", "done.");
 878 //        }
 879 //        }
 880 //      }
 881 //      fCpvImpacts->Clear();
 882 //    }
 883
 884 //    Info("AnalyzeCPV1", " PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits.");
 885
 886 }
 887
 888
 889 //  void AliPHOSIhepAnalyze::ChargedHitsCPV(TClonesArray* hits, Int_t ievent, Int_t iModule)
 890 //  {
 891 //    // Cumulative list of charged CPV hits in event ievent
 892 //    // in PHOS module iModule.
 893
 894 //    HitsCPV(hits,ievent,iModule);
 895 //    TIter next(hits);
 896
 897 //    while(AliPHOSCPVHit* cpvhit = (AliPHOSCPVHit*)next())
 898 //      {
 899 //        if(!IsCharged(cpvhit->GetIpart()))
 900 //      {
 901 //        hits->Remove(cpvhit);
 902 //        delete cpvhit;
 903 //        hits->Compress();
 904 //      }
 905 //      }
 906
 907 //    Info("AnalyzeCPV1", " PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits.");
 908 //  }
 909
 910 Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdgCode)
 911 {
 912   // For HIJING
 913   Info("AnalyzeCPV1", "pdgCode %d", pdgCode);
 914   if(pdgCode==211 || pdgCode==-211 || pdgCode==321 || pdgCode==-321 || pdgCode==11 || pdgCode==-11 || pdgCode==2212 || pdgCode==-2212) return kTRUE;
 915   else
 916     return kFALSE;
 917 }
 918 //---------------------------------------------------------------------------
 919
 920
 921
 922
 923
 924