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