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