PHOS/AliPHOSIhepAnalyze.cxx

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