1 #include "AliAnalysisTaskSDMGammaMC.h"
17 #include <TLorentzVector.h>
21 #include "AliAnalysisTaskSE.h"
22 #include "AliAnalysisManager.h"
24 #include "AliESDtrackCuts.h"
25 #include "AliESDEvent.h"
26 #include "AliESDInputHandler.h"
27 #include "AliAODEvent.h"
28 #include "AliMCEvent.h"
29 #include "AliEMCALGeometry.h"
30 #include "AliInputEventHandler.h"
31 #include "AliESDInputHandler.h"
32 #include "AliAODInputHandler.h"
34 #include "AliEMCALRecoUtils.h"
35 #include "AliExternalTrackParam.h"
38 #include <TGeoManager.h>
39 #include <TGeoMatrix.h>
44 #include <TObjArray.h>
47 #include "AliVCluster.h"
48 #include "AliVCaloCells.h"
51 #include "AliESDEvent.h"
52 #include "AliAODEvent.h"
53 #include "AliESDtrack.h"
54 #include "AliAODTrack.h"
55 #include "AliExternalTrackParam.h"
56 #include "AliESDfriendTrack.h"
57 #include "AliTrackerBase.h"
60 #include "AliEMCALRecoUtils.h"
61 #include "AliEMCALGeometry.h"
62 #include "AliTrackerBase.h"
63 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
64 #include "AliEMCALPIDUtils.h"
69 ClassImp(AliAnalysisTaskSDMGammaMC)
71 //________________________________________________________________________
72 AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC() :
92 h3_gE_RecTruth_ncellscut(0),
95 h1_Pi0TruthPtEmcal(0),
96 h1_PriPi0TruthPtEmcal(0),
98 h2_PriPi0TruthPhiEta(0),
99 h2_Pi0TruthPhiEtaEmcal(0),
100 h2_PriPi0TruthPhiEtaEmcal(0),
101 h1_TruthPhotonsEmcal(0),
102 h2_TruthPhotonsPhiEta(0),
104 h1_PhotonsNCellsCut(0),
105 h1_PhotonsTrackMatchCut(0),
107 h2_PhotonsPhiEtaIsEmcal(0),
122 // Dummy constructor ALWAYS needed for I/O.
125 //________________________________________________________________________
126 AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC(const char *name) :
127 AliAnalysisTaskSE(name),
146 h3_gE_RecTruth_ncellscut(0),
149 h1_Pi0TruthPtEmcal(0),
150 h1_PriPi0TruthPtEmcal(0),
151 h2_Pi0TruthPhiEta(0),
152 h2_PriPi0TruthPhiEta(0),
153 h2_Pi0TruthPhiEtaEmcal(0),
154 h2_PriPi0TruthPhiEtaEmcal(0),
155 h1_TruthPhotonsEmcal(0),
156 h2_TruthPhotonsPhiEta(0),
158 h1_PhotonsNCellsCut(0),
159 h1_PhotonsTrackMatchCut(0),
161 h2_PhotonsPhiEtaIsEmcal(0),
177 // Define input and output slots here (never in the dummy constructor)
178 // Input slot #0 works with a TChain - it is connected to the default input container
179 // Output slot #1 writes into a TH1 container
182 DefineOutput(1, TList::Class()); // for output list
185 //________________________________________________________________________
186 AliAnalysisTaskSDMGammaMC::~AliAnalysisTaskSDMGammaMC()
188 // Destructor. Clean-up the output list, but not the histograms that are put inside
189 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
190 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
196 //________________________________________________________________________
197 void AliAnalysisTaskSDMGammaMC::UserCreateOutputObjects()
200 // Called once (on the worker node)
202 fOutput = new TList();
203 fOutput->SetOwner(); // IMPORTANT!
205 fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
207 cout << "__________AliAnalysisTaskSDMGammaMC: Input settings__________" << endl;
208 cout << " fMcMode: " << fMcMode << endl;
209 cout << " fRecalibrator: " << fRecalibrator << endl;
210 cout << " phi range: " << fPhimin << ", " << fPhimax << endl;
211 cout << " eta range: " << fEtamin << ", " << fEtamax << endl;
212 cout << " number of zvtx bins: " << zvtx_bins << endl;
213 cout << " number of mult bins: " << mult_bins << endl;
214 cout << " poolDepth: " << poolDepth << endl;
218 double TotalNBins = 0.0;
221 Int_t nClustersbins = 501;
222 Float_t nClusterslow = -0.5, nClustersup = 500.5;
223 h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup);
224 h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt");
225 h1_nClusters->GetYaxis()->SetTitle("counts");
226 h1_nClusters->SetMarkerStyle(kFullCircle);
227 TotalNBins+=nClustersbins;
229 Int_t nZvertexbins = 501;
230 Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
231 h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
232 h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
233 h1_zvtx->GetYaxis()->SetTitle("counts");
234 h1_zvtx->SetMarkerStyle(kFullCircle);
235 TotalNBins+=nZvertexbins;
237 h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
241 Float_t ptlow = 0.0, ptup = 20.0;
243 Float_t Elow = 0.0, Eup = 20.0;
244 h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup);
245 h1_E->GetXaxis()->SetTitle("E [GeV]");
246 h1_E->GetYaxis()->SetTitle("counts");
247 h1_E->SetMarkerStyle(kFullCircle);
250 h1_Phi = new TH1F("h1_Phi", "phi distribution", 1000, -7, 7);
251 h1_Phi->GetXaxis()->SetTitle("#phi [rad]");
252 h1_Phi->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
253 h1_Phi->SetMarkerStyle(kFullCircle);
256 h2_PiMotherID = new TH2F("h2_PiMotherID", "Mother ID for Truth Pi0's", 100001, -0.5,100000.5, 2, 0.5,2.5);
257 h2_PiMotherID->GetXaxis()->SetTitle("#pi^{0} Mother Particle ID");
258 h2_PiMotherID->GetYaxis()->SetTitle("primary or non-primary");
259 h2_PiMotherID->GetZaxis()->SetTitle("counts");
260 h2_PiMotherID->SetMarkerStyle(kFullCircle);
261 TotalNBins+=2*100001;
263 h2_GaMotherID = new TH2F("h2_GaMotherID", "Mother ID for Truth #gamma's", 100001, -0.5,100000.5, 2, 0.5,2.5);
264 h2_GaMotherID->GetXaxis()->SetTitle("#gamma Mother Particle ID");
265 h2_GaMotherID->GetYaxis()->SetTitle("primary or non-primary");
266 h2_GaMotherID->GetZaxis()->SetTitle("counts");
267 h2_GaMotherID->SetMarkerStyle(kFullCircle);
268 TotalNBins+=2*100001;
270 h3_gE_RecTruth = new TH3F("h3_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
271 h3_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
272 h3_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
273 h3_gE_RecTruth->GetZaxis()->SetTitle("category");
274 h3_gE_RecTruth->SetMarkerStyle(kFullCircle);
275 TotalNBins+=Ebins*1000*4;
277 h3_gE_RecTruth_ncellscut = new TH3F("h3_gE_RecTruth_ncellscut", "#gamma E_{truth}/E_{clust} vs E_{clust} (for nCells>1)", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
278 h3_gE_RecTruth_ncellscut->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
279 h3_gE_RecTruth_ncellscut->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
280 h3_gE_RecTruth_ncellscut->GetZaxis()->SetTitle("category");
281 h3_gE_RecTruth_ncellscut->SetMarkerStyle(kFullCircle);
282 TotalNBins+=Ebins*1000*4;
284 h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup);
285 h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
286 h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
287 h1_Pi0TruthPt->SetMarkerStyle(kFullCircle);
290 h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup);
291 h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
292 h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
293 h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle);
296 h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup);
297 h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
298 h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
299 h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
302 h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
303 h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
304 h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
305 h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
308 h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
309 h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
310 h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
313 h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
314 h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
315 h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
318 h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
319 h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
320 h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
323 h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
324 h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
325 h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
328 h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
329 h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
330 h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
331 h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle);
334 h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
335 "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
336 h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
337 h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
338 h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
341 h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
342 h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
343 h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
344 h1_PhotonsEmcal->SetMarkerStyle(kFullCircle);
347 h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup);
348 h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
349 h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
350 h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle);
353 h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup);
354 h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
355 h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
356 h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle);
359 h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup);
360 h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
361 h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
362 h1_PhotonsAllCut->SetMarkerStyle(kFullCircle);
365 h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
366 "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 150,-1.5,1.5);
367 h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
368 h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
371 h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10);
372 h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
373 h1_dR_RealMC->GetYaxis()->SetTitle("N");
374 h1_dR_RealMC->SetMarkerStyle(kFullCircle);
378 Float_t etalow = -1.5, etaup = 1.5;
379 h1_Eta = new TH1F("h1_Eta","#eta distribution for reconstructed",etabins, etalow, etaup);
380 h1_Eta->GetXaxis()->SetTitle("#eta");
381 h1_Eta->GetYaxis()->SetTitle("counts");
384 Int_t chi2bins = 100;
385 Float_t chi2low = -2, chi2up = 2;
386 h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up);
387 h1_Chi2->GetXaxis()->SetTitle("#chi^{2}");
388 h1_Chi2->GetYaxis()->SetTitle("counts");
389 TotalNBins+=chi2bins;
391 h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5);
392 h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched");
393 h1_nTrkMatch->GetYaxis()->SetTitle("counts");
396 h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3);
397 h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
398 h1_ClusterDisp->GetYaxis()->SetTitle("counts");
401 h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1);
402 h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()");
403 h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()");
404 h2_Ellipse->GetZaxis()->SetTitle("counts");
407 h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup);
408 h2_EtaPt->GetXaxis()->SetTitle("E [GeV]");
409 h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
410 TotalNBins+=etabins*ptbins;
412 h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0);
413 h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi");
414 h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta");
417 h2_dphi_deta_mix = new TH2F("h2_dphi_deta_mix","#Delta#phi vs #Delta#eta (mixed events)", 349,-1.5,5, 400,-2.0,2.0);
418 h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi");
419 h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta");
422 h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2);
423 h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
424 h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
425 h2_DispRes->GetZaxis()->SetTitle("counts");
428 h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5);
429 h2_cells_M02->GetXaxis()->SetTitle("nCells");
430 h2_cells_M02->GetYaxis()->SetTitle("M02");
431 h2_cells_M02->GetZaxis()->SetTitle("counts");
434 cout << endl << "Total number of bins in booked histograms: " << TotalNBins << endl << endl;
436 //TFile *f = OpenFile(1);
437 //TDirectory::TContext context(f);
439 fOutput->Add(h1_nClusters);
440 fOutput->Add(h1_zvtx);
441 fOutput->Add(h1_trigger);
443 fOutput->Add(h1_Phi);
444 fOutput->Add(h2_PiMotherID);
445 fOutput->Add(h2_GaMotherID);
446 fOutput->Add(h3_gE_RecTruth);
447 fOutput->Add(h3_gE_RecTruth_ncellscut);
448 fOutput->Add(h1_Pi0TruthPt);
449 fOutput->Add(h1_PriPi0TruthPt);
450 fOutput->Add(h1_Pi0TruthPtEmcal);
451 fOutput->Add(h1_PriPi0TruthPtEmcal);
452 fOutput->Add(h2_Pi0TruthPhiEta);
453 fOutput->Add(h2_PriPi0TruthPhiEta);
454 fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
455 fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
456 fOutput->Add(h1_TruthPhotonsEmcal);
457 fOutput->Add(h2_TruthPhotonsPhiEta);
458 fOutput->Add(h1_PhotonsEmcal);
459 fOutput->Add(h1_PhotonsNCellsCut);
460 fOutput->Add(h1_PhotonsTrackMatchCut);
461 fOutput->Add(h1_PhotonsAllCut);
462 fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
463 fOutput->Add(h1_dR_RealMC);
464 fOutput->Add(h1_Eta);
465 fOutput->Add(h1_Chi2);
466 fOutput->Add(h1_nTrkMatch);
467 fOutput->Add(h1_ClusterDisp);
468 fOutput->Add(h2_Ellipse);
469 fOutput->Add(h2_EtaPt);
470 fOutput->Add(h2_dphi_deta);
471 fOutput->Add(h2_dphi_deta_mix);
472 fOutput->Add(h2_DispRes);
473 fOutput->Add(h2_cells_M02);
475 // Post data for ALL output slots >0 here,
476 // To get at least an empty histogram
477 // 1 is the outputnumber of a certain weg of task 1
478 PostData(1, fOutput);
481 //________________________________________________________________________
482 void AliAnalysisTaskSDMGammaMC::UserExec(Option_t *)
484 // Main loop Called for each event
486 AliMCEvent *mcEvent = MCEvent();
487 Bool_t isMC = bool(mcEvent);//is this the right way to do this?
489 cout << "no MC event" << endl;
493 TRandom3 randy; randy.SetSeed(0);
495 double recalScale = 1.0;
497 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
499 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
500 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
501 if (!aodH && !esdH) Printf("ERROR: Could not get ESD or AODInputHandler");
503 if(esdH) fEsdEv = esdH->GetEvent();
504 else if(aodH) fAodEv = aodH->GetEvent();
506 AliFatal("Neither ESD nor AOD event found");
511 // get pointer to reconstructed event
512 AliVEvent *event = InputEvent();
514 AliError("Pointer == 0, this can not happen!"); return;}
515 //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
516 //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
518 //AliError("Cannot get the ESD event"); return;}
521 if (fEsdEv) iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
522 else if (fAodEv) iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
523 //h1_trigger->Fill(iTrigger);
526 Int_t iTriggerBin = 0;
527 for(unsigned long j=0; j<TriggerList.size(); j++){
528 if(iTrigger==TriggerList[j])
532 TriggerList.push_back(iTrigger);
533 iTriggerBin=TriggerList.size();
536 h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
537 sprintf(saythis,"%d",iTrigger);
538 h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
541 TString trigClasses = fEsdEv->GetFiredTriggerClasses();
542 // remove "fast cluster events":
543 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
547 TString trigClasses = fAodEv->GetFiredTriggerClasses();
548 // remove "fast cluster events":
549 if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL"))
554 if(!(fEsdEv->GetPrimaryVertex()->GetStatus())) return;
557 //if(!(fAodEv->GetPrimaryVertex()->GetStatus())) return;
560 Double_t vertDisp=0.0;
561 Double_t vertZres=0.0;
562 Bool_t vertIsfromZ=0;
564 vertDisp = fEsdEv->GetPrimaryVertex()->GetDispersion();
565 vertZres = fEsdEv->GetPrimaryVertex()->GetZRes();
566 vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
574 h2_DispRes->Fill(vertDisp, vertZres);
575 // if vertex is from spd vertexZ, require more stringent cut
577 if (vertDisp>0.02 || vertZres>0.25 )
578 return; // bad vertex from VertexerZ
583 //Int_t evtN = fEsdEv->GetEventNumberInFile();
584 //Int_t ntracks = fEsdEv->GetNumberOfTracks();
585 nclusters = fEsdEv->GetNumberOfCaloClusters();
588 //Int_t evtN = fAodEv->GetEventNumberInFile();
589 //Int_t ntracks = fAodEv->GetNumberOfTracks();
590 nclusters = fAodEv->GetNumberOfCaloClusters();
593 // EMCal cluster loop for reconstructed event
594 //numberofclusters set above!
596 if (fEsdEv) vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
597 else if (fAodEv) vertZ = fAodEv->GetPrimaryVertex()->GetZ();
599 h1_zvtx->Fill(vertZ);
604 h1_nClusters->Fill(nclusters);
606 //cout << iskip << " " << izvtx << " " << imult << endl;
607 //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
611 //######################### ~~~~~~~~~~~ ##################################
612 //######################### STARTING MC ##################################
613 //######################### ~~~~~~~~~~~ ##################################
620 cout << "no MC event" << endl;
624 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
628 mcEvent->PreReadAll();
630 Int_t nTracksMC = mcEvent->GetNumberOfTracks();
631 Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
633 for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
634 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
638 if(iTrack<nPTracksMC) isPrimary = 1;
641 if(mcP->PdgCode() == 22){
643 if(mcP->GetMother()>-1)
644 h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
646 h2_GaMotherID->Fill(0.0,1);
649 h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
653 if(mcP->PdgCode() != 111)
657 if(mcP->GetMother()>-1){
658 if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 310 ||
659 ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310 )
664 //Double_t dR_vtx = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
665 // (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
666 //if(dR_vtx <= 0.01) isPrimary = 1;
667 //else isPrimary = 0;
671 if(mcP->GetMother()>-1)
672 h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
674 h2_PiMotherID->Fill(0.0,1);
677 h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
679 h1_Pi0TruthPt ->Fill(mcP->Pt());
680 h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
683 h1_PriPi0TruthPt ->Fill(mcP->Pt());
684 h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
687 if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
690 Int_t DecayPhotonLabel[2] = {mcP->GetFirstDaughter(),
691 mcP->GetLastDaughter() };
693 if (DecayPhotonLabel[0]<0) continue;
694 if (DecayPhotonLabel[1]<0) DecayPhotonLabel[1]=DecayPhotonLabel[0];
695 if (DecayPhotonLabel[1]-DecayPhotonLabel[0] != 1) continue;
699 bool isConv[2] = {1,1};
700 Int_t convIndices[2][2] = { {-1,-1},{-1,-1} };
701 Double_t eta_d[2] = {0.0,0.0};
702 Double_t phi_d[2] = {0.0,0.0};
703 Int_t daughter_index = -1;
704 for (Int_t iPhoton=DecayPhotonLabel[0];iPhoton<=DecayPhotonLabel[1];++iPhoton){
705 if(iPhoton==DecayPhotonLabel[0]) daughter_index=0;
706 else daughter_index=1;
707 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(iPhoton));
708 eta_d[daughter_index] = dmc->Eta();
709 phi_d[daughter_index] = dmc->Phi();
710 if(!(dmc->PdgCode()==22)) binp = false;
711 if(!(eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax &&
712 phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax )) bacc = false;
714 if( ((TParticle*)dmc->Particle())->GetNDaughters() != 2 ) isConv[daughter_index] = 0;
715 else{//if photon has 2 daughters.
717 Int_t dd1 = dmc->GetFirstDaughter();
718 Int_t dd2 = dmc->GetLastDaughter();
719 if (dd2-dd1 != 1) cout << "How can this happen???? " << endl;
720 const AliMCParticle *dd1mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd1));
721 const AliMCParticle *dd2mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd2));
722 if( dd1mc->PdgCode() != -dd2mc->PdgCode() )
723 isConv[daughter_index] = 0;
724 else if( TMath::Abs(dd1mc->PdgCode())!=11 )
725 isConv[daughter_index] = 0;
726 if(isConv[daughter_index]==1){//store the e+e- indices...
727 convIndices[daughter_index][0] = dd1;
728 convIndices[daughter_index][1] = dd2;
729 }//if this photon converted.
730 }//close else-if photon has 2 daughters.
731 }//loop over 2 decay photons (iPhoton)
733 if(binp && bacc){// 2 Photons hit the EMCAL!
735 h1_Pi0TruthPtEmcal ->Fill(mcP->Pt());
736 h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
739 h1_PriPi0TruthPtEmcal ->Fill(mcP->Pt());
740 h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
743 Int_t PhotonClusterMatch[2][3] = { {0,-1,-1},
745 Int_t PhotonElectronMatch[2][6] = { {0,-1,-1,-1,-1,-1},
746 {0,-1,-1,-1,-1,-1} };
748 for(int iCluster=0; iCluster<nclusters; iCluster++) {
750 AliESDCaloCluster* esdCluster=NULL;
751 AliAODCaloCluster* aodCluster=NULL;
752 if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
753 else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
755 Double_t clustMC_phi, clustMC_eta;
757 if(esdCluster->IsEMCAL()){
759 if(!isGoodEsdCluster(esdCluster))
762 Float_t pos[3] = {0,0,0};
763 esdCluster->GetPosition(pos);
765 h1_Phi->Fill(vpos.Phi());
766 clustMC_phi = vpos.Phi();
767 clustMC_eta = vpos.Eta();
769 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
770 (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
771 h1_dR_RealMC->Fill(dR);
772 //matches_pion_photon = 0;
773 //if(dR<=0.04) matches_pion_photon = 1;
775 TArrayI *TruthLabelsA = esdCluster->GetLabelsArray();
777 for(int itl=0; itl<TruthLabelsA->GetSize(); itl++){
779 for(int iPhoton=0; iPhoton<2; iPhoton++){
780 if(TruthLabelsA->At(itl)==DecayPhotonLabel[iPhoton]){
781 PhotonClusterMatch[iPhoton][0] = 1;
782 PhotonClusterMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
783 PhotonClusterMatch[iPhoton][2] = iCluster;
785 }//loop over truth labels.
787 AliMCParticle *elecCandidate = (AliMCParticle*)(mcEvent->GetTrack(TruthLabelsA->At(itl)));
788 if(TMath::Abs(elecCandidate->PdgCode())==11){//if we have an electron...
789 Int_t elecMother_index = elecCandidate->GetMother();
790 if(elecMother_index>1 && elecMother_index<nTracksMC){
791 AliMCParticle *elecMother = (AliMCParticle*)(mcEvent->GetTrack(elecMother_index));
792 if( TMath::Abs(elecMother->PdgCode())==22 ){//if the e's mother is a photon...
793 Int_t elecGrandMother_index = elecMother->GetMother();
794 if(elecGrandMother_index==iTrack){//if the e's gMother is THE pi0 in question...
795 AliMCParticle *elecGrandMother = (AliMCParticle*)(mcEvent->GetTrack(elecGrandMother_index));
796 if( TMath::Abs(elecGrandMother->PdgCode())!=111 ) cout << "|| This can't happen!! A pion is a pion is a pion is a pion... ||" << endl;
798 for(int iPhoton=0; iPhoton<2; iPhoton++){
799 //if(convIndices[iPhoton][0]==elecMother_index){
800 if(convIndices[iPhoton][0]==TruthLabelsA->At(itl) || convIndices[iPhoton][1]==TruthLabelsA->At(itl)){
801 if(PhotonElectronMatch[iPhoton][1] == DecayPhotonLabel[iPhoton]) PhotonElectronMatch[iPhoton][0] = 2;
802 else PhotonElectronMatch[iPhoton][0] = 1;
803 PhotonElectronMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
804 if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][2] = iCluster;//first cluster
805 else PhotonElectronMatch[iPhoton][3] = iCluster;//second cluster
806 if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][4] = elecCandidate->PdgCode();
807 else PhotonElectronMatch[iPhoton][5] = elecCandidate->PdgCode();
809 }//loop over both decay photons
812 }//if we have a photon.
813 }//if we have an electron.
814 }//if the candidate has a real mother.
816 }//itl (TruthLabel loop)
817 }//if(TruthLabelsA exists)
820 }//loop over clusters.
823 for(int iPhoton=0; iPhoton<2; iPhoton++){
824 AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(DecayPhotonLabel[iPhoton]));
827 if(PhotonClusterMatch[iPhoton][0]==1){
828 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonClusterMatch[iPhoton][2]);
829 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
830 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
831 if(esdCluster->GetNCells()>=2)
832 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
834 else if(PhotonElectronMatch[iPhoton][0]==2 && PhotonElectronMatch[iPhoton][2] == PhotonElectronMatch[iPhoton][3]){//merged conv photon
835 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
836 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
837 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
838 if(esdCluster->GetNCells()>=2)
839 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
841 else if(PhotonElectronMatch[iPhoton][0]==2){//non-merged conv photon (but both hit emcal)
842 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
843 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
844 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
845 if(esdCluster->GetNCells()>=2)
846 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
847 esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][3]);
848 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
849 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
850 if(esdCluster->GetNCells()>=2)
851 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
853 else if(PhotonElectronMatch[iPhoton][0]==1){//non-merged conv photon (one missed emcal)
854 AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
855 recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
856 h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
857 if(esdCluster->GetNCells()>=2)
858 h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
860 }//loop over decay photons (iPhoton).
862 }// 2 Photons pointed at the EMCAL!
863 }//for(nTracksMC) ie. Truth Pion loop.
867 //######################### ~~~~~~~~~~~~ ##################################
868 //######################### DONE WITH MC ##################################
869 //######################### ~~~~~~~~~~~~ ##################################
873 // NEW HISTO should be filled before this point, as PostData puts the
874 // information for this iteration of the UserExec in the container
875 PostData(1, fOutput);
878 //________________________________________________________________________
879 void AliAnalysisTaskSDMGammaMC::Terminate(Option_t *) //specify what you want to have done
881 // Called once at the end of the query.
885 //________________________________________________________________________
886 Int_t AliAnalysisTaskSDMGammaMC::GetZvtxBin(Double_t vertZ)
927 //________________________________________________________________________
928 Int_t AliAnalysisTaskSDMGammaMC::GetMultBin(Int_t mult){
942 //________________________________________________________________________
943 Int_t AliAnalysisTaskSDMGammaMC::isGoodEsdCluster(AliESDCaloCluster* esdclust){
948 //double MinErat = 0;
953 if (!esdclust->IsEMCAL())
954 pass = 0;//removes ~70% of clusters.
955 if (esdclust->E()<MinE)
956 pass = 0;//does nothing
957 if (esdclust->GetNCells()<nMinCells)
958 pass = 0;//does nothing
959 //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
961 //if (esdclust->Chi2()<MinEcc) // eccentricity cut
962 //pass = 0;//this is always -1.
965 //This cuts out more than just 1 cell clusters
966 //and drains the statistics badly.
967 //haven't figured out what it does yet.
968 if(esdclust->GetM20()<0.02)
971 //if(esdclust->GetM02()<0.1)
973 //if(esdclust->GetM02()>0.5)
975 //if(esdclust->GetNCells()<2)
978 Float_t pos[3] = {0,0,0};
979 esdclust->GetPosition(pos);
980 TVector3 clusterPosition(pos);
981 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
982 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
984 clusterPosition.Delete();
986 //doing this by hand now...
987 //if(!esdclust->GetNTracksMatched()==0)
993 //________________________________________________________________________
994 Int_t AliAnalysisTaskSDMGammaMC::isGoodAodCluster(AliAODCaloCluster* aodclust){
999 //double MinErat = 0;
1000 //double MinEcc = 0;
1004 if (!aodclust->IsEMCAL())
1005 pass = 0;//removes ~70% of clusters.
1006 if (aodclust->E()<MinE)
1007 pass = 0;//does nothing
1008 if (aodclust->GetNCells()<nMinCells)
1009 pass = 0;//does nothing
1010 //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1012 //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1013 //pass = 0;//this is always -1.
1016 //This cuts out more than just 1 cell clusters
1017 //and drains the statistics badly.
1018 //haven't figured out what it does yet.
1019 if(aodclust->GetM20()<0.02)
1021 if(aodclust->GetM02()<0.02)
1024 //if(aodclust->GetM02()<0.1)
1026 //if(aodclust->GetM02()>0.5)
1028 //if(aodclust->GetNCells()<2)
1031 Float_t pos[3] = {0,0,0};
1032 aodclust->GetPosition(pos);
1033 TVector3 clusterPosition(pos);
1034 if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1035 clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax )
1037 clusterPosition.Delete();
1039 //if(!aodclust->GetNTracksMatched()==0)
1045 //________________________________________________________________________
1046 Double_t AliAnalysisTaskSDMGammaMC::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1048 double dphi = p1.Phi() - p2.Phi();
1050 if(dphi<0.5*TMath::Pi())
1051 dphi = dphi + 2.0*TMath::Pi();
1053 if(dphi>1.5*TMath::Pi())
1054 dphi = dphi - 2.0*TMath::Pi();
1059 //________________________________________________________________________
1060 Double_t AliAnalysisTaskSDMGammaMC::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1062 double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1068 //________________________________________________________________________
1069 Double_t AliAnalysisTaskSDMGammaMC::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1071 double recalibfactor = 0.0;
1073 if(iCalib==0){// no recalibration!
1074 recalibfactor = 1.0;
1076 else if(iCalib==1){// just a scale factor:
1077 recalibfactor = 0.984;
1079 else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1080 Double_t p[3] = {0.96968, -2.68720, -0.831607};
1081 recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1083 else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1084 Double_t p[7] = {1.00000e+00, 3.04925e-02, 4.69043e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.00046e+00};
1085 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1087 else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1088 Double_t p[7] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01};
1089 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1091 else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1092 Double_t p[10] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01, 0.96968, -2.68720, -0.831607};
1093 recalibfactor = ( (p[6]/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) ) / ( p[7] + exp(p[8] + p[9]*energy/2.0) );
1095 else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1096 Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1097 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1099 else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1100 Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1101 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1103 else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1104 Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.964, -3.132, -0.435};
1105 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) * (p[7] + exp(p[8]+p[9]*energy*2.0));
1107 else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1108 Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1109 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1111 else if(iCalib==10){// Jason played with test beam data
1112 Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1113 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1115 else if(iCalib==11){// Jason played with test beam MC
1116 Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1117 recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1120 return recalibfactor;