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