| 02f0b8e7 |
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 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame