10 #include "TParticle.h"
15 #include "AliAODMCParticle.h"
16 #include "AliAnalysisManager.h"
17 #include "AliMCEventHandler.h"
18 #include "AliMCEvent.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"
32 #include "AliCDBManager.h"
33 #include "AliCentrality.h"
35 // Analysis task to fill histograms with PHOS ESD clusters and cells
36 // Authors: Yuri Kharlov
39 ClassImp(AliAnalysisTaskPi0DiffEfficiency)
41 //________________________________________________________________________
42 AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name)
43 : AliAnalysisTaskSE(name),
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 ;
68 // Output slots #0 write into a TH1 container
69 DefineOutput(1,TList::Class());
71 // Set bad channel map
73 for(Int_t i=0; i<6; i++){
74 snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
75 fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
77 // Initialize the PHOS geometry
78 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
83 //________________________________________________________________________
84 void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
90 if(fOutputContainer != NULL){
91 delete fOutputContainer;
93 fOutputContainer = new TList();
94 fOutputContainer->SetOwner(kTRUE);
97 fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
100 fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
103 fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
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));
118 for(Int_t cent=0; cent<6; cent++){
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() ;
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() ;
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() ;
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() ;
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));
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));
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() ;
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.)) ;
238 PostData(1, fOutputContainer);
242 //________________________________________________________________________
243 void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *)
245 // Main loop, called for each event
248 FillHistogram("hSelEvents",0.5) ;
250 AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
252 Printf("ERROR: Could not retrieve event");
253 PostData(1, fOutputContainer);
257 FillHistogram("hSelEvents",1.5) ;
258 AliAODHeader *header = event->GetHeader() ;
260 // Checks if we have a primary vertex
261 // Get primary vertices form ESD
262 const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
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.};
268 vtx5[0] = esdVertex5->GetX();
269 vtx5[1] = esdVertex5->GetY();
270 vtx5[2] = esdVertex5->GetZ();
273 FillHistogram("hZvertex",esdVertex5->GetZ());
274 if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
275 PostData(1, fOutputContainer);
278 FillHistogram("hSelEvents",2.5) ;
281 // Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
282 // if(zvtx<0)zvtx=0 ;
283 // if(zvtx>9)zvtx=9 ;
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() ;
290 if( fCentrality < 0. || fCentrality>80.){
291 PostData(1, fOutputContainer);
294 FillHistogram("hSelEvents",3.5) ;
295 Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
297 while(fCenBin<6 && fCentrality > bins[fCenBin+1])
302 fRPfull= header->GetZDCN1Energy() ;
303 if(fRPfull==999){ //reaction plain was not defined
304 PostData(1, fOutputContainer);
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)
313 Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
316 if(!fPHOSEvents[zvtx][fCenBin][irp])
317 fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
318 TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
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);
333 fPHOSEvent1->Clear() ;
334 fPHOSEvent2->Clear() ;
337 fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
338 fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
341 TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
342 TClonesArray * clustersOld = event->GetCaloClusters() ;
343 TVector3 vertex(vtx0);
346 Int_t multClustOld = clustersOld->GetEntriesFast();
347 Int_t multClustEmb = clustersEmb->GetEntriesFast();
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;
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);
362 clu->GetPosition(position);
363 TVector3 global(position) ;
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) )
371 if(clu->GetNCells()<3)
374 snprintf(key,55,"hCluM%d",mod) ;
375 FillHistogram(key,cellX,cellZ,1.);
378 clu->GetMomentum(pv1 ,vtx0);
380 if(inPHOSold>=fPHOSEvent1->GetSize()){
381 fPHOSEvent1->Expand(inPHOSold+50) ;
383 AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
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
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;
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);
406 clu->GetPosition(position);
407 TVector3 global(position) ;
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) )
415 if(clu->GetNCells()<3)
418 snprintf(key,55,"hCluM%d",mod) ;
419 FillHistogram(key,cellX,cellZ,1.);
422 clu->GetMomentum(pv1 ,vtx0);
424 if(inPHOSemb>=fPHOSEvent2->GetSize()){
425 fPHOSEvent2->Expand(inPHOSemb+50) ;
427 AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
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
440 for (Int_t i1=0; i1<inPHOSold; i1++) {
441 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
444 snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
445 FillHistogram(key,ph1->Pt(),-1.) ;
447 snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
448 FillHistogram(key,ph1->Pt(),-1.) ;
451 snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
452 FillHistogram(key,ph1->Pt(),-1.) ;
454 snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
455 FillHistogram(key,ph1->Pt(),-1.) ;
460 for (Int_t i1=0; i1<inPHOSemb; i1++) {
461 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
464 snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
465 FillHistogram(key,ph1->Pt(),1.) ;
467 snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
468 FillHistogram(key,ph1->Pt(),1.) ;
471 snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
472 FillHistogram(key,ph1->Pt(),1.) ;
474 snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
475 FillHistogram(key,ph1->Pt(),1.) ;
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) ;
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) ;
499 if(ph1->IsDispOK() && ph2->IsDispOK()){
500 snprintf(key,55,"hOldMassPtDisp_cen%d",fCenBin) ;
501 FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
503 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
504 snprintf(key,55,"hOldMassPtBoth_cen%d",fCenBin) ;
505 FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
509 //Now fill mail histograms with negative contributions
510 if(!(ph1->IsTagged() || ph2->IsTagged()) )
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) ;
518 if(ph1->IsDispOK() && ph2->IsDispOK()){
519 snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
520 FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
522 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
523 snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
524 FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
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) ;
539 TLorentzVector p12 = *ph1 + *ph2;
540 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
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) ;
549 if(ph1->IsDispOK() && ph2->IsDispOK()){
550 snprintf(key,55,"hNewMassPtDisp_cen%d",fCenBin) ;
551 FillHistogram(key,p12.M() ,p12.Pt(),1) ;
553 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
554 snprintf(key,55,"hNewMassPtBoth_cen%d",fCenBin) ;
555 FillHistogram(key,p12.M() ,p12.Pt(),1) ;
559 //Now fill main histogamm
560 //new clusters with positive contribution
561 if(!(ph1->IsTagged() || ph2->IsTagged()) )
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) ;
569 if(ph1->IsDispOK() && ph2->IsDispOK()){
570 snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
571 FillHistogram(key,p12.M() ,p12.Pt(),1) ;
573 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
574 snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
575 FillHistogram(key,p12.M() ,p12.Pt(),1) ;
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());
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.) ;
598 if(ph1->IsDispOK() && ph2->IsDispOK()){
599 snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
600 FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
602 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
603 snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
604 FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
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) ;
619 if(prevPHOS->GetSize()>100){//Remove redundant events
620 TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
621 prevPHOS->RemoveLast() ;
626 PostData(1, fOutputContainer);
630 //________________________________________________________________________
631 void AliAnalysisTaskPi0DiffEfficiency::Terminate(Option_t *)
633 // Draw result to the screen
634 // Called once at the end of the query
638 //________________________________________________________________________
639 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
641 //Check if this channel belogs to the good ones
643 if(strcmp(det,"PHOS")==0){
645 AliError(Form("No bad map for PHOS module %d ",mod)) ;
648 if(!fPHOSBadMap[mod]){
649 AliError(Form("No Bad map for PHOS module %d",mod)) ;
652 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
658 AliError(Form("Can not find bad channels for detector %s ",det)) ;
662 //_____________________________________________________________________________
663 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
665 TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
670 TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
675 TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
680 AliInfo(Form("can not find histogram <%s> ",key)) ;
682 //_____________________________________________________________________________
683 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
685 TObject * tmp = fOutputContainer->FindObject(key) ;
687 AliInfo(Form("can not find histogram <%s> ",key)) ;
690 if(tmp->IsA() == TClass::GetClass("TH1F")){
691 ((TH1F*)tmp)->Fill(x,y) ;
694 if(tmp->IsA() == TClass::GetClass("TH2F")){
695 ((TH2F*)tmp)->Fill(x,y) ;
698 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
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) ;
706 AliInfo(Form("can not find histogram <%s> ",key)) ;
709 if(tmp->IsA() == TClass::GetClass("TH2F")){
710 ((TH2F*)tmp)->Fill(x,y,z) ;
713 if(tmp->IsA() == TClass::GetClass("TH3F")){
714 ((TH3F*)tmp)->Fill(x,y,z) ;
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 ;
725 Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
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-----------------------------
737 AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
739 TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
740 for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
741 AliAODMCParticle* particle = (AliAODMCParticle*) mcArray->At(i);
742 if(particle->GetPdgCode() == 111)
743 snprintf(partName,10,"pi0") ;
745 if(particle->GetPdgCode() == 221)
746 snprintf(partName,10,"eta") ;
748 if(particle->GetPdgCode() == 22)
749 snprintf(partName,10,"gamma") ;
754 Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
758 Double_t pt = particle->Pt() ;
759 //Total number of pi0 with creation radius <1 cm
760 snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
761 FillHistogram(hkey,pt) ;
762 if(TMath::Abs(particle->Y())<0.12){
763 snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
764 FillHistogram(hkey,pt) ;
767 snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
768 FillHistogram(hkey,particle->Y()) ;
770 Double_t phi=particle->Phi() ;
771 while(phi<0.)phi+=TMath::TwoPi() ;
772 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
773 snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
774 FillHistogram(hkey,phi) ;
777 //Check if one of photons converted
778 if(particle->GetNDaughters()!=2)
779 continue ; //Do not account Dalitz decays
782 TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
783 TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
784 //Number of pi0s decayed into acceptance
787 Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
788 Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
790 Bool_t goodPair=kFALSE ;
791 if( hitPHOS1 && hitPHOS2){
792 sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
793 FillHistogram(hkey,pt) ;
800 //Now calculate "Real" distribution of clusters with primary
801 TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
802 Int_t multClust = event->GetNumberOfCaloClusters();
804 Double_t vtx0[3] = {0,0,0};
805 for (Int_t i=0; i<multClust; i++) {
806 AliVCluster *clu = event->GetCaloCluster(i);
807 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
808 if(clu->GetLabel()<0) continue ;
811 clu->GetPosition(position);
812 TVector3 global(position) ;
814 fPHOSGeo->GlobalPos2RelId(global,relId) ;
815 Int_t mod = relId[0] ;
816 Int_t cellX = relId[2];
817 Int_t cellZ = relId[3] ;
818 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
820 if(clu->GetNCells()<3)
824 clu->GetMomentum(pv1 ,vtx0);
826 if(inPHOS>=cluPrim.GetSize()){
827 cluPrim.Expand(inPHOS+50) ;
829 AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
830 //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
833 ph->SetNCells(clu->GetNCells());
834 ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
835 ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
843 for (Int_t i1=0; i1<inPHOS; i1++) {
844 AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
845 snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
846 FillHistogram(key,ph1->Pt()) ;
848 snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
849 FillHistogram(key,ph1->Pt()) ;
853 snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
854 FillHistogram(key,ph1->Pt()) ;
856 snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
857 FillHistogram(key,ph1->Pt()) ;
862 // Fill Real disribution
863 for (Int_t i1=0; i1<inPHOS-1; i1++) {
864 AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
865 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
866 AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
867 TLorentzVector p12 = *ph1 + *ph2;
868 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
870 snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
871 FillHistogram(key,p12.M() ,p12.Pt()) ;
872 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
873 snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
874 FillHistogram(key,p12.M() ,p12.Pt()) ;
876 if(ph1->IsDispOK() && ph2->IsDispOK()){
877 snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
878 FillHistogram(key,p12.M() ,p12.Pt()) ;
880 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
881 snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
882 FillHistogram(key,p12.M() ,p12.Pt()) ;
890 //___________________________________________________________________________
891 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
892 //Compare clusters before and after embedding
893 //clusters are the same if
894 // - Energy changed less than 0.1% (numerical accuracy in reconstruction)
895 // - lists of digits are the same
897 if(c1->GetNCells() != c2->GetNCells())
900 if(TMath::Abs(c1->E()-c2->E())>0.001*c1->E())
903 UShort_t *list1 = c1->GetCellsAbsId() ;
904 UShort_t *list2 = c2->GetCellsAbsId() ;
905 for(Int_t i=0; i<c1->GetNCells(); i++){
906 if(list1[i] != list2[i])