PWG4/UserTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TObjArray.h"
   4 #include "TF1.h"
   5 #include "TFile.h"
   6 #include "TH1F.h"
   7 #include "TH2F.h"
   8 #include "TH2I.h"
   9 #include "TH3F.h"
  10 #include "TParticle.h"
  11 #include "TCanvas.h"
  12 #include "TStyle.h"
  13 #include "TRandom.h"
  14
  15 #include "AliAODMCParticle.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliMCEventHandler.h"
  18 #include "AliMCEvent.h"
  19 #include "AliStack.h"
  20 #include "AliAnalysisTaskSE.h"
  21 #include "AliAnalysisTaskPi0DiffEfficiency.h"
  22 #include "AliCaloPhoton.h"
  23 #include "AliPHOSGeometry.h"
  24 #include "AliPHOSEsdCluster.h"
  25 #include "AliPHOSCalibData.h"
  26 #include "AliAODEvent.h"
  27 #include "AliAODCaloCluster.h"
  28 #include "AliAODVertex.h"
  29 #include "AliESDtrackCuts.h"
  30 #include "AliLog.h"
  31 #include "AliPID.h"
  32 #include "AliCDBManager.h"
  33 #include "AliCentrality.h"
  34
  35 // Analysis task to fill histograms with PHOS ESD clusters and cells
  36 // Authors: Yuri Kharlov
  37 // Date   : 28.05.2009
  38
  39 ClassImp(AliAnalysisTaskPi0DiffEfficiency)
  40
  41 //________________________________________________________________________
  42 AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name)
  43 : AliAnalysisTaskSE(name),
  44   fStack(0),
  45   fOutputContainer(0),
  46   fPHOSEvent1(0),
  47   fPHOSEvent2(0),
  48   fPHOSCalibData(0),
  49   fNonLinCorr(0),
  50   fRPfull(0),
  51   fRPA(0),
  52   fRPC(0),
  53   fRPFar(0),
  54   fRPAFar(0),
  55   fRPCFar(0),
  56   fCentrality(0),
  57   fCenBin(0),
  58   fPHOSGeo(0),
  59   fEventCounter(0)
  60 {
  61   // Constructor
  62   for(Int_t i=0;i<1;i++){
  63     for(Int_t j=0;j<10;j++)
  64       for(Int_t k=0;k<11;k++)
  65         fPHOSEvents[i][j][k]=0 ;
  66   }
  67
  68   // Output slots #0 write into a TH1 container
  69   DefineOutput(1,TList::Class());
  70
  71   // Set bad channel map
  72   char key[55] ;
  73   for(Int_t i=0; i<6; i++){
  74     sprintf(key,"PHOS_BadMap_mod%d",i) ;
  75     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  76   }
  77   // Initialize the PHOS geometry
  78   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  79
  80
  81 }
  82
  83 //________________________________________________________________________
  84 void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
  85 {
  86   // Create histograms
  87   // Called once
  88
  89   // ESD histograms
  90   if(fOutputContainer != NULL){
  91     delete fOutputContainer;
  92   }
  93   fOutputContainer = new TList();
  94   fOutputContainer->SetOwner(kTRUE);
  95
  96   //Event selection
  97   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
  98
  99   //vertex distribution
 100   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
 101
 102   //Centrality
 103   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
 104
 105   //QA histograms
 106   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 107   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 108   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 109
 110   Int_t nM       = 500;
 111   Double_t mMin  = 0.0;
 112   Double_t mMax  = 1.0;
 113   Int_t nPt      = 200;
 114   Double_t ptMin = 0;
 115   Double_t ptMax = 20;
 116
 117   char key[55] ;
 118   for(Int_t cent=0; cent<6; cent++){
 119     //Single photon
 120     snprintf(key,55,"hPhotAll_cen%d",cent) ;
 121     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 122     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 123     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
 124     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 125     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 126     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
 127     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 128     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 129     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
 130     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 131     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 132
 133     snprintf(key,55,"hOldMassPtAll_cen%d",cent) ;
 134     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 135     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 136     snprintf(key,55,"hOldMassPtCPV_cen%d",cent) ;
 137     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 138     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 139     snprintf(key,55,"hOldMassPtDisp_cen%d",cent) ;
 140     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 141     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 142     snprintf(key,55,"hOldMassPtBoth_cen%d",cent) ;
 143     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 144     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 145
 146     snprintf(key,55,"hNewMassPtAll_cen%d",cent) ;
 147     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 148     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 149     snprintf(key,55,"hNewMassPtCPV_cen%d",cent) ;
 150     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 151     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 152     snprintf(key,55,"hNewMassPtDisp_cen%d",cent) ;
 153     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 154     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 155     snprintf(key,55,"hNewMassPtBoth_cen%d",cent) ;
 156     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 157     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 158
 159     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
 160     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 161     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 162     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
 163     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 164     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 165     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
 166     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 167     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 168     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
 169     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 170     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 171
 172     //Mixed
 173     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
 174     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 175     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
 176     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 177     snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
 178     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 179     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
 180     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 181
 182     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
 183     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 184     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
 185     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 186     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
 187     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 188     snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
 189     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 190
 191     //Single photon
 192     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
 193     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 194     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 195     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
 196     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 197     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 198     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
 199     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 200     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 201     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
 202     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 203     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 204
 205
 206
 207   }
 208
 209
 210   //MC
 211   for(Int_t cent=0; cent<6; cent++){
 212     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 213     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 214     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 215     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 216     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 217     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
 218     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 219     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 220     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 221     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 222     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 223     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 224     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 225     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 226     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 227     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 228     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 229     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 230     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 231     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 232     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 233     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 234     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 235     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 236   }
 237
 238   PostData(1, fOutputContainer);
 239
 240 }
 241
 242 //________________________________________________________________________
 243 void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *)
 244 {
 245   // Main loop, called for each event
 246   // Analyze ESD/AOD
 247
 248   FillHistogram("hSelEvents",0.5) ;
 249
 250   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 251   if (!event) {
 252     Printf("ERROR: Could not retrieve event");
 253     PostData(1, fOutputContainer);
 254     return;
 255   }
 256
 257   FillHistogram("hSelEvents",1.5) ;
 258   AliAODHeader *header = event->GetHeader() ;
 259
 260   // Checks if we have a primary vertex
 261   // Get primary vertices form ESD
 262   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 263
 264  // don't rely on ESD vertex, assume (0,0,0)
 265   Double_t vtx0[3] ={0.,0.,0.};
 266   Double_t vtx5[3] ={0.,0.,0.};
 267
 268   vtx5[0] = esdVertex5->GetX();
 269   vtx5[1] = esdVertex5->GetY();
 270   vtx5[2] = esdVertex5->GetZ();
 271
 272
 273   FillHistogram("hZvertex",esdVertex5->GetZ());
 274   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 275     PostData(1, fOutputContainer);
 276     return;
 277   }
 278   FillHistogram("hSelEvents",2.5) ;
 279
 280   //Vtx class z-bin
 281   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 282   //  if(zvtx<0)zvtx=0 ;
 283   //  if(zvtx>9)zvtx=9 ;
 284   Int_t zvtx=0 ;
 285
 286 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 287 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 288    fCentrality=header->GetZDCN2Energy() ;
 289
 290   if( fCentrality < 0. || fCentrality>80.){
 291     PostData(1, fOutputContainer);
 292     return;
 293   }
 294   FillHistogram("hSelEvents",3.5) ;
 295   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 296   fCenBin=0 ;
 297   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 298     fCenBin++ ;
 299
 300
 301   //reaction plain
 302   fRPfull= header->GetZDCN1Energy() ;
 303   if(fRPfull==999){ //reaction plain was not defined
 304     PostData(1, fOutputContainer);
 305     return;
 306   }
 307
 308   FillHistogram("hSelEvents",4.5) ;
 309   //All event selections done
 310   FillHistogram("hCentrality",fCentrality) ;
 311   //Reaction plain is defined in the range (-pi/2;pi/2)
 312   //We have 10 bins
 313   Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
 314   if(irp>9)irp=9 ;
 315
 316   if(!fPHOSEvents[zvtx][fCenBin][irp])
 317     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 318   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 319
 320   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 321   if(fEventCounter == 0) {
 322     for(Int_t mod=0; mod<5; mod++) {
 323       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 324       fPHOSGeo->SetMisalMatrix(m,mod) ;
 325       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 326     }
 327     fEventCounter++ ;
 328   }
 329
 330   ProcessMC() ;
 331
 332   if(fPHOSEvent1){
 333     fPHOSEvent1->Clear() ;
 334     fPHOSEvent2->Clear() ;
 335   }
 336   else{
 337     fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
 338     fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
 339   }
 340
 341   TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 342   TClonesArray * clustersOld = event->GetCaloClusters() ;
 343   TVector3 vertex(vtx0);
 344   char key[55] ;
 345   //Before Embedding
 346   Int_t multClustOld = clustersOld->GetEntriesFast();
 347   Int_t multClustEmb = clustersEmb->GetEntriesFast();
 348   Int_t inPHOSold=0 ;
 349   Int_t inPHOSemb=0 ;
 350   for (Int_t i=0; i<multClustOld; i++) {
 351     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
 352     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 353
 354     Bool_t survive=kFALSE ;
 355     for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
 356        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
 357        survive=IsSameCluster(clu,clu2);
 358     }
 359
 360
 361     Float_t  position[3];
 362     clu->GetPosition(position);
 363     TVector3 global(position) ;
 364     Int_t relId[4] ;
 365     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 366     Int_t mod  = relId[0] ;
 367     Int_t cellX = relId[2];
 368     Int_t cellZ = relId[3] ;
 369     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 370       continue ;
 371     if(clu->GetNCells()<3)
 372       continue ;
 373
 374     sprintf(key,"hCluM%d",mod) ;
 375     FillHistogram(key,cellX,cellZ,1.);
 376
 377     TLorentzVector pv1 ;
 378     clu->GetMomentum(pv1 ,vtx0);
 379
 380     if(inPHOSold>=fPHOSEvent1->GetSize()){
 381       fPHOSEvent1->Expand(inPHOSold+50) ;
 382     }
 383     AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 384     ph->SetModule(mod) ;
 385     ph->SetMomV2(&pv1) ;
 386     ph->SetNCells(clu->GetNCells());
 387     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
 388     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
 389     if(!survive) //this cluster found in list after embedding, skipping it
 390      ph->SetTagged(1) ;
 391
 392     inPHOSold++ ;
 393   }
 394
 395   for (Int_t i=0; i<multClustEmb; i++) {
 396     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
 397     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 398
 399     Bool_t survive=kFALSE ;
 400     for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
 401        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
 402        survive=IsSameCluster(clu,clu2);
 403     }
 404
 405     Float_t  position[3];
 406     clu->GetPosition(position);
 407     TVector3 global(position) ;
 408     Int_t relId[4] ;
 409     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 410     Int_t mod  = relId[0] ;
 411     Int_t cellX = relId[2];
 412     Int_t cellZ = relId[3] ;
 413     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 414       continue ;
 415     if(clu->GetNCells()<3)
 416       continue ;
 417
 418     sprintf(key,"hCluM%d",mod) ;
 419     FillHistogram(key,cellX,cellZ,1.);
 420
 421     TLorentzVector pv1 ;
 422     clu->GetMomentum(pv1 ,vtx0);
 423
 424     if(inPHOSemb>=fPHOSEvent2->GetSize()){
 425       fPHOSEvent2->Expand(inPHOSemb+50) ;
 426     }
 427     AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 428     ph->SetModule(mod) ;
 429     ph->SetMomV2(&pv1) ;
 430     ph->SetNCells(clu->GetNCells());
 431     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
 432     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
 433     if(!survive) //this cluster found in list after embedding, skipping it
 434      ph->SetTagged(1) ;
 435
 436     inPHOSemb++ ;
 437   }
 438
 439   //Single photon
 440   for (Int_t i1=0; i1<inPHOSold; i1++) {
 441     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 442     if(!ph1->IsTagged())
 443       continue ;
 444     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 445     FillHistogram(key,ph1->Pt(),-1.) ;
 446     if(ph1->IsCPVOK() ){
 447       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 448       FillHistogram(key,ph1->Pt(),-1.) ;
 449     }
 450     if(ph1->IsDispOK()){
 451       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 452       FillHistogram(key,ph1->Pt(),-1.) ;
 453       if(ph1->IsCPVOK()){
 454         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 455         FillHistogram(key,ph1->Pt(),-1.) ;
 456       }
 457     } // end of loop i2
 458   } // end of loop i1
 459
 460   for (Int_t i1=0; i1<inPHOSemb; i1++) {
 461     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
 462     if(!ph1->IsTagged())
 463       continue ;
 464     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 465     FillHistogram(key,ph1->Pt(),1.) ;
 466     if(ph1->IsCPVOK() ){
 467       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 468       FillHistogram(key,ph1->Pt(),1.) ;
 469     }
 470     if(ph1->IsDispOK()){
 471       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 472       FillHistogram(key,ph1->Pt(),1.) ;
 473       if(ph1->IsCPVOK()){
 474         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 475         FillHistogram(key,ph1->Pt(),1.) ;
 476       }
 477     } // end of loop i2
 478   } // end of loop i1
 479
 480
 481
 482   // Fill Real disribution:
 483   // Disappeared clusters enter with negative contribution
 484   // In addition fill control histogam with Real before embedding
 485   for (Int_t i1=0; i1<inPHOSold-1; i1++) {
 486     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 487     for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
 488       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
 489
 490       TLorentzVector p12  = *ph1  + *ph2;
 491       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 492       //Fill Controll histogram: Real before embedding
 493       snprintf(key,55,"hOldMassPtAll_cen%d",fCenBin) ;
 494       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 495       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 496         snprintf(key,55,"hOldMassPtCPV_cen%d",fCenBin) ;
 497         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 498       }
 499       if(ph1->IsDispOK() && ph2->IsDispOK()){
 500         snprintf(key,55,"hOldMassPtDisp_cen%d",fCenBin) ;
 501         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 502
 503         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 504           snprintf(key,55,"hOldMassPtBoth_cen%d",fCenBin) ;
 505           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 506         }
 507       }
 508
 509       //Now fill mail histograms with negative contributions
 510       if(!(ph1->IsTagged() || ph2->IsTagged()) )
 511         continue ;
 512       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
 513       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 514       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 515         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
 516         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 517       }
 518       if(ph1->IsDispOK() && ph2->IsDispOK()){
 519         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
 520         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 521
 522         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 523           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
 524           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 525         }
 526       }
 527     } // end of loop i2
 528   } // end of loop i1
 529
 530
 531   // Further fill Real disribution
 532   // now with positive contribution from new clusters
 533   // ass well fill controll histogram
 534   for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
 535     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
 536     for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
 537       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
 538
 539       TLorentzVector p12  = *ph1  + *ph2;
 540       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 541
 542       // Controll histogram: Real after embedding
 543       snprintf(key,55,"hNewMassPtAll_cen%d",fCenBin) ;
 544       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 545       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 546         snprintf(key,55,"hNewMassPtCPV_cen%d",fCenBin) ;
 547         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 548       }
 549       if(ph1->IsDispOK() && ph2->IsDispOK()){
 550         snprintf(key,55,"hNewMassPtDisp_cen%d",fCenBin) ;
 551         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 552
 553         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 554           snprintf(key,55,"hNewMassPtBoth_cen%d",fCenBin) ;
 555           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 556         }
 557       }
 558
 559       //Now fill main histogamm
 560       //new clusters with positive contribution
 561       if(!(ph1->IsTagged() || ph2->IsTagged()) )
 562         continue ;
 563       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
 564       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 565       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 566         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
 567         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 568       }
 569       if(ph1->IsDispOK() && ph2->IsDispOK()){
 570         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
 571         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 572
 573         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 574           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
 575           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
 576         }
 577       }
 578     } // end of loop i2
 579   } // end of loop i1
 580
 581
 582   //now mixed, does not really matter old or new list of clusters
 583   for (Int_t i1=0; i1<inPHOSemb; i1++) {
 584     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
 585     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 586       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 587       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 588         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 589         TLorentzVector p12  = *ph1  + *ph2;
 590         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 591
 592         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
 593         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 594         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 595           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
 596           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 597         }
 598         if(ph1->IsDispOK() && ph2->IsDispOK()){
 599           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
 600           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 601
 602           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 603             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
 604             FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
 605           }
 606         }
 607       } // end of loop i2
 608     }
 609   } // end of loop i1
 610
 611
 612   //Now we either add current events to stack or remove
 613   //If no photons in current event - no need to add it to mixed
 614   if(fPHOSEvent2->GetEntriesFast()>0){
 615     prevPHOS->AddFirst(fPHOSEvent2) ;
 616     fPHOSEvent2=0;
 617     delete fPHOSEvent1;
 618     fPHOSEvent1=0;
 619     if(prevPHOS->GetSize()>100){//Remove redundant events
 620       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 621       prevPHOS->RemoveLast() ;
 622       delete tmp ;
 623     }
 624   }
 625   // Post output data.
 626   PostData(1, fOutputContainer);
 627   fEventCounter++;
 628 }
 629
 630 //________________________________________________________________________
 631 void AliAnalysisTaskPi0DiffEfficiency::Terminate(Option_t *)
 632 {
 633   // Draw result to the screen
 634   // Called once at the end of the query
 635
 636 }
 637
 638 //________________________________________________________________________
 639 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 640 {
 641   //Check if this channel belogs to the good ones
 642
 643   if(strcmp(det,"PHOS")==0){
 644     if(mod>5 || mod<1){
 645       AliError(Form("No bad map for PHOS module %d ",mod)) ;
 646       return kTRUE ;
 647     }
 648     if(!fPHOSBadMap[mod]){
 649       AliError(Form("No Bad map for PHOS module %d",mod)) ;
 650       return kTRUE ;
 651     }
 652     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 653       return kFALSE ;
 654     else
 655       return kTRUE ;
 656   }
 657   else{
 658     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 659   }
 660   return kTRUE ;
 661 }
 662 //_____________________________________________________________________________
 663 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
 664   //FillHistogram
 665   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
 666   if(tmpI){
 667     tmpI->Fill(x) ;
 668     return ;
 669   }
 670   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
 671   if(tmpF){
 672     tmpF->Fill(x) ;
 673     return ;
 674   }
 675   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
 676   if(tmpD){
 677     tmpD->Fill(x) ;
 678     return ;
 679   }
 680   AliInfo(Form("can not find histogram <%s> ",key)) ;
 681 }
 682 //_____________________________________________________________________________
 683 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
 684   //FillHistogram
 685   TObject * tmp = fOutputContainer->FindObject(key) ;
 686   if(!tmp){
 687     AliInfo(Form("can not find histogram <%s> ",key)) ;
 688     return ;
 689   }
 690   if(tmp->IsA() == TClass::GetClass("TH1F")){
 691     ((TH1F*)tmp)->Fill(x,y) ;
 692     return ;
 693   }
 694   if(tmp->IsA() == TClass::GetClass("TH2F")){
 695     ((TH2F*)tmp)->Fill(x,y) ;
 696     return ;
 697   }
 698   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
 699 }
 700
 701 //_____________________________________________________________________________
 702 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
 703   //Fills 1D histograms with key
 704   TObject * tmp = fOutputContainer->FindObject(key) ;
 705   if(!tmp){
 706     AliInfo(Form("can not find histogram <%s> ",key)) ;
 707     return ;
 708   }
 709   if(tmp->IsA() == TClass::GetClass("TH2F")){
 710     ((TH2F*)tmp)->Fill(x,y,z) ;
 711     return ;
 712   }
 713   if(tmp->IsA() == TClass::GetClass("TH3F")){
 714     ((TH3F*)tmp)->Fill(x,y,z) ;
 715     return ;
 716   }
 717 }
 718 //_____________________________________________________________________________
 719 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda(Double_t l1,Double_t l2){
 720   Double_t l1Mean=1.22 ;
 721   Double_t l2Mean=2.0 ;
 722   Double_t l1Sigma=0.42 ;
 723   Double_t l2Sigma=0.71 ;
 724   Double_t c=-0.59 ;
 725   Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 726   return (R2<9.) ;
 727
 728 }
 729 //___________________________________________________________________________
 730 void AliAnalysisTaskPi0DiffEfficiency::ProcessMC(){
 731   //fill histograms for efficiensy etc. calculation
 732   const Double_t rcut = 1. ; //cut for primary particles
 733   //---------First pi0/eta-----------------------------
 734   char partName[10] ;
 735   char hkey[55] ;
 736
 737   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 738   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 739   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
 740      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
 741     if(particle->GetPdgCode() == 111)
 742       sprintf(partName,"pi0") ;
 743     else
 744       if(particle->GetPdgCode() == 221)
 745         sprintf(partName,"eta") ;
 746       else
 747         if(particle->GetPdgCode() == 22)
 748            sprintf(partName,"gamma") ;
 749         else
 750            continue ;
 751
 752     //Primary particle
 753     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
 754     if(r >rcut)
 755       continue ;
 756
 757     Double_t pt = particle->Pt() ;
 758     //Total number of pi0 with creation radius <1 cm
 759     sprintf(hkey,"hMC_all_%s_cen%d",partName,fCenBin) ;
 760     FillHistogram(hkey,pt) ;
 761     if(TMath::Abs(particle->Y())<0.12){
 762       sprintf(hkey,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
 763       FillHistogram(hkey,pt) ;
 764     }
 765
 766     sprintf(hkey,"hMC_rap_%s_cen%d",partName,fCenBin) ;
 767     FillHistogram(hkey,particle->Y()) ;
 768
 769     Double_t phi=particle->Phi() ;
 770     while(phi<0.)phi+=TMath::TwoPi() ;
 771     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 772     sprintf(hkey,"hMC_phi_%s_cen%d",partName,fCenBin) ;
 773     FillHistogram(hkey,phi) ;
 774
 775
 776     //Check if one of photons converted
 777     if(particle->GetNDaughters()!=2)
 778      continue ; //Do not account Dalitz decays
 779
 780 /*
 781     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
 782     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
 783     //Number of pi0s decayed into acceptance
 784     Int_t mod1,mod2 ;
 785     Double_t x=0.,z=0. ;
 786     Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
 787     Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
 788
 789     Bool_t goodPair=kFALSE ;
 790     if( hitPHOS1 && hitPHOS2){
 791       sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
 792       FillHistogram(hkey,pt) ;
 793       goodPair=kTRUE ;
 794     }
 795
 796 */
 797   }
 798
 799   //Now calculate "Real" distribution of clusters with primary
 800   TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
 801   Int_t multClust = event->GetNumberOfCaloClusters();
 802   Int_t inPHOS=0 ;
 803   Double_t vtx0[3] = {0,0,0};
 804   for (Int_t i=0; i<multClust; i++) {
 805     AliVCluster *clu = event->GetCaloCluster(i);
 806     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 807     if(clu->GetLabel()<0) continue ;
 808
 809     Float_t  position[3];
 810     clu->GetPosition(position);
 811     TVector3 global(position) ;
 812     Int_t relId[4] ;
 813     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 814     Int_t mod  = relId[0] ;
 815     Int_t cellX = relId[2];
 816     Int_t cellZ = relId[3] ;
 817     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 818       continue ;
 819     if(clu->GetNCells()<3)
 820       continue ;
 821
 822     TLorentzVector pv1 ;
 823     clu->GetMomentum(pv1 ,vtx0);
 824
 825     if(inPHOS>=cluPrim.GetSize()){
 826       cluPrim.Expand(inPHOS+50) ;
 827     }
 828     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 829     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 830     ph->SetModule(mod) ;
 831     ph->SetMomV2(&pv1) ;
 832     ph->SetNCells(clu->GetNCells());
 833     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
 834     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
 835
 836     inPHOS++ ;
 837
 838   }
 839
 840   //Single photon
 841   char key[55] ;
 842   for (Int_t i1=0; i1<inPHOS; i1++) {
 843     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
 844     snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
 845     FillHistogram(key,ph1->Pt()) ;
 846     if(ph1->IsCPVOK() ){
 847       snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
 848       FillHistogram(key,ph1->Pt()) ;
 849
 850     }
 851     if(ph1->IsDispOK()){
 852       snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
 853       FillHistogram(key,ph1->Pt()) ;
 854       if(ph1->IsCPVOK()){
 855         snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
 856         FillHistogram(key,ph1->Pt()) ;
 857       }
 858     } // end of loop i2
 859   } // end of loop i1
 860
 861   // Fill Real disribution
 862   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 863     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
 864     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 865       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
 866       TLorentzVector p12  = *ph1  + *ph2;
 867       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 868
 869       snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
 870       FillHistogram(key,p12.M() ,p12.Pt()) ;
 871       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 872         snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
 873         FillHistogram(key,p12.M() ,p12.Pt()) ;
 874       }
 875       if(ph1->IsDispOK() && ph2->IsDispOK()){
 876         snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
 877         FillHistogram(key,p12.M() ,p12.Pt()) ;
 878
 879         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 880           snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
 881           FillHistogram(key,p12.M() ,p12.Pt()) ;
 882         }
 883       }
 884     } // end of loop i2
 885   } // end of loop i1
 886
 887
 888 }
 889 //___________________________________________________________________________
 890 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
 891  //Compare clusters before and after embedding
 892  //clusters are the same if
 893  // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
 894  // - lists of digits are the same
 895
 896  if(c1->GetNCells() != c2->GetNCells())
 897    return kFALSE ;
 898
 899  if(TMath::Abs(c1->E()-c2->E())>0.001*c1->E())
 900    return kFALSE ;
 901
 902  UShort_t *list1 = c1->GetCellsAbsId() ;
 903  UShort_t *list2 = c2->GetCellsAbsId() ;
 904  for(Int_t i=0; i<c1->GetNCells(); i++){
 905   if(list1[i] != list2[i])
 906     return kFALSE ;
 907  }
 908  return kTRUE ;
 909
 910 }
 911
 912
 913