Passing AliAODEvent to IsSelected (ChiaraB)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSESignificance.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /////////////////////////////////////////////////////////////
17 //
18 // AliAnalysisTaskSESignificane to calculate effects on 
19 // significance of D mesons  cut 
20 // Authors: G. Ortona, ortona@to.infn.it
21 // F. Prino, prino@to.infn.it
22 // Renu Bala, bala@to.infn.it
23 // Chiara Bianchin, cbianchi@pd.infn.it
24 /////////////////////////////////////////////////////////////
25
26 #include <Riostream.h>
27 #include <TClonesArray.h>
28 #include <TCanvas.h>
29 #include <TList.h>
30 #include <TFile.h>
31 #include <TString.h>
32 #include <TH1F.h>
33 #include <TDatabasePDG.h>
34
35 #include <AliLog.h>
36 #include "AliAnalysisManager.h"
37 #include "AliAODHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODVertex.h"
40 #include "AliAODTrack.h"
41 #include "AliAODMCHeader.h"
42 #include "AliAODMCParticle.h"
43 #include "AliAODRecoDecayHF3Prong.h"
44 #include "AliAODRecoDecayHF.h"
45 #include "AliAODRecoDecayHF2Prong.h"
46 #include "AliAODRecoDecayHF4Prong.h"
47 #include "AliAODRecoCascadeHF.h"
48
49 #include "AliAnalysisTaskSE.h"
50 #include "AliRDHFCutsDplustoKpipi.h"
51 #include "AliRDHFCutsD0toKpipipi.h"
52 #include "AliRDHFCutsDstoKKpi.h"
53 #include "AliRDHFCutsDStartoKpipi.h"
54 #include "AliRDHFCutsD0toKpi.h"
55 #include "AliRDHFCutsLctopKpi.h"
56 #include "AliMultiDimVector.h"
57
58 #include "AliAnalysisTaskSESignificance.h"
59
60 ClassImp(AliAnalysisTaskSESignificance)
61
62
63 //________________________________________________________________________
64 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
65   AliAnalysisTaskSE(),
66   fOutput(0),
67   fCutList(0),
68   fHistNEvents(0),
69   fUpmasslimit(1.965),
70   fLowmasslimit(1.765),
71   fRDCuts(0),
72   fNPtBins(0),
73   fReadMC(kFALSE),
74   fDecChannel(0),
75   fSelectionlevel(0)
76 {
77   // Default constructor
78 }
79
80 //________________________________________________________________________
81 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
82   AliAnalysisTaskSE(name),
83   fOutput(0),
84   fCutList(listMDV),
85   fHistNEvents(0),
86   fUpmasslimit(0),
87   fLowmasslimit(0),
88   fRDCuts(rdCuts),
89   fNPtBins(0),
90   fReadMC(kFALSE),
91   fDecChannel(decaychannel),
92   fSelectionlevel(selectionlevel)
93 {
94
95   Int_t pdg=421;
96   switch(fDecChannel){
97   case 0:
98     pdg=411;
99     break;
100   case 1:
101     pdg=421;
102     break;
103   case 2:
104     pdg=413;
105     break;
106   case 3:
107     pdg=431;
108     break;
109   case 4:
110     pdg=421;
111     break;
112   case 5:
113     pdg=4122;
114     break;
115   }
116
117   SetMassLimits(0.1,pdg); //check range
118   fNPtBins=fRDCuts->GetNPtBins();
119
120   if(fDebug>1)fRDCuts->PrintAll();
121    // Output slot #1 writes into a TList container
122   DefineOutput(1,TList::Class());  //My private output
123   DefineOutput(2,TList::Class());
124   CheckConsistency();
125 }
126
127  //________________________________________________________________________
128 AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
129 {
130   // Destructor
131   if (fOutput) {
132     delete fOutput;
133     fOutput = 0;
134   }
135   if (fCutList) {
136     delete fCutList;
137     fCutList = 0;
138   }
139   if(fHistNEvents){
140     delete fHistNEvents;
141     fHistNEvents=0;
142   }
143
144 /*
145   if(fRDCuts) {
146     delete fRDCuts;
147     fRDCuts= 0;
148   } 
149 */
150   
151 }
152 //_________________________________________________________________
153 Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
154
155   Bool_t result = kTRUE;
156
157   const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
158   //Float_t *vars = new Float_t[nvars];
159   Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
160   Bool_t *uppervars = fRDCuts->GetIsUpperCut();
161   TString *names = fRDCuts->GetVarNames();
162
163   for(Int_t i=0;i<fNPtBins;i++){
164     TString mdvname=Form("multiDimVectorPtBin%d",i);
165     Int_t ic=0;
166  
167     for(Int_t ivar=0;ivar<nvars;ivar++){
168       if(varsforopt[ivar]){
169         Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
170         Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
171         if(min==max){
172           printf("AliAnalysisTaskSESignificance::CheckConsistency: ERROR! \n tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i);
173           result = kFALSE;
174         }
175         Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
176         if(uppervars[ivar]&&lowermdv){
177           printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n ",names[ivar].Data(),i);
178           ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
179           result = kTRUE;
180         }
181         if(!uppervars[ivar]&&!lowermdv){
182           printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n",names[ivar].Data(),i);
183           ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
184           result = kTRUE;
185         }
186         ic++;
187       }
188     }
189   }
190   return result;
191 }
192 //_________________________________________________________________
193 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
194   Float_t mass=0;
195   Int_t abspdg=TMath::Abs(pdg);
196   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
197   fUpmasslimit = mass+range;
198   fLowmasslimit = mass-range;
199 }
200 //_________________________________________________________________
201 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
202   if(uplimit>lowlimit)
203     {
204       fUpmasslimit = lowlimit;
205       fLowmasslimit = uplimit;
206     }
207 }
208
209
210
211 //________________________________________________________________________
212 void AliAnalysisTaskSESignificance::LocalInit()
213 {
214   // Initialization
215
216   if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
217
218   TList *mdvList =  new TList();
219   mdvList->SetOwner();
220   mdvList = fCutList;
221   
222   PostData(2,mdvList);
223
224   return;
225 }
226 //________________________________________________________________________
227 void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
228 {
229   // Create the output container
230  
231   if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
232
233   // Several histograms are more conveniently managed in a TList
234   fOutput = new TList();
235   fOutput->SetOwner();
236   fOutput->SetName("OutputHistos");
237
238   //same number of steps in each multiDimVectorPtBin%d !
239   Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
240   cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
241   nHist=nHist*fNPtBins;
242   cout<<"Total = "<<nHist<<endl;
243   for(Int_t i=0;i<nHist;i++){
244
245     TString hisname;
246     TString signame;
247     TString bkgname;
248     TString rflname;
249     TString title;
250     
251     hisname.Form("hMass_%d",i);
252     signame.Form("hSig_%d",i);
253     bkgname.Form("hBkg_%d",i);
254     rflname.Form("hRfl_%d",i);
255     
256     title.Form("Invariant mass;M[GeV/c^{2}];Entries");
257
258     fMassHist[i]=new TH1F(hisname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
259     fMassHist[i]->Sumw2();
260     fOutput->Add(fMassHist[i]);
261
262     if(fReadMC){
263       fSigHist[i]=new TH1F(signame.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
264       fSigHist[i]->Sumw2();
265       fOutput->Add(fSigHist[i]);
266       
267       fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
268       fBkgHist[i]->Sumw2();
269       fOutput->Add(fBkgHist[i]);
270
271       fRflHist[i]=new TH1F(rflname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
272       fRflHist[i]->Sumw2();
273       fOutput->Add(fRflHist[i]);
274
275     }
276   }
277
278   fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",4,0,4.);
279   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
280   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvNotSelected");
281   fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
282   fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
283   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
284   fOutput->Add(fHistNEvents);
285
286  
287   return;
288 }
289
290 //________________________________________________________________________
291 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
292 {
293   // Execute analysis for current event:
294   // heavy flavor candidates association to MC truth
295    
296   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
297   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
298   // Post the data already here
299   PostData(1,fOutput);
300   TClonesArray *arrayProng =0;
301
302   if(!aod && AODEvent() && IsStandardAOD()) { 
303     // In case there is an AOD handler writing a standard AOD, use the AOD 
304     // event in memory rather than the input (ESD) event.    
305     aod = dynamic_cast<AliAODEvent*> (AODEvent());
306     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
307     // have to taken from the AOD event hold by the AliAODExtension
308     AliAODHandler* aodHandler = (AliAODHandler*) 
309       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
310     if(aodHandler->GetExtensions()) {
311       
312       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
313       AliAODEvent *aodFromExt = ext->GetAOD();
314    
315       
316       switch(fDecChannel){
317       case 0:
318         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
319         break; 
320       case 1:
321         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
322         break; 
323       case 2:
324         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
325         break; 
326       case 3:
327         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
328         break; 
329       case 4:
330         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
331         break; 
332       case 5:
333         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
334         break; 
335       }
336     }
337   } else {
338     switch(fDecChannel){
339     case 0:
340       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
341       break; 
342     case 1:
343       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
344       break; 
345     case 2:
346       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
347       break; 
348     case 3:
349       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
350       break; 
351     case 4:
352       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
353       break; 
354     case 5:
355       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
356       break; 
357     }
358   }
359   if(!arrayProng) {
360     printf("AliAnalysisTaskSESignificance::UserExec: branch not found!\n");
361     return;
362   }
363   
364   // fix for temporary bug in ESDfilter 
365   // the AODs with null vertex pointer didn't pass the PhysSel
366   if(!aod->GetPrimaryVertex()) return;
367
368   TClonesArray *arrayMC=0;
369   AliAODMCHeader *mcHeader=0;
370
371   // load MC particles
372   if(fReadMC){
373     
374     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
375     if(!arrayMC) {
376       printf("AliAnalysisTaskSESignificance::UserExec: MC particles branch not found!\n");
377       //    return;
378     }
379     
380     // load MC header
381     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
382     if(!mcHeader) {
383       printf("AliAnalysisTaskSESignificance::UserExec: MC header branch not found!\n");
384       return;
385     }
386   }
387
388   fHistNEvents->Fill(0); // count event
389
390   AliAODRecoDecayHF *d=0;
391   Int_t nprongs=1;
392   Int_t *pdgdaughters=0x0;
393   Int_t prongpdg=211;
394   
395
396   switch(fDecChannel){
397   case 0:
398     //D+
399     pdgdaughters =new Int_t[3];
400     pdgdaughters[0]=211;//pi
401     pdgdaughters[1]=321;//K
402     pdgdaughters[2]=211;//pi
403     nprongs=3;
404     prongpdg=411;
405     break;
406   case 1:
407     //D0
408     pdgdaughters =new Int_t[2];
409     pdgdaughters[0]=211;//pi 
410     pdgdaughters[1]=321;//K
411     nprongs=2;
412     prongpdg=421;
413     break;
414   case 2:
415     //D*
416     pdgdaughters =new Int_t[3];
417     pdgdaughters[1]=211;//pi
418     pdgdaughters[0]=321;//K
419     pdgdaughters[2]=211;//pi (soft?)
420     nprongs=3;
421     prongpdg=413;
422     break;
423   case 3:
424     //Ds
425     pdgdaughters =new Int_t[3];
426     pdgdaughters[0]=211;//pi
427     pdgdaughters[1]=321;//K
428     pdgdaughters[2]=321;//K
429     nprongs=3;
430     prongpdg=431;
431     break;
432   case 4:
433     //D0
434     pdgdaughters =new Int_t[4];
435     pdgdaughters[0]=321;
436     pdgdaughters[1]=211;
437     pdgdaughters[2]=211;
438     pdgdaughters[3]=211;
439     nprongs=4;
440     prongpdg=421;
441     break;
442   case 5:
443     //Lambda_c
444     pdgdaughters =new Int_t[3];
445     //check order
446     pdgdaughters[0]=2212;//p
447     pdgdaughters[1]=321;//K
448     pdgdaughters[2]=211;//pi
449     nprongs=3;
450     prongpdg=4122;
451     fReadMC=kFALSE;
452     printf("Warning! For the moment it is not allowed to use MC for the Lambdac\n");
453     break;
454   }
455   Int_t prongPdgPlus=prongpdg;
456   Int_t prongPdgMinus=TMath::Abs(prongPdgPlus)*(-1);
457
458   Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
459   Int_t nProng = arrayProng->GetEntriesFast();
460   if(fDebug>1) printf("Number of D2H: %d\n",nProng);
461
462   if(!fRDCuts->IsEventSelected(aod)){
463     fHistNEvents->Fill(1);
464     return;
465   }
466
467   for (Int_t iProng = 0; iProng < nProng; iProng++) {
468     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
469     Double_t invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
470     Double_t invMassC=0;
471     
472     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
473     
474     if(isSelected) {
475       fHistNEvents->Fill(2); // count selected with loosest cuts
476       if(fDebug>1) printf("+++++++Is Selected\n");
477     
478       const Int_t nvars=fRDCuts->GetNVarsForOpt();
479       Float_t *vars = new Float_t[nvars];
480       Int_t nVals=0;
481       
482       fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
483       Int_t ptbin=fRDCuts->PtBin(d->Pt());
484       if(ptbin==-1) continue;
485       TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
486       ULong64_t *addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
487       if(fDebug>1)printf("nvals = %d\n",nVals);
488       for(Int_t ivals=0;ivals<nVals;ivals++){
489         if(addresses[ivals]>=((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetNTotCells()){
490           if (fDebug>1) printf("Overflow!!\n");
491           return;
492         }
493         fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
494         fHistNEvents->Fill(3);
495         if(fReadMC){
496           Int_t lab=-1;
497           lab = d->MatchToMC(prongPdgPlus,arrayMC,nprongs,pdgdaughters);
498           if(lab>=0){ //signal
499             AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
500             Int_t pdgMC = dMC->GetPdgCode();
501             
502             if(pdgMC==prongPdgPlus) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
503             else if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi) fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass); //temporary! an implementation of the reflections for 3 prongs is in progress
504           }
505           else{ //background
506             fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
507           }
508         
509         }
510       }
511
512
513       if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi || fDecChannel==AliAnalysisTaskSESignificance::kLambdactopKpi){ //repeat the cycle to checks cuts passed as D0bar
514         if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi){
515           pdgdaughters[0]=321;
516           pdgdaughters[1]=211; 
517         }
518         if(fDecChannel==AliAnalysisTaskSESignificance::kLambdactopKpi){
519           pdgdaughters[0]=211;
520           pdgdaughters[1]=321;
521           pdgdaughters[2]=2212;
522         }
523
524         invMassC=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
525         fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
526         addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
527         for(Int_t ivals=0;ivals<nVals;ivals++){
528           fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
529           if(fReadMC){
530             Int_t lab=-1;
531             lab = d->MatchToMC(prongPdgMinus,arrayMC,nprongs,pdgdaughters);
532             if(lab>=0){ //signal
533               AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
534               Int_t pdgMC = dMC->GetPdgCode();
535               if(pdgMC==prongPdgMinus) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
536               else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
537             }
538             else{ //background
539               fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
540             }
541             
542           }
543         }
544       }
545   
546     }// end if selected
547
548     
549   }
550
551   PostData(1,fOutput);    
552   return;
553 }
554
555
556 //________________________________________________________________________
557 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
558 {
559   // Terminate analysis
560   //
561   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
562
563
564   fOutput = dynamic_cast<TList*> (GetOutputData(1));
565   if (!fOutput) {     
566     printf("ERROR: fOutput not available\n");
567     return;
568   }
569  
570   fCutList = dynamic_cast<TList*> (GetOutputData(2));
571   if (!fCutList) {     
572     printf("ERROR: fCutList not available\n");
573     return;
574   }
575
576   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
577   if (!mdvtmp){
578     cout<<"multidimvec not found in TList"<<endl;
579     fCutList->ls();
580     return;
581   }
582   Int_t nHist=mdvtmp->GetNTotCells();
583   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
584   Bool_t drawn=kFALSE;
585   for(Int_t i=0;i<nHist;i++){
586
587     TString hisname;
588     TString signame;
589     TString bkgname;
590     TString rflname;
591     
592     hisname.Form("hMass_%d",i);
593     signame.Form("hSig_%d",i);
594     bkgname.Form("hBkg_%d",i);
595     rflname.Form("hRfl_%d",i);
596
597     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
598
599     if (!drawn && fMassHist[i]->GetEntries() > 0){
600       c1->cd();
601       fMassHist[i]->Draw();
602       drawn=kTRUE;
603     }
604     
605     if(fReadMC){
606       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
607       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
608       fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
609     }
610     
611   }
612   
613
614  
615   
616   return;
617 }
618 //-------------------------------------------
619