PWG4/AliAnaGamma.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 /* $Id$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.4  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  * Revision 1.2  2007/08/17 12:40:04  schutz
  24  * New analysis classes by Gustavo Conesa
  25  *
  26  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  27  * new analysis classes in the the new analysis framework
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Base class for gamma and correlation analysis
  34 // It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
  35 // got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes
  36 // AliAnaGammaDirect, AliAnaGammaCorrelation ....
  37 //
  38 //*-- Author: Gustavo Conesa (INFN-LNF)
  39
  40 // --- ROOT system ---
  41
  42 #include <TParticle.h>
  43 #include <TH2.h>
  44
  45 //---- AliRoot system ----
  46 #include "AliAnaGamma.h"
  47 #include "AliGammaReader.h"
  48 #include "AliAnaGammaDirect.h"
  49 #include "AliAnaGammaCorrelation.h"
  50 #include "AliAnaGammaSelection.h"
  51 #include "AliNeutralMesonSelection.h"
  52 #include "AliAODCaloCluster.h"
  53 #include "AliAODEvent.h"
  54 #include "Riostream.h"
  55 #include "AliLog.h"
  56
  57 ClassImp(AliAnaGamma)
  58
  59
  60 //____________________________________________________________________________
  61   AliAnaGamma::AliAnaGamma() :
  62     TObject(),
  63     fOutputContainer(0x0),
  64     fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0),
  65     fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
  66     fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
  67 {
  68   //Default Ctor
  69
  70   //Initialize parameters, pointers and histograms
  71   if(!fReader)
  72     fReader = new AliGammaReader();
  73   if(!fGammaDirect)
  74     fGammaDirect = new AliAnaGammaDirect();
  75   if(!fGammaCorrelation)
  76     fGammaCorrelation = new AliAnaGammaCorrelation();
  77   if(!fGammaSelection)
  78     fGammaSelection = new AliAnaGammaSelection();
  79   if(!fNeutralMesonSelection)
  80     fNeutralMesonSelection = new AliNeutralMesonSelection();
  81
  82   fAODclusters = 0;
  83
  84   InitParameters();
  85
  86 }
  87
  88 //____________________________________________________________________________
  89 AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :
  90   TObject(),
  91   fOutputContainer(g. fOutputContainer),
  92   fAnaType(g.fAnaType),  fCalorimeter(g.fCalorimeter),
  93   fData(g.fData), fKine(g.fKine),fReader(g.fReader),
  94   fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
  95   fGammaSelection(g.fGammaSelection),
  96   fNeutralMesonSelection(g.fNeutralMesonSelection),
  97   fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
  98 {
  99   // cpy ctor
 100
 101 }
 102
 103 //_________________________________________________________________________
 104 AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
 105 {
 106   // assignment operator
 107
 108   if(this == &source)return *this;
 109   ((TObject *)this)->operator=(source);
 110
 111   fOutputContainer = source.fOutputContainer ;
 112   fAnaType = source.fAnaType;
 113   fCalorimeter = source.fCalorimeter ;
 114   fData = source.fData ;
 115   fKine = source.fKine ;
 116   fReader = source.fReader ;
 117   fGammaDirect = source.fGammaDirect ;
 118   fGammaCorrelation = source.fGammaCorrelation ;
 119   fGammaSelection = source.fGammaSelection ;
 120   fNeutralMesonSelection = source.fNeutralMesonSelection ;
 121
 122   return *this;
 123
 124 }
 125
 126 //____________________________________________________________________________
 127 AliAnaGamma::~AliAnaGamma()
 128 {
 129   // Remove all pointers.
 130
 131   // Protection added in case of NULL pointers (MG)
 132   if (fOutputContainer) {
 133      fOutputContainer->Clear();
 134      delete fOutputContainer ;
 135   }
 136
 137   if (fData) delete fData ;
 138   if (fKine) delete fKine ;
 139   if (fReader) delete fReader ;
 140   if (fGammaDirect) delete fGammaDirect ;
 141   if (fGammaCorrelation) delete fGammaCorrelation ;
 142   if (fGammaSelection) delete fGammaSelection ;
 143   if (fNeutralMesonSelection) delete fNeutralMesonSelection ;
 144
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAnaGamma::Init()
 149 {
 150
 151   //Init container histograms and other common variables
 152
 153   //Histograms container
 154   fOutputContainer = new TList ;
 155
 156   //Fill container with appropriate histograms
 157
 158   //Set selection  analysis histograms
 159   TList * selectcontainer =  fGammaSelection->GetCreateOutputObjects();
 160   for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
 161     Bool_t  add = kTRUE ;
 162     TString name = (selectcontainer->At(i))->GetName();
 163     if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
 164     if(!fReader->IsPHOSOn() && name.Contains("PHOS"))   add = kFALSE;
 165     if(!fReader->IsCTSOn() &&  !fGammaSelection->FillCTS() && name.Contains("CTS"))   add = kFALSE;
 166     if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
 167   }  //Set selection  analysis histograms
 168
 169
 170   //Set prompt photon analysis histograms
 171   TList * promptcontainer =  fGammaDirect->GetCreateOutputObjects();
 172   for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
 173     fOutputContainer->Add(promptcontainer->At(i)) ;
 174
 175   //Check if selected options are correct or set them when necessary
 176   if(fReader->GetDataType() == AliGammaReader::kMCData){
 177     fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
 178     fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
 179   }
 180
 181   if(fAnaType == kCorrelation){
 182
 183     //Check if selected options are correct
 184     if (fGammaDirect->GetICMethod()==AliAnaGammaDirect::kSeveralIC)
 185       AliFatal("Correlation not allowed with multiple isolation cuts, kCorrelation and kSeveralIC do not go together");
 186
 187     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kParton &&
 188        fReader->GetDataType() != AliGammaReader::kMC)
 189       AliFatal("kParton must be analyzed with data kMC");
 190
 191     //Set the parameters for the neutral pair selection depending on the analysis,
 192     fNeutralMesonSelection->SetDeltaPhiCutRange(fGammaCorrelation->GetDeltaPhiMinCut(),
 193                                                 fGammaCorrelation->GetDeltaPhiMaxCut());
 194
 195
 196     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kHadron){
 197       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiMinPt);
 198       fNeutralMesonSelection->SetMinPt(fGammaCorrelation->GetMinPtHadron());
 199
 200     }
 201
 202     if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kJetLeadCone){
 203       fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiPtRatio);
 204       fNeutralMesonSelection->SetRatioCutRange(fGammaCorrelation->GetRatioMinCut(),
 205                                                fGammaCorrelation->GetRatioMaxCut());
 206     }
 207
 208     //Set the neutral mesosn selection histograms
 209     TList * neutralmesoncontainer =  fNeutralMesonSelection->GetCreateOutputObjects();
 210     if(fNeutralMesonSelection->AreNeutralMesonSelectionHistosKept()){
 211       for(Int_t i = 0; i < neutralmesoncontainer->GetEntries(); i++)
 212         fOutputContainer->Add(neutralmesoncontainer->At(i)) ;
 213     }
 214
 215     //Set correlation histograms
 216     TList * correlationcontainer =  fGammaCorrelation->GetCreateOutputObjects();
 217     for(Int_t i = 0; i < correlationcontainer->GetEntries(); i++)
 218       fOutputContainer->Add(correlationcontainer->At(i)) ;
 219     fGammaCorrelation->SetOutputContainer(fOutputContainer);
 220     fGammaCorrelation->SetNeutralMesonSelection(fNeutralMesonSelection);
 221
 222   }//kCorrelation
 223
 224 }
 225
 226 //____________________________________________________________________________
 227 void AliAnaGamma::InitParameters()
 228 {
 229
 230   //Init data members
 231
 232   fAnaType = kPrompt;
 233   fCalorimeter = "EMCAL";
 234
 235 }
 236
 237 //__________________________________________________________________
 238 void AliAnaGamma::Print(const Option_t * opt) const
 239 {
 240
 241   //Print some relevant parameters set for the analysis
 242   if(! opt)
 243     return;
 244
 245   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 246   printf("Analysis type           =     %d\n", fAnaType) ;
 247   printf("Calorimeter           =     %s\n", fCalorimeter.Data()) ;
 248
 249   switch(fAnaType)
 250     {
 251     case kPrompt:
 252       {
 253         fGammaDirect->Print("");
 254       }// case kIsolationCut
 255       break;
 256
 257     case kCorrelation:
 258       {
 259         fGammaCorrelation->Print("");
 260       }//  case kCorrelation
 261       break;
 262
 263     }//switch
 264 }
 265
 266
 267 //____________________________________________________________________________
 268 Bool_t AliAnaGamma::ProcessEvent(){
 269   //Process analysis for this event
 270
 271   if(!fOutputContainer)
 272     AliFatal("Histograms not initialized");
 273
 274   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
 275
 276   TClonesArray * plCTS      = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
 277   TClonesArray * plEMCAL    = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
 278   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured  Gamma calorimeter
 279   TClonesArray * plPrimCTS      = new TClonesArray("TParticle",1000); // primary tracks
 280   TClonesArray * plPrimEMCAL    = new TClonesArray("TParticle",1000);   // primary emcal clusters
 281   TClonesArray * plPrimPHOS     = new TClonesArray("TParticle",1000);  // primary phos clusters
 282   TClonesArray * plParton   = new TClonesArray("TParticle",1000);  // All partons
 283   //Fill lists with photons, neutral particles and charged particles
 284   //look for the highest energy photon in the event inside fCalorimeter
 285
 286   //Fill particle lists
 287   if(fReader->GetDataType() == AliGammaReader::kData){
 288     AliDebug(1,"Data analysis");
 289     fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL);
 290   }
 291   else if( fReader->GetDataType()== AliGammaReader::kMC){
 292     AliDebug(1,"Kinematics analysis");
 293     fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL);
 294   }
 295   else if(fReader->GetDataType() == AliGammaReader::kMCData) {
 296    AliDebug(1,"Data + Kinematics analysis");
 297    fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS);
 298   }
 299   else
 300     AliError("Option not implemented");
 301
 302   //Fill AOD with calorimeter
 303   //Temporal solution, just for testing
 304   //FillAODs(plPHOS,plEMCAL);
 305
 306   //Select particles to do the final analysis.
 307   if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
 308   if(fReader->IsPHOSOn())  fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
 309   if(fReader->IsCTSOn() &&  fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
 310
 311
 312   //Search highest energy prompt gamma in calorimeter
 313   if(fCalorimeter == "PHOS")
 314     MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ;
 315   else if (fCalorimeter == "EMCAL")
 316     MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ;
 317   else
 318     AliFatal("Wrong calorimeter name");
 319
 320   plCTS->Clear() ;
 321   plEMCAL->Clear() ;
 322   plPHOS->Clear() ;
 323   plParton->Clear() ;
 324   plPrimCTS->Clear() ;
 325   plPrimEMCAL->Clear() ;
 326   plPrimPHOS->Clear() ;
 327
 328   delete plCTS ;
 329   delete plPHOS ;
 330   delete plEMCAL ;
 331   delete plParton ;
 332   delete plPrimCTS ;
 333   delete plPrimPHOS ;
 334   delete plPrimEMCAL ;
 335
 336   return kTRUE;
 337
 338 }
 339
 340 //____________________________________________________________________________
 341 void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo)  {
 342
 343   TParticle * pGamma = new TParticle ;
 344   Bool_t isInCalo = kFALSE ;
 345
 346   switch(fAnaType)
 347     {
 348
 349       //Only Prompt photon analysis
 350     case kPrompt:
 351       {
 352         AliDebug(1,"kPrompt analysis");
 353         switch(fGammaDirect->GetICMethod())
 354           {
 355
 356           case AliAnaGammaDirect::kSeveralIC:
 357             {
 358               fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
 359               AliDebug(1,"kSeveralIC analysis");
 360             }
 361             break;
 362
 363           default :
 364             {
 365               fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
 366               if(!isInCalo)
 367                 AliDebug(1,"Prompt gamma not found");
 368             }
 369             break;
 370           }//IC method
 371       }// case kPrompt:
 372       break;
 373
 374       //Correlate prompt photon with something: parton, hadron, jet.
 375     case kCorrelation:
 376       {
 377         AliDebug(1,"kCorrelation analysis");
 378         //Find prompt photon
 379         fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
 380
 381         if(isInCalo){//If prompt photon found, do correlation
 382
 383           switch(fGammaCorrelation->GetCorrelationType())
 384             {
 385             case AliAnaGammaCorrelation::kParton:
 386               {
 387                 AliDebug(1,"kParton correlation");
 388                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
 389               }//  case kParton
 390               break;
 391
 392             case AliAnaGammaCorrelation::kHadron:
 393               {
 394                 AliDebug(1,"kHadron correlation");
 395                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
 396               }//  case kHadron
 397               break;
 398
 399             case AliAnaGammaCorrelation::kJetLeadCone:
 400               {
 401                 AliDebug(1,"kJetLeadCone correlation");
 402                 fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
 403               }//  case kJetLeadCone
 404               break;
 405
 406             case AliAnaGammaCorrelation::kJetFinder:
 407               {
 408                 AliDebug(1,"kJetFinder correlation");
 409                 printf("Analysis not implemented \n");
 410               }//  case kJetFinder
 411               break;
 412             }// switch correlation
 413         }// is in calo
 414         else  AliDebug(2,"Prompt gamma not found");
 415       }// case kCorrelation
 416       break;
 417
 418     } //switch(fAnaType)
 419
 420   delete pGamma ;
 421
 422 }
 423
 424 //____________________________________________________
 425 void AliAnaGamma::AddCluster(AliAODCaloCluster p)
 426 {
 427 // Add new jet to the list
 428   cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
 429   cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
 430   new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
 431 }
 432
 433 //___________________________________________________
 434 void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
 435 {
 436 // Connect to the AOD
 437   fAODclusters = aod->GetCaloClusters();
 438 }
 439
 440 //____________________________________________________________________________
 441 void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
 442   //Fill AOD caloClusters
 443   //Temporal method, just for testing AOD creation
 444
 445   //Fill AOD with PHOS clusters
 446   Int_t nphos =  plPHOS->GetEntries() ;
 447   cout<<"PHOS entries "<<nphos<<endl;
 448   for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
 449     TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
 450     Float_t pos []= {0,0,0};
 451     AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
 452     AddCluster(phos);
 453   }
 454
 455   //Fill AOD with EMCAL clusters
 456   cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
 457   for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
 458     TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
 459     Float_t pos []= {0,0,0};
 460     AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
 461     AddCluster(emcal);
 462   }
 463 }