PWG4/UserTasks/PHOS_PbPbQA/AliAnalysisTaskPHOSPbPbQA.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 #include "AliAnalysisManager.h"
  18 #include "AliAnalysisTaskSE.h"
  19 #include "AliAnalysisTaskPHOSPbPbQA.h"
  20 #include "AliCaloPhoton.h"
  21 #include "AliPHOSGeometry.h"
  22 #include "AliESDEvent.h"
  23 #include "AliESDCaloCells.h"
  24 #include "AliCentrality.h"
  25 #include "AliLog.h"
  26 #include "TObjArray.h"
  27 #include "TList.h"
  28 #include "TH1.h"
  29 #include "TH2.h"
  30
  31 // Stripped-down version of Dmitri Peressounko' AliAnalysisTaskPi0Flow class
  32 // used for the fast QA of PbPb data.
  33 //...
  34 // Author: Boris Polishchuk
  35 // Date   : 19.10.2011
  36
  37 ClassImp(AliAnalysisTaskPHOSPbPbQA)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskPHOSPbPbQA::AliAnalysisTaskPHOSPbPbQA() : AliAnalysisTaskSE(),
  41   fOutputContainer(0),fPHOSEvent(0),fCentrality(0),fCenBin(0),
  42   fPHOSGeo(0),fEventCounter(0)
  43 {
  44   //Default constructor
  45
  46   for(Int_t i=0;i<1;i++){
  47     for(Int_t j=0;j<2;j++)
  48       fPHOSEvents[i][j]=0 ;
  49   }
  50
  51   // Initialize the PHOS geometry
  52   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  53
  54 }
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskPHOSPbPbQA::AliAnalysisTaskPHOSPbPbQA(const char *name)
  58 : AliAnalysisTaskSE(name),
  59   fOutputContainer(0),
  60   fPHOSEvent(0),
  61   fCentrality(0),fCenBin(0),
  62   fPHOSGeo(0),
  63   fEventCounter(0)
  64 {
  65   // Constructor
  66   for(Int_t i=0;i<1;i++){
  67     for(Int_t j=0;j<2;j++)
  68         fPHOSEvents[i][j]=0 ;
  69   }
  70
  71   // Output slots #0 write into a TH1 container
  72   DefineOutput(1,TList::Class());
  73
  74   // Initialize the PHOS geometry
  75   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  76
  77 }
  78
  79 //________________________________________________________________________
  80 void AliAnalysisTaskPHOSPbPbQA::UserCreateOutputObjects()
  81 {
  82   // Create histograms
  83   // Called once
  84
  85   // ESD histograms
  86   if(fOutputContainer != NULL){
  87     delete fOutputContainer;
  88   }
  89
  90   fOutputContainer = new TList();
  91   fOutputContainer->SetOwner(kTRUE);
  92
  93
  94   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
  95   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
  96   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
  97
  98   //pi0 spectrum
  99   Int_t nPtPhot = 300 ;
 100   Double_t ptPhotMax = 30 ;
 101   Int_t nM       = 500;
 102   Double_t mMin  = 0.0;
 103   Double_t mMax  = 1.0;
 104   char key[55] ;
 105
 106   for(Int_t cent=0; cent<2; cent++){
 107
 108     snprintf(key,55,"hPi0All_cen%d",cent) ;
 109     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 110
 111     snprintf(key,55,"hMiPi0All_cen%d",cent) ;
 112     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 113   }
 114
 115   //per module
 116   for(Int_t cent=0; cent<2; cent++){
 117     for(Int_t sm=1; sm<4; sm++) {
 118
 119       snprintf(key,55,"hPi0AllSM%d_cen%d",sm,cent) ;
 120       fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 121
 122       snprintf(key,55,"hMiPi0AllSM%d_cen%d",sm,cent) ;
 123       fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 124     }
 125   }
 126
 127   PostData(1, fOutputContainer);
 128
 129 }
 130
 131 //________________________________________________________________________
 132 void AliAnalysisTaskPHOSPbPbQA::UserExec(Option_t *)
 133 {
 134   // Main loop, called for each event
 135   // Analyze ESD/AOD
 136
 137   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 138
 139   if (!event) {
 140     Printf("ERROR: Could not retrieve event");
 141     PostData(1, fOutputContainer);
 142     return;
 143   }
 144
 145   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 146   char key[55] ;
 147
 148   if(fEventCounter == 0) {
 149     for(Int_t mod=0; mod<5; mod++) {
 150       if(!event->GetPHOSMatrix(mod)) continue;
 151       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 152     }
 153     fEventCounter++ ;
 154   }
 155
 156   Int_t zvtx=0 ;
 157
 158   AliCentrality *centrality = event->GetCentrality();
 159   fCentrality=centrality->GetCentralityPercentile("V0M");
 160
 161   if( fCentrality <= 0. ){
 162     PostData(1, fOutputContainer);
 163     return;
 164   }
 165
 166   if(fCentrality < 20.) fCenBin = 0;
 167   else fCenBin = 1;
 168
 169   printf("centrality %.3f [%d]\n",fCentrality,fCenBin);
 170
 171   if(!fPHOSEvents[zvtx][fCenBin])
 172     fPHOSEvents[zvtx][fCenBin]=new TList() ;
 173
 174   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin] ;
 175
 176   if(fPHOSEvent)
 177     fPHOSEvent->Clear() ;
 178   else
 179     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
 180
 181   Int_t multClust = event->GetNumberOfCaloClusters();
 182   AliESDCaloCells * cells = event->GetPHOSCells() ;
 183
 184   FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
 185   FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
 186
 187   Int_t inPHOS = 0;
 188   Double_t vtx0[3] = {0,0,0}; // vertex
 189
 190   for (Int_t i=0; i<multClust; i++) {
 191
 192     AliESDCaloCluster *clu = event->GetCaloCluster(i);
 193
 194     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 195     if(clu->GetNCells()<3) continue;
 196
 197     Float_t  position[3];
 198     clu->GetPosition(position);
 199     TVector3 global(position) ;
 200     Int_t relId[4] ;
 201     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 202     Int_t mod  = relId[0] ;
 203     Int_t cellX = relId[2];
 204     Int_t cellZ = relId[3] ;
 205
 206     TLorentzVector pv1 ;
 207     clu->GetMomentum(pv1 ,vtx0);
 208
 209     if(inPHOS>=fPHOSEvent->GetSize()){
 210       fPHOSEvent->Expand(inPHOS+50) ;
 211     }
 212
 213     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 214     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 215     ph->SetModule(mod) ;
 216     ph->SetMomV2(&pv1) ;
 217     ph->SetNCells(clu->GetNCells());
 218
 219     ph->SetEMCx(float(cellX)) ;
 220     ph->SetEMCz(float(cellZ)) ;
 221
 222     inPHOS++ ;
 223   }
 224
 225   FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
 226
 227   //pi0
 228   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 229     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 230     Int_t sm1 = ph1->Module();
 231
 232     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 233       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 234
 235       Int_t sm2 = ph2->Module();
 236       TLorentzVector p12  = *ph1  + *ph2;
 237
 238       snprintf(key,55,"hPi0All_cen%d",fCenBin) ;
 239       FillHistogram(key,p12.M() ,p12.Pt()) ;
 240
 241       if(sm1==sm2) {
 242         snprintf(key,55,"hPi0AllSM%d_cen%d",sm1,fCenBin) ;
 243          FillHistogram(key,p12.M() ,p12.Pt()) ;
 244       }
 245
 246     } // end of loop i2
 247   } // end of loop i1
 248
 249   //now mixed
 250   for (Int_t i1=0; i1<inPHOS; i1++) {
 251     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 252     Int_t sm1 = ph1->Module();
 253
 254     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 255       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 256
 257       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 258         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 259
 260         Int_t sm2 = ph2->Module();
 261         TLorentzVector p12  = *ph1  + *ph2;
 262
 263         snprintf(key,55,"hMiPi0All_cen%d",fCenBin) ;
 264         FillHistogram(key,p12.M() ,p12.Pt()) ;
 265
 266         if(sm1==sm2) {
 267           snprintf(key,55,"hMiPi0AllSM%d_cen%d",sm1,fCenBin) ;
 268           FillHistogram(key,p12.M() ,p12.Pt()) ;
 269         }
 270
 271       } // end of loop i2
 272     }
 273   } // end of loop i1
 274
 275
 276   //Now we either add current events to stack or remove
 277   //If no photons in current event - no need to add it to mixed
 278   if(fPHOSEvent->GetEntriesFast()>0){
 279     prevPHOS->AddFirst(fPHOSEvent) ;
 280     fPHOSEvent=0;
 281     if(prevPHOS->GetSize()>100){//Remove redundant events
 282       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 283       prevPHOS->RemoveLast() ;
 284       delete tmp ;
 285     }
 286   }
 287   // Post output data.
 288   PostData(1, fOutputContainer);
 289   fEventCounter++;
 290 }
 291
 292 //_____________________________________________________________________________
 293 void AliAnalysisTaskPHOSPbPbQA::FillHistogram(const char * key,Double_t x)const{
 294   //FillHistogram
 295   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
 296   if(tmpI){
 297     tmpI->Fill(x) ;
 298     return ;
 299   }
 300   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
 301   if(tmpF){
 302     tmpF->Fill(x) ;
 303     return ;
 304   }
 305   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
 306   if(tmpD){
 307     tmpD->Fill(x) ;
 308     return ;
 309   }
 310   AliInfo(Form("can not find histogram <%s> ",key)) ;
 311 }
 312 //_____________________________________________________________________________
 313 void AliAnalysisTaskPHOSPbPbQA::FillHistogram(const char * key,Double_t x,Double_t y)const{
 314   //FillHistogram
 315   TObject * tmp = fOutputContainer->FindObject(key) ;
 316   if(!tmp){
 317     AliInfo(Form("can not find histogram <%s> ",key)) ;
 318     return ;
 319   }
 320   if(tmp->IsA() == TClass::GetClass("TH1F")){
 321     ((TH1F*)tmp)->Fill(x,y) ;
 322     return ;
 323   }
 324   if(tmp->IsA() == TClass::GetClass("TH2F")){
 325     ((TH2F*)tmp)->Fill(x,y) ;
 326     return ;
 327   }
 328   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
 329 }