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