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 "AliLog.h"
  47 #include "AliPHOSIhepAnalyze.h"
  48 #include "AliPHOSDigit.h"
  49 #include "AliPHOSRecParticle.h"
  50 #include "AliPHOSLoader.h"
  51 #include "AliPHOSHit.h"
  52 #include "AliPHOSImpact.h"
  53 #include "AliPHOSvImpacts.h"
  54 #include "AliPHOSCpvRecPoint.h"
  55 #include "AliRun.h"
  56 #include "AliPHOSGeometry.h"
  57 #include "AliPHOSEvalRecPoint.h"
  58
  59 ClassImp(AliPHOSIhepAnalyze)
  60
  61 //____________________________________________________________________________
  62 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze():
  63   fRunLoader(0),
  64   fFileName()
  65 {
  66 }
  67
  68 //____________________________________________________________________________
  69 AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) :
  70   fRunLoader(0),
  71   fFileName(name)
  72 {
  73   // Constructor: open a header file
  74   fRunLoader = AliRunLoader::Open(fFileName);
  75   if (fRunLoader == 0x0)
  76    {
  77      AliFatal(Form("Can not load event from file %s",name));
  78    }
  79 }
  80
  81 //____________________________________________________________________________
  82 void AliPHOSIhepAnalyze::AnalyzeCPV1(Int_t Nevents)
  83 {
  84   //
  85   // Analyzes CPV characteristics: resolutions, cluster multiplicity,
  86   // cluster lengths along Z and Phi.
  87   // Author: Yuri Kharlov
  88   // 9 October 2000
  89   // Modified by Boris Polichtchouk, 3.07.2001
  90   //
  91
  92   // Book histograms
  93
  94   TH1F *hDx   = new TH1F("hDx"  ,"CPV x-resolution@reconstruction",100,-5. , 5.);
  95   TH1F *hDz   = new TH1F("hDz"  ,"CPV z-resolution@reconstruction",100,-5. , 5.);
  96   TH1F *hChi2 = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 10.);
  97   TH1S *hNrp  = new TH1S("hNrp" ,"CPV rec.point multiplicity",      21,-0.5,20.5);
  98   TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length"  ,      21,-0.5,20.5);
  99   TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length"    ,      21,-0.5,20.5);
 100
 101   TH1F *hEg   = new TH1F("hEg"  ,"Energies of impacts",30,0.,6.);
 102   TH1F *hEr   = new TH1F("hEr"  ,"Amplitudes of rec. points",50,0.,20.);
 103
 104   TList * fCpvImpacts ;
 105   TBranch * branchCPVimpacts;
 106
 107
 108
 109   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 110   if ( please == 0 )
 111    {
 112      AliError(Form("Could not obtain the Loader object !"));
 113      return ;
 114    }
 115
 116   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 117
 118   AliInfo(Form("Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths")) ;
 119   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 120
 121     Int_t nTotalGen = 0;
 122     Int_t nChargedGen = 0;
 123
 124     Int_t ntracks = gAlice->GetEvent(ievent);
 125     AliInfo(Form(">>>>>>>Event %d .<<<<<<<", ievent)) ;
 126
 127   /******************************************************************/
 128       TTree* treeH = please->TreeH();
 129       if (treeH == 0x0)
 130        {
 131         AliError(Form("Can not get TreeH"));
 132          return;
 133        }
 134 /******************************************************************/
 135
 136     // Get branch of CPV impacts
 137     if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) ) {
 138       AliWarning(Form("Couldn't find branch PHOSCpvImpacts. Exit.")) ;
 139       return;
 140     }
 141
 142     // Create and fill arrays of hits for each CPV module
 143
 144     Int_t nOfModules = phosgeom->GetNModules();
 145     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 146     Int_t iModule = 0;
 147     for (iModule=0; iModule < nOfModules; iModule++)
 148       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 149
 150     TClonesArray    *impacts;
 151     AliPHOSImpact   *impact;
 152     TLorentzVector   p;
 153
 154     // First go through all primary tracks and fill the arrays
 155     // of hits per each CPV module
 156
 157     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 158       branchCPVimpacts ->SetAddress(&fCpvImpacts);
 159       branchCPVimpacts ->GetEntry(itrack,0);
 160
 161       for (iModule=0; iModule < nOfModules; iModule++) {
 162         impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 163         // Do loop over impacts in the module
 164         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 165           impact=(AliPHOSImpact*)impacts->At(iImpact);
 166           hEg->Fill(impact->GetMomentum().E());
 167           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 168           if(IsCharged(impact->GetPid())) nChargedGen++;
 169           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 170         }
 171       }
 172       fCpvImpacts->Clear();
 173     }
 174     for (iModule=0; iModule < nOfModules; iModule++) {
 175       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 176       printf("CPV module %d has %d impacts\n",iModule,nsum);
 177       nTotalGen += nsum;
 178     }
 179
 180     // Then go through reconstructed points and for each find
 181     // the closeset hit
 182     // The distance from the rec.point to the closest hit
 183     // gives the coordinate resolution of the CPV
 184
 185     fRunLoader->GetEvent(ievent);
 186     TIter nextRP(please->CpvRecPoints()) ;
 187     AliPHOSCpvRecPoint *cpvRecPoint ;
 188     Float_t xgen, ygen, zgen;
 189     while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
 190
 191       Float_t chi2dof = ((AliPHOSEvalRecPoint*)cpvRecPoint)->Chi2Dof();
 192       hChi2->Fill(chi2dof);
 193       hEr->Fill(cpvRecPoint->GetEnergy());
 194
 195       TVector3  locpos;
 196       cpvRecPoint->GetLocalPosition(locpos);
 197       Int_t phosModule = cpvRecPoint->GetPHOSMod();
 198       Int_t rpMult     = cpvRecPoint->GetMultiplicity();
 199       Int_t rpMultX, rpMultZ;
 200       cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
 201       Float_t xrec  = locpos.X();
 202       Float_t zrec  = locpos.Z();
 203       Float_t dxmin = 1.e+10;
 204       Float_t dzmin = 1.e+10;
 205       Float_t r2min = 1.e+10;
 206       Float_t r2;
 207
 208       Int_t nCPVhits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
 209       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
 210       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
 211       for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
 212         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
 213         xgen   = impact->X();
 214         zgen   = impact->Z();
 215         ygen   = impact->Y();
 216
 217         //Transform to the local ref.frame
 218         Float_t phig = phosgeom->GetPHOSAngle(phosModule);
 219         Float_t phi = TMath::Pi()/180*phig;
 220         Float_t distanceIPtoCPV = phosgeom->GetIPtoOuterCoverDistance() -
 221           (phosgeom->GetFTPosition(1)+
 222            phosgeom->GetFTPosition(2)+
 223            phosgeom->GetCPVTextoliteThickness())/2;
 224         Float_t xoL,yoL,zoL ;
 225 //      xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
 226 //      yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoCPV;
 227         xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
 228         yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoCPV;
 229         zoL = zgen;
 230
 231         r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
 232         if ( r2 < r2min ) {
 233           r2min = r2;
 234           dxmin = xoL - xrec;
 235           dzmin = zoL - zrec;
 236           locImpX = xoL;
 237           locImpZ = zoL;
 238           gImpX = xgen;
 239           gImpZ = zgen;
 240           gImpY = ygen;
 241         }
 242       }
 243       AliInfo(Form("Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY));
 244       AliInfo(Form("Impact local (X,Z) = %f %f", locImpX, locImpZ));
 245       AliInfo(Form("Reconstructed (X,Z) = %f %f", xrec, zrec));
 246       AliInfo(Form("dxmin %f dzmin %f", dxmin, dzmin));
 247       hDx  ->Fill(dxmin);
 248       hDz  ->Fill(dzmin);
 249 //        hDr  ->Fill(TMath::Sqrt(r2min));
 250       hNrp ->Fill(rpMult);
 251       hNrpX->Fill(rpMultX);
 252       hNrpZ->Fill(rpMultZ);
 253     }
 254     delete [] hitsPerModule;
 255
 256     AliInfo(Form("++++ Event %d : total %d impacts, %d charged impacts and %d  rec. points.",
 257           ievent, nTotalGen, nChargedGen, please->CpvRecPoints()->GetEntries())) ;
 258   }
 259   // Save histograms
 260
 261   Text_t outputname[80] ;
 262   snprintf(outputname,80,"%s.analyzed",GetFileName().Data());
 263   TFile output(outputname,"RECREATE");
 264   output.cd();
 265
 266   // Plot histograms
 267
 268   TCanvas *cpvCanvas = new TCanvas("Cpv1","CPV analysis-I",20,20,800,600);
 269   gStyle->SetOptStat(111111);
 270   gStyle->SetOptFit(1);
 271   gStyle->SetOptDate(1);
 272   cpvCanvas->Divide(3,3);
 273
 274   cpvCanvas->cd(1);
 275   gPad->SetFillColor(10);
 276   hNrp->SetFillColor(16);
 277   hNrp->Draw();
 278
 279   cpvCanvas->cd(2);
 280   gPad->SetFillColor(10);
 281   hNrpX->SetFillColor(16);
 282   hNrpX->Draw();
 283
 284   cpvCanvas->cd(3);
 285   gPad->SetFillColor(10);
 286   hNrpZ->SetFillColor(16);
 287   hNrpZ->Draw();
 288
 289   cpvCanvas->cd(4);
 290   gPad->SetFillColor(10);
 291   hDx->SetFillColor(16);
 292   hDx->Fit("gaus");
 293   hDx->Draw();
 294
 295   cpvCanvas->cd(5);
 296   gPad->SetFillColor(10);
 297   hDz->SetFillColor(16);
 298   hDz->Fit("gaus");
 299   hDz->Draw();
 300
 301   cpvCanvas->cd(6);
 302   gPad->SetFillColor(10);
 303   hChi2->SetFillColor(16);
 304   hChi2->Draw();
 305
 306   cpvCanvas->cd(7);
 307   gPad->SetFillColor(10);
 308   hEg->SetFillColor(16);
 309   hEg->Draw();
 310
 311   cpvCanvas->cd(8);
 312   gPad->SetFillColor(10);
 313   hEr->SetFillColor(16);
 314   hEr->Draw();
 315
 316   cpvCanvas->Write(0,kOverwrite);
 317
 318 }
 319
 320
 321 void AliPHOSIhepAnalyze::AnalyzeEMC1(Int_t Nevents)
 322 {
 323   //
 324   // Analyzes Emc characteristics: resolutions, cluster multiplicity,
 325   // cluster lengths along Z and Phi.
 326   // Author: Boris Polichtchouk, 24.08.2001
 327   //
 328
 329   // Book histograms
 330
 331   TH1F *hDx   = new TH1F("hDx"  ,"EMC x-resolution@reconstruction",100,-5. , 5.);
 332   TH1F *hDz   = new TH1F("hDz"  ,"EMC z-resolution@reconstruction",100,-5. , 5.);
 333   TH1F *hChi2   = new TH1F("hChi2"  ,"Chi2/dof of one-gamma fit",30, 0. , 3.);
 334   TH1S *hNrp  = new TH1S("hNrp" ,"EMC rec.point multiplicity",      21,-0.5,20.5);
 335   TH1S *hNrpX = new TH1S("hNrpX","EMC rec.point Phi-length"  ,      21,-0.5,20.5);
 336   TH1S *hNrpZ = new TH1S("hNrpZ","EMC rec.point Z-length"    ,      21,-0.5,20.5);
 337
 338   TList * fEmcImpacts ;
 339   TBranch * branchEMCimpacts;
 340
 341   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 342   if ( please == 0 )
 343    {
 344      AliError(Form("Could not obtain the Loader object !"));
 345      return ;
 346    }
 347
 348   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 349
 350   AliInfo(Form("Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths"));
 351   for ( Int_t ievent=0; ievent<Nevents; ievent++) {
 352
 353     Int_t nTotalGen = 0;
 354
 355     Int_t ntracks = gAlice->GetEvent(ievent);
 356
 357     AliInfo(Form(" >>>>>>>Event %d .<<<<<<<", ievent)) ;
 358
 359     TTree* treeH = please->TreeH();
 360     if (treeH == 0x0)
 361      {
 362       AliError(Form("Can not get TreeH"));
 363        return;
 364      }
 365
 366
 367     // Get branch of EMC impacts
 368     if (! (branchEMCimpacts =treeH->GetBranch("PHOSEmcImpacts")) ) {
 369       AliWarning(Form(" Couldn't find branch PHOSEmcImpacts. Exit."));
 370       return;
 371     }
 372
 373     // Create and fill arrays of hits for each EMC module
 374
 375     Int_t nOfModules = phosgeom->GetNModules();
 376     TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 377     Int_t iModule = 0;
 378     for (iModule=0; iModule < nOfModules; iModule++)
 379       hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 380
 381     TClonesArray    *impacts;
 382     AliPHOSImpact   *impact;
 383     TLorentzVector   p;
 384
 385     // First go through all primary tracks and fill the arrays
 386     // of hits per each EMC module
 387
 388     for (Int_t itrack=0; itrack<ntracks; itrack++) {
 389       branchEMCimpacts ->SetAddress(&fEmcImpacts);
 390       branchEMCimpacts ->GetEntry(itrack,0);
 391
 392       for (iModule=0; iModule < nOfModules; iModule++) {
 393         impacts = (TClonesArray *)fEmcImpacts->At(iModule);
 394         // Do loop over impacts in the module
 395         for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 396           impact=(AliPHOSImpact*)impacts->At(iImpact);
 397           TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 398           new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 399         }
 400       }
 401       fEmcImpacts->Clear();
 402     }
 403     for (iModule=0; iModule < nOfModules; iModule++) {
 404       Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 405       printf("EMC module %d has %d hits\n",iModule,nsum);
 406       nTotalGen += nsum;
 407     }
 408
 409     // Then go through reconstructed points and for each find
 410     // the closeset hit
 411     // The distance from the rec.point to the closest hit
 412     // gives the coordinate resolution of the EMC
 413
 414     fRunLoader->GetEvent(ievent);
 415     TIter nextRP(please->EmcRecPoints()) ;
 416     AliPHOSEmcRecPoint *emcRecPoint ;
 417     Float_t xgen, ygen, zgen;
 418     while( ( emcRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
 419
 420       Float_t chi2dof = ((AliPHOSEvalRecPoint*)emcRecPoint)->Chi2Dof();
 421       hChi2->Fill(chi2dof);
 422
 423       TVector3  locpos;
 424       emcRecPoint->GetLocalPosition(locpos);
 425       Int_t phosModule = emcRecPoint->GetPHOSMod();
 426       Int_t rpMult     = emcRecPoint->GetMultiplicity();
 427       Int_t rpMultX, rpMultZ;
 428       ((AliPHOSEvalRecPoint*)emcRecPoint)->GetClusterLengths(rpMultX,rpMultZ);
 429       Float_t xrec  = locpos.X();
 430       Float_t zrec  = locpos.Z();
 431       Float_t dxmin = 1.e+10;
 432       Float_t dzmin = 1.e+10;
 433       Float_t r2min = 1.e+10;
 434       Float_t r2;
 435
 436       Int_t nEMChits  = (hitsPerModule[phosModule-1])->GetEntriesFast();
 437       Float_t locImpX=1e10,locImpZ=1e10;         // local coords of closest impact
 438       Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
 439       for (Int_t ihit=0; ihit<nEMChits; ihit++) {
 440         impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
 441         xgen   = impact->X();
 442         zgen   = impact->Z();
 443         ygen   = impact->Y();
 444
 445
 446         //Transform to the local ref.frame
 447         Float_t phig = phosgeom->GetPHOSAngle(phosModule);
 448         Float_t phi = TMath::Pi()/180*phig;
 449         Float_t distanceIPtoEMC = phosgeom->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       AliInfo(Form(" Impact global (X,Z,Y) = %f %f %f", gImpX, gImpZ, gImpY));
 470       AliInfo(Form(" Impact local (X,Z) = %f %f", locImpX, locImpZ));
 471       AliInfo(Form(" Reconstructed (X,Z) = %f %f", xrec, zrec));
 472       AliInfo(Form(" 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     AliInfo(Form("++++ 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   snprintf(outputname,80,"%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      AliError(Form("Could not obtain the Loader object !"));
 563      return ;
 564    }
 565   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 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          AliError(Form("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 = phosgeom->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         AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
 600         branchCPVimpacts ->GetEntry(itrack,0);
 601
 602         for (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 //                      AliInfo(Form("(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       AliInfo(Form(" 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      AliError(Form("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         AliInfo(Form("+++++ Event %d . (Mult,MultX,MultZ) = %d %d %d +++++",
 729              ievent, rpMult, rpMultX, rpMultZ)) ;
 730
 731       }
 732
 733     }
 734
 735   Text_t outputname[80] ;
 736   snprintf(outputname,80,"%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   AliPHOSLoader* please = dynamic_cast<AliPHOSLoader*>(fRunLoader->GetLoader("PHOSLoader"));
 786   if ( please == 0 )
 787    {
 788      AliError(Form("Could not obtain the Loader object !"));
 789      return ;
 790    }
 791
 792
 793   printf("\n=================== Event %10d ===================\n",nev);
 794   //16.03.2011: DP. Code below seems to be obsollete
 795   //Comment it to sutisfy Coverity
 796 /*
 797   TList * fCpvImpacts ;
 798   TBranch * branchCPVimpacts;
 799
 800   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 801
 802   fRunLoader->GetEvent(nev);
 803   Int_t ntracks = fRunLoader->GetHeader()->GetNtrack();
 804
 805 //    Int_t nRecCPV = 0; // Reconstructed points in event // 01.10.2001
 806 //    Int_t nGenCPV = 0; // Impacts in event
 807
 808   // Get branch of CPV impacts
 809    TTree* treeH = please->TreeH();
 810    if (treeH == 0x0)
 811     {
 812       AliError(Form("Can not get TreeH"));
 813       return;
 814     }
 815
 816   if (! (branchCPVimpacts =treeH->GetBranch("PHOSCpvImpacts")) )  return;
 817
 818   // Create and fill arrays of hits for each CPV module
 819   Int_t nOfModules = phosgeom->GetNModules();
 820   TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
 821   Int_t iModule = 0;
 822   for (iModule=0; iModule < nOfModules; iModule++)
 823     hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
 824
 825   TClonesArray    *impacts;
 826   AliPHOSImpact   *impact;
 827
 828   for (Int_t itrack=0; itrack<ntracks; itrack++) {
 829     branchCPVimpacts ->SetAddress(&fCpvImpacts);
 830     AliInfo(Form(" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."));
 831     branchCPVimpacts ->GetEntry(itrack,0);
 832
 833     for (iModule=0; iModule < nOfModules; iModule++) {
 834       impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 835       // Do loop over impacts in the module
 836       for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 837         impact=(AliPHOSImpact*)impacts->At(iImpact);
 838         TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
 839         new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
 840       }
 841     }
 842     fCpvImpacts->Clear();
 843   }
 844
 845   for (iModule=0; iModule < nOfModules; iModule++) {
 846     Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
 847     printf("CPV module %d has %d hits\n",iModule,nsum);
 848   }
 849
 850 */
 851
 852
 853 //    TList * fCpvImpacts ;
 854 //    TBranch * branchCPVimpacts;
 855 //    AliPHOSImpact* impact;
 856 //    TClonesArray* impacts;
 857
 858 //    AliPHOSLoader * please = AliPHOSLoader::GetInstance(GetFileName().Data(),"PHOS");
 859 //    const AliPHOSGeometry *  fGeom  = please->PHOSGeometry();
 860
 861 //    Int_t ntracks = gAlice->GetEvent(ievent);
 862 //    Int_t nOfModules = fGeom->GetNModules();
 863 //    AliInfo(Form(" Tracks: "<<ntracks<<"  Modules: "<<nOfModules));
 864
 865 //    if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) )  return;
 866
 867 //    for (Int_t itrack=0; itrack<ntracks; itrack++) {
 868 //      branchCPVimpacts ->SetAddress(&fCpvImpacts);
 869 //      Info("AnalyzeCPV1", " branchCPVimpacts ->SetAddress(&fCpvImpacts) OK.");
 870 //      branchCPVimpacts ->GetEntry(itrack,0);
 871 //      Info(Form(" branchCPVimpacts ->GetEntry(itrack,0) OK."));
 872
 873 //      for (Int_t iModule=0; iModule < nOfModules; iModule++) {
 874 //        impacts = (TClonesArray *)fCpvImpacts->At(iModule);
 875 //        Info(Form(" fCpvImpacts->At(iModule) OK."));
 876 //        // Do loop over impacts in the module
 877 //        for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
 878 //      impact=(AliPHOSImpact*)impacts->At(iImpact);
 879 //      impact->Print();
 880 //      if(IsCharged(impact->GetPid()))
 881 //        {
 882 //          Info(Form(" Add charged hit.."));
 883 //          new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
 884 //          Info(Form("done."));
 885 //        }
 886 //        }
 887 //      }
 888 //      fCpvImpacts->Clear();
 889 //    }
 890
 891 //    Info(Form(" PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits."));
 892
 893 }
 894
 895
 896 //  void AliPHOSIhepAnalyze::ChargedHitsCPV(TClonesArray* hits, Int_t ievent, Int_t iModule)
 897 //  {
 898 //    // Cumulative list of charged CPV hits in event ievent
 899 //    // in PHOS module iModule.
 900
 901 //    HitsCPV(hits,ievent,iModule);
 902 //    TIter next(hits);
 903
 904 //    while(AliPHOSCPVHit* cpvhit = (AliPHOSCPVHit*)next())
 905 //      {
 906 //        if(!IsCharged(cpvhit->GetIpart()))
 907 //      {
 908 //        hits->Remove(cpvhit);
 909 //        delete cpvhit;
 910 //        hits->Compress();
 911 //      }
 912 //      }
 913
 914 //    Info(Form(" PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits."));
 915 //  }
 916
 917 Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdgCode)
 918 {
 919   // For HIJING
 920   AliInfo(Form("pdgCode %d", pdgCode));
 921   if(pdgCode==211 || pdgCode==-211 || pdgCode==321 || pdgCode==-321 || pdgCode==11 || pdgCode==-11 || pdgCode==2212 || pdgCode==-2212) return kTRUE;
 922   else
 923     return kFALSE;
 924 }
 925 //---------------------------------------------------------------------------
 926
 927
 928
 929
 930
 931