fix compilation warnings, coverity and add new functionalities (Laurent)
authorpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Mar 2013 12:51:42 +0000 (12:51 +0000)
committerpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Mar 2013 12:51:42 +0000 (12:51 +0000)
15 files changed:
PWG/muon/AliAnalysisMuMuBinning.cxx
PWG/muon/AliAnalysisMuMuBinning.h
PWG/muon/AliAnalysisTaskMuMu.cxx
PWG/muon/AliAnalysisTaskMuMu.h
PWG/muondep/AccEffTemplates/Config.C
PWG/muondep/AccEffTemplates/rec.C
PWG/muondep/AliAnalysisMuMu.cxx
PWG/muondep/AliAnalysisMuMu.h
PWG/muondep/AliAnalysisMuMuResult.cxx
PWG/muondep/AliAnalysisMuMuResult.h
PWG/muondep/AliAnalysisMuMuSpectra.cxx
PWG/muondep/AliAnalysisMuMuSpectra.h
PWG/muondep/AliAnalysisTriggerScalers.cxx
PWG/muondep/AliMuonAccEffSubmitter.cxx
PWG/muondep/AliMuonAccEffSubmitter.h

index b26debc..00451fe 100644 (file)
@@ -16,6 +16,7 @@
 #include "TMath.h"
 #include "TObjArray.h"
 #include "TObjString.h"
+#include <cassert>
 
 ClassImp(AliAnalysisMuMuBinning::Range)
 ClassImp(AliAnalysisMuMuBinning)
@@ -76,14 +77,14 @@ void AliAnalysisMuMuBinning::AddBin(const AliAnalysisMuMuBinning::Range& bin)
   /// add one bin
   AddBin(bin.Particle().Data(),bin.Type().Data(),
          bin.Xmin(),bin.Xmax(),
-         bin.Ymin(),bin.Ymax());
+         bin.Ymin(),bin.Ymax(),bin.Flavour());
 }
-                                    
+
 //______________________________________________________________________________
 void AliAnalysisMuMuBinning::AddBin(const char* particle, const char* type,
                                     Double_t xmin, Double_t xmax,
                                     Double_t ymin, Double_t ymax,
-                                    Bool_t warn)
+                                    const char* flavour)
 {
   /// Add a bin
   /// Note that particle and type are not case sensitive.
@@ -107,14 +108,11 @@ void AliAnalysisMuMuBinning::AddBin(const char* particle, const char* type,
     fBins->Add(new TObjString(sparticle),b);
   }
 
-  Range* r = new Range(sparticle.Data(),stype.Data(),xmin,xmax,ymin,ymax);
+  Range* r = new Range(sparticle.Data(),stype.Data(),xmin,xmax,ymin,ymax,flavour);
   
   if ( b->FindObject(r) )
   {
-    if (warn)
-    {
-      AliWarning("Trying to add an already existing bin. Not doing it.");      
-    }
+    AliDebug(1,"Trying to add an already existing bin. Not doing it.");
     delete r;
   }
   else
@@ -234,7 +232,7 @@ TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle) const
 
 
 //______________________________________________________________________________
-TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle, const char* type) const
+TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle, const char* type, const char* flavour) const
 {
   /// Get the list of bins for a given particle and given type
   /// The returned array must be deleted by the user
@@ -255,6 +253,8 @@ TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle, const
   TString stype(type);
   stype.ToUpper();
 
+  TString sflavour(flavour);
+  
   TObjArray* types = stype.Tokenize(",");
   TObjString* onetype;
   TIter nextType(types);
@@ -264,7 +264,8 @@ TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle, const
     nextType.Reset();
     while ( ( onetype = static_cast<TObjString*>(nextType()) ) )
     {
-      if ( r->Type() == onetype->String() )
+      if ( r->Type() == onetype->String() &&
+          ( ( sflavour.Length() > 0 && r->Flavour() == sflavour.Data() ) || sflavour.Length()==0 ) )
       {
         a->Add(r->Clone());
       }
@@ -284,23 +285,24 @@ TObjArray* AliAnalysisMuMuBinning::CreateBinObjArray(const char* particle, const
 //______________________________________________________________________________
 void AliAnalysisMuMuBinning::CreateMesh(const char* particle,
                                         const char* type1, const char* type2,
+                                        const char* flavour,
                                         Bool_t remove12)
 {
   /// Create 2D bins from existing 1d ones of type1 and type2
-  TObjArray* a1 = CreateBinObjArray(particle,type1);
+  TObjArray* a1 = CreateBinObjArray(particle,type1,flavour);
   if (!a1)
   {
     AliError(Form("No bin for type %s. Done nothing.",type1));
     return;
   }
-  TObjArray* a2 = CreateBinObjArray(particle,type2);
+  TObjArray* a2 = CreateBinObjArray(particle,type2,flavour);
   if (!a2)
   {
     AliError(Form("No bin for type %s. Done nothing.",type2));
     return;
   }
   
-  TString meshType(Form("%s VS %s",type1,type2));
+  TString meshType(Form("%s VS %s - %s",type1,type2,flavour));
   
   for ( Int_t i1 = 0; i1 <= a1->GetLast(); ++i1 )
   {
@@ -310,7 +312,7 @@ void AliAnalysisMuMuBinning::CreateMesh(const char* particle,
     {
       Range* r2 = static_cast<Range*>(a2->At(i2));
       
-      AddBin(particle,meshType,r2->Xmin(),r2->Xmax(),r1->Xmin(),r1->Xmax(),kFALSE);
+      AddBin(particle,meshType,r2->Xmin(),r2->Xmax(),r1->Xmin(),r1->Xmax(),Form("%s VS %s",r1->Flavour().Data(),r2->Flavour().Data()));
     }
   }
   
@@ -487,11 +489,11 @@ Long64_t AliAnalysisMuMuBinning::Merge(TCollection* list)
 
 //______________________________________________________________________________
 AliAnalysisMuMuBinning*
-AliAnalysisMuMuBinning::Project(const char* particle, const char* type) const
+AliAnalysisMuMuBinning::Project(const char* particle, const char* type, const char* flavour) const
 {
   /// Create a sub-binning object with only the bins pertaining to (particle,type)
   
-  TObjArray* bins = CreateBinObjArray(particle,type);
+  TObjArray* bins = CreateBinObjArray(particle,type,flavour);
   if (!bins) return 0x0;
   AliAnalysisMuMuBinning* p = new AliAnalysisMuMuBinning;
   TIter next(bins);
@@ -501,9 +503,9 @@ AliAnalysisMuMuBinning::Project(const char* particle, const char* type) const
   
   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
   {
-    if  (bin->Type()==stype)
+    assert  (bin->Type()==stype && bin->Flavour()==flavour);
     {
-      p->AddBin(particle,bin->Type(),bin->Xmin(),bin->Xmax(),bin->Ymin(),bin->Ymax());
+      p->AddBin(particle,bin->Type(),bin->Xmin(),bin->Xmax(),bin->Ymin(),bin->Ymax(),bin->Flavour().Data());
     }
   }
   
@@ -550,8 +552,11 @@ void AliAnalysisMuMuBinning::Print(Option_t* /*opt*/) const
 //______________________________________________________________________________
 AliAnalysisMuMuBinning::Range::Range(const char* particle, const char* type,
                                      Double_t xmin, Double_t xmax,
-                                     Double_t ymin, Double_t ymax)
-: TObject(), fParticle(particle), fType(type), fXmin(xmin), fXmax(xmax), fYmin(ymin), fYmax(ymax)
+                                     Double_t ymin, Double_t ymax,
+                                     const char* flavour)
+: TObject(), fParticle(particle), fType(type),
+  fXmin(xmin), fXmax(xmax), fYmin(ymin), fYmax(ymax),
+  fFlavour(flavour)
 {
   /// ctor
   fParticle.ToUpper();
@@ -567,7 +572,14 @@ TString AliAnalysisMuMuBinning::Range::AsString() const
   
   TString s;
   
-  s.Form("%s_%05.2f_%05.2f",Type().Data(),Xmin(),Xmax());
+  if ( fFlavour.Length() > 0 )
+  {
+    s.Form("%s_%s_%05.2f_%05.2f",Type().Data(),Flavour().Data(),Xmin(),Xmax());
+  }
+  else
+  {
+    s.Form("%s_%05.2f_%05.2f",Type().Data(),Xmin(),Xmax());
+  }
   
   if (Is2D())
   {
@@ -596,6 +608,12 @@ Int_t      AliAnalysisMuMuBinning::Range::Compare(const TObject* obj) const
   
   if (s) return s;
   
+  s = strcmp(Flavour().Data(),other->Flavour().Data());
+  
+  if (s) return s;
+  
+  if ( IsNullObject() ) return 0;
+  
   if ( Xmin() < other->Xmin() )
   {
     return -1;
@@ -693,6 +711,11 @@ void AliAnalysisMuMuBinning::Range::Print(Option_t* /*opt*/) const
     std::cout << Form(" ; %5.2f : %5.2f",Ymin(),Ymax());
   }
   
+  if (Flavour().Length()>0)
+  {
+    std::cout << " - " << Flavour().Data();
+  }
+  
   std::cout << "->" << AsString().Data() << std::endl;
   
 }
index fbc1f54..ccf111b 100644 (file)
@@ -35,7 +35,7 @@ public:
               Double_t xmax=TMath::Limits<Double_t>::Max(),
               Double_t ymin=TMath::Limits<Double_t>::Max(),
               Double_t ymax=TMath::Limits<Double_t>::Max(),
-              Bool_t warn=kTRUE);
+              const char* flavour="");
 
   TObjArray* CreateParticleArray() const;
 
@@ -45,13 +45,13 @@ public:
 
   TObjArray* CreateBinObjArray() const;
   TObjArray* CreateBinObjArray(const char* particle) const;
-  TObjArray* CreateBinObjArray(const char* particle, const char* type) const;
+  TObjArray* CreateBinObjArray(const char* particle, const char* type, const char* flavour) const;
   
-  AliAnalysisMuMuBinning* Project(const char* particle, const char* type) const;
+  AliAnalysisMuMuBinning* Project(const char* particle, const char* type, const char* flavour="") const;
   
   virtual void Print(Option_t* opt="") const;
   
-  void CreateMesh(const char* particle, const char* type1, const char* type2, Bool_t remove12=kFALSE);
+  void CreateMesh(const char* particle, const char* type1, const char* type2, const char* flavour="", Bool_t remove12=kFALSE);
 
   Long64_t Merge(TCollection* list);
 
@@ -65,7 +65,8 @@ public:
           Double_t xmin=TMath::Limits<Double_t>::Max(),
           Double_t xmax=TMath::Limits<Double_t>::Max(),
           Double_t ymin=TMath::Limits<Double_t>::Max(),
-          Double_t ymax=TMath::Limits<Double_t>::Max());
+          Double_t ymax=TMath::Limits<Double_t>::Max(),
+          const char* version="");
     
     virtual Int_t      Compare(const TObject* obj) const;
     Bool_t IsEqual(const TObject* obj) const { return Compare(obj)==0; }
@@ -100,6 +101,8 @@ public:
     
     Bool_t IsInRange(Double_t x, Double_t y=TMath::Limits<Double_t>::Max()) const;
     
+    TString Flavour() const { return fFlavour; }
+    
   private:
     TString fParticle; // particle
     TString fType; // binning type (e.g. pt, y, phi)
@@ -107,8 +110,9 @@ public:
     Double_t fXmax; // x-max of the range
     Double_t fYmin; // x-min of the range
     Double_t fYmax; // x-max of the range
+    TString fFlavour; // flavour (if any) this range, e.g. coarse, fine, etc...
     
-    ClassDef(AliAnalysisMuMuBinning::Range,1)
+    ClassDef(AliAnalysisMuMuBinning::Range,2)
   };
   
 
@@ -116,7 +120,7 @@ public:
 
   TMap* fBins; // list of bins (particle -> list of bins) = (TObjString -> TObjArray)
   
-  ClassDef(AliAnalysisMuMuBinning,1)
+  ClassDef(AliAnalysisMuMuBinning,1) // custom binning for MuMu analysis
 };
 
 #endif
index 6af9405..67ed3de 100644 (file)
@@ -257,10 +257,11 @@ AliAnalysisTaskMuMu::~AliAnalysisTaskMuMu()
 //_____________________________________________________________________________
 void AliAnalysisTaskMuMu::AddBin(const char* particle, const char* type,
                                  Double_t xmin, Double_t xmax,
-                                 Double_t ymin, Double_t ymax)
+                                 Double_t ymin, Double_t ymax,
+                                 const char* flavour)
 {
   /// Add one bin
-  fBinning->AddBin(particle,type,xmin,xmax,ymin,ymax);
+  fBinning->AddBin(particle,type,xmin,xmax,ymin,ymax,flavour);
 
   // invalidate cached bins, if any
   if (fBinArray)
@@ -414,7 +415,7 @@ void AliAnalysisTaskMuMu::ComputeTrackMask(const AliVParticle& track, Int_t trac
 }
 
 //_____________________________________________________________________________
-void AliAnalysisTaskMuMu::CreateMesh(const char* particle, const char* type1, const char* type2, Bool_t remove12)
+void AliAnalysisTaskMuMu::CreateMesh(const char* particle, const char* type1, const char* type2, const char* flavour, Bool_t remove12)
 {
   /// Create a 2d binning from 2 1d binning.
   /// WARNING : not fully tested yet
@@ -424,7 +425,7 @@ void AliAnalysisTaskMuMu::CreateMesh(const char* particle, const char* type1, co
     AliError("Cannot create a mesh as I have no bin at all !");
     return;
   }
-  fBinning->CreateMesh(particle,type1,type2,remove12);
+  fBinning->CreateMesh(particle,type1,type2,flavour,remove12);
 }
 
 //_____________________________________________________________________________
@@ -440,7 +441,7 @@ AliAnalysisTaskMuMu::CreateMinvHistograms(const char* physics,
 
   Int_t nMCMinvBins = GetNbins(minvMin,minvMax,0.1);
 
-  TObjArray* bins = fBinning->CreateBinObjArray("psi","pt vs y,pt,y,phi");
+  TObjArray* bins = fBinning->CreateBinObjArray("psi","y vs pt,integrated,pt,y,phi","");
   
   CreatePairHisto(physics,triggerClassName,"Pt","#mu+#mu- Pt distribution",
                   200,0,20);
@@ -653,7 +654,7 @@ void AliAnalysisTaskMuMu::DefineCentralityClasses(TArrayF* centralities)
 void AliAnalysisTaskMuMu::DefineDefaultBinning()
 {
   fBinning = new AliAnalysisMuMuBinning("BIN");
-  fBinning->AddBin("psi","pt vs y");
+  fBinning->AddBin("psi","integrated");
 }
 
 //_____________________________________________________________________________
@@ -863,7 +864,7 @@ void AliAnalysisTaskMuMu::FillMC()
 
   if (!fBinArray)
   {
-    fBinArray = fBinning->CreateBinObjArray("psi","pt vs y,pt,y");
+    fBinArray = fBinning->CreateBinObjArray("psi","y vs pt,integrated,pt,y","");
   }
 
   TIter nextBin(fBinArray);
@@ -898,9 +899,24 @@ void AliAnalysisTaskMuMu::FillMC()
       {
         Bool_t ok(kFALSE);
         
-        if ( r->Is2D() || r->IsNullObject() )
+        if ( r->IsNullObject() )
         {
-          ok = r->IsInRange(part->Y(),part->Pt());
+          ok = kTRUE;
+        }
+        else if ( r->Is2D() )
+        {
+          if ( r->AsString().BeginsWith("PTVSY") )
+          {
+            ok = r->IsInRange(part->Y(),part->Pt());
+          }
+          else if ( r->AsString().BeginsWith("YVSPT") )
+          {
+            ok = r->IsInRange(part->Pt(),part->Y());
+          }
+          else
+          {
+            AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
+          }
         }
         else
         {
@@ -912,6 +928,10 @@ void AliAnalysisTaskMuMu::FillMC()
           {
             ok = r->IsInRange(part->Y());
           }
+          else if ( r->Type() == "PHI" )
+          {
+            ok = r->IsInRange(part->Phi());
+          }
         }
         
         if ( ok )
@@ -979,6 +999,8 @@ void AliAnalysisTaskMuMu::FillEventHistos(const char* physics, const char* trigg
 {
   // Fill event-wise histograms
   
+//  AliCodeTimerAuto("",0);
+  
   if (!IsHistogramDisabled("BCX"))
   {
     Histo(physics,triggerClassName,centrality,"BCX")->Fill(1.0*Event()->GetBunchCrossNumber());
@@ -1094,6 +1116,7 @@ void AliAnalysisTaskMuMu::FillEventHistos(const char* physics, const char* trigg
 
 }
 
+
 //_____________________________________________________________________________
 void AliAnalysisTaskMuMu::FillHistogramCollection(const char* physics, const char* triggerClassName)
 {
@@ -1243,6 +1266,8 @@ void AliAnalysisTaskMuMu::FillHistosForTrack(const char* physics,
 {
   /// Fill histograms for one track
   
+//  AliCodeTimerAuto("",0);
+  
   TLorentzVector p(track.Px(),track.Py(),track.Pz(),
                    TMath::Sqrt(MuonMass2()+track.P()*track.P()));
   
@@ -1353,13 +1378,15 @@ void AliAnalysisTaskMuMu::FillHistos(const char* physics, const char* triggerCla
 {
   /// Fill histograms for /physics/triggerClassName/centrality
   
+//  AliCodeTimerAuto("",0);
+
   FillEventHistos(physics,triggerClassName,centrality);
   
   // Track loop
   
   if (!fBinArray)
   {
-    fBinArray = fBinning->CreateBinObjArray("psi","pt vs y,pt,y,phi");
+    fBinArray = fBinning->CreateBinObjArray("psi","y vs pt,integrated,pt,y,phi","");
   }
   
   Int_t nMuonTracks = AliAnalysisMuonUtility::GetNTracks(Event());
@@ -1456,14 +1483,34 @@ void AliAnalysisTaskMuMu::FillHistos(const char* physics, const char* triggerCla
             
             while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
             {
-              AliDebug(2,Form("bin %s pt %e rad %e = %d",r->AsString().Data(),
-                              pj.Pt(),pj.Rapidity(),r->IsInRange(pj.Pt(),pj.Rapidity())));
+//              AliInfo(Form("%s y %e pt %e ok1 %d ok2 %d",
+//                           r->AsString().Data(),
+//                           pj.Rapidity(),
+//                           pj.Pt(),
+//                           r->IsInRange(pj.Rapidity(),pj.Pt()),
+//                           r->IsInRange(pj.Pt(),pj.Rapidity())));
               
+                           
               Bool_t ok(kFALSE);
               
-              if ( r->Is2D() || r->IsNullObject() )
+              if ( r->IsNullObject() )
               {
-                ok = r->IsInRange(pj.Rapidity(),pj.Pt());
+                ok = kTRUE;
+              }
+              else if ( r->Is2D() )
+              {
+                if ( r->AsString().BeginsWith("PTVSY") )
+                {
+                  ok = r->IsInRange(pj.Rapidity(),pj.Pt());                  
+                }
+                else if ( r->AsString().BeginsWith("YVSPT") )
+                {
+                  ok = r->IsInRange(pj.Pt(),pj.Rapidity());
+                }
+                else
+                {
+                  AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
+                }
               }
               else
               {
@@ -1493,10 +1540,6 @@ void AliAnalysisTaskMuMu::FillHistos(const char* physics, const char* triggerCla
                   {
                     AliError(Form("Could not get %s",hname.Data()));
                   }
-                  else
-                  {
-                    AliDebug(1,Form("filling %s",hname.Data()));
-                  }
                   h->Fill(pj.M());
                 }
               }
@@ -1542,6 +1585,8 @@ UInt_t AliAnalysisTaskMuMu::GetEventMask() const
    
    */
   
+//  AliCodeTimerAuto("",0);
+  
   UInt_t eMask = EventCuts()->GetSelectionMask(fInputHandler);
 
   UInt_t m(AliAnalysisTaskMuMu::kEventAll);
@@ -1620,21 +1665,21 @@ UInt_t AliAnalysisTaskMuMu::GetEventMask() const
     m |= AliAnalysisTaskMuMu::kEventZ7;
   }
   
-  AliVVZERO* vzero = Event()->GetVZEROData();
-  
-  if (vzero)
-  {
-    Float_t v0a = vzero->GetV0ATime();
-    Float_t v0c = vzero->GetV0CTime();
-    
-    Float_t x = v0a-v0c;
-    Float_t y = v0a+v0c;
-    
-    if ( ( x > 6 && x < 10 ) && y > 20 )
-    {
-      m |= AliAnalysisTaskMuMu::kEventV0UP;
-    }
-  }
+//  AliVVZERO* vzero = Event()->GetVZEROData();
+//  
+//  if (vzero)
+//  {
+//    Float_t v0a = vzero->GetV0ATime();
+//    Float_t v0c = vzero->GetV0CTime();
+//    
+//    Float_t x = v0a-v0c;
+//    Float_t y = v0a+v0c;
+//    
+//    if ( ( x > 6 && x < 10 ) && y > 20 )
+//    {
+//      m |= AliAnalysisTaskMuMu::kEventV0UP;
+//    }
+//  }
   
   Bool_t backgroundFlag(kFALSE);
   Bool_t pileupFlag(kFALSE);
@@ -2041,23 +2086,22 @@ AliAnalysisTaskMuMu::Print(Option_t* /*opt*/) const
       str->Print();
     }
   }
+
+  if ( !fEventCutNames )
+  {
+    cout << "No event cut defined yet" << endl;
+  }
+  else
+  {
+    TIter nextEventCutName(fEventCutNames);
+    TObjString* eventCutName;
+    
+    while ( ( eventCutName = static_cast<TObjString*>(nextEventCutName()) ) )
+    {
+      cout << Form("EVENT  CUT %s MASK %x",eventCutName->String().Data(),eventCutName->GetUniqueID()) << endl;
+    }
+  }
   
-//  if ( !fTriggerClasses || !fTriggerClasses->First() ) 
-//  {
-//    cout << "No trigger classes defined yet" << endl;
-//  }
-//  else
-//  {
-//    cout << "Trigger classes that will be considered:" << endl;
-//    TIter next(fTriggerClasses);
-//    TObjString* s;
-//    
-//    while ( ( s = static_cast<TObjString*>(next()) ) )
-//    {
-//      cout << s->String().Data() << endl;
-//    }
-//  }
-//  
   if ( fBinning )
   {
     cout << "Binning for Minv plots" << endl;
@@ -2131,6 +2175,8 @@ void AliAnalysisTaskMuMu::UserExec(Option_t* /*opt*/)
 {
   /// Executed at each event
   
+//  AliCodeTimerAuto("",0);
+  
   fHasMC = (MCEvent()!=0x0);
 
   TString firedTriggerClasses(AliAnalysisMuonUtility::GetFiredTriggerClasses(Event()));
@@ -2183,8 +2229,6 @@ void AliAnalysisTaskMuMu::UserExec(Option_t* /*opt*/)
         fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), "HASMC", fCurrentRunNumber));
       }
       
-//      fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), eventType.Data(), fCurrentRunNumber));
-
       if ( firedTriggerClasses == "" )
       {
         fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), "EMPTY", fCurrentRunNumber));
@@ -2208,22 +2252,11 @@ void AliAnalysisTaskMuMu::UserExec(Option_t* /*opt*/)
       {
         fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), "ATLEASTONEMBTRIGGER", fCurrentRunNumber));
       }
-
-      if ( AtLeastOneEmcalTrigger(firedTriggerClasses) )
-      {
-        fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), "ATLEASTONEEMCALTRIGGER", fCurrentRunNumber));
-
-        if ( AtLeastOneMBTrigger(firedTriggerClasses) )
-        {
-          fEventCounters->Count(Form("event:%s/trigger:%s/run:%d", et->String().Data(), "ATLEASTONEEMCALORMBTRIGGER", fCurrentRunNumber));
-        }
-      }
-
     }
   }
 
   // second loop to count only the triggers we're interested in
-  
+
   TIter next(EventCuts()->GetSelectedTrigClassesInEvent(Event()));
   TObjString* tname;
   
index a5be6be..31aef24 100644 (file)
 #  include "TArrayI.h"
 #endif
 
+#ifndef ROOT_TMath
+#  include "TMath.h"
+#endif
+
 class AliAnalysisMuMuBinning;
 class AliCounterCollection;
 class AliMergeableCollection;
@@ -76,12 +80,18 @@ public:
   AliAnalysisTaskMuMu(Bool_t fromESD, TList* triggerClassesToConsider, const char* beamYear=0x0, TArrayF* centralities=0x0);
   AliAnalysisTaskMuMu(Bool_t fromESD, const char* beamYear, TArrayF* centralities=0x0);
   virtual ~AliAnalysisTaskMuMu();
-  
+
   void AddBin(const char* particle, const char* type,
               Double_t xmin, Double_t xmax,
-              Double_t ymin=0.0, Double_t ymax=0.0);
+              Double_t ymin,
+              Double_t ymax,
+              const char* flavour="");
 
-  void CreateMesh(const char* particle, const char* type1, const char* type2, Bool_t remove12=kFALSE);
+  void AddBin(const char* particle, const char* type,
+              Double_t xmin, Double_t xmax,
+              const char* flavour="") { AddBin(particle,type,xmin,xmax,TMath::Limits<Double_t>::Max(),TMath::Limits<Double_t>::Max(),flavour); }
+  
+  void CreateMesh(const char* particle, const char* type1, const char* type2, const char* flavour="", Bool_t remove12=kFALSE);
   
   virtual void AddEventCut(const char* cutName, UInt_t mask);
   
index b9e0351..70e66e2 100644 (file)
@@ -47,7 +47,6 @@
 
 //--- Functions ---
 class AliGenPythia;
-AliGenerator* CreateGenerator();
 
 void Config()
 {
@@ -144,8 +143,15 @@ void Config()
   //=========================//
   // Generator Configuration //
   //=========================//
-  AliGenerator* gener = CreateGenerator();
 
+  //AliGenerator* gener = CreateGenerator();
+
+  std::cout << "VAR_GENERATOR settings " << std::endl;
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+  gSystem->AddIncludePath("-I$ALICE_ROOT/EVGEN");
+  gROOT->LoadMacro("VAR_GENERATOR.C+");
+  AliGenerator* gener = VAR_GENERATOR();
+  
   gener->SetOrigin(0., 0., 0.); // Taken from OCDB
   gener->SetSigma(0., 0., 0.);      // Sigma in (X,Y,Z) (cm) on IP position, sigmaz taken from OCDB
 //  gener->SetVertexSmear(kPerEvent);
@@ -267,34 +273,3 @@ void Config()
         AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
     }
 }
-
-//
-//           generators
-//
-AliGenerator* CreateGenerator()
-{
-  AliGenMC* generator = 0x0;
-  TString genType("VAR_GEN_TYPE");
-  genType.ToUpper();
-  if ( genType == "GENLIB" ) {
-    generator = new AliGenParam(1, VAR_GENLIB_TYPE,VAR_GENLIB_PARNAME);
-    generator->SetMomentumRange(0,999);
-    generator->SetPtRange(0,50.);
-    generator->SetYRange(-4.1,-2.4);
-    generator->SetPhiRange(0., 360.);
-    generator->SetChildThetaRange(0.,180.);
-    generator->SetForceDecay(kDiMuon);
-  }
-  else if ( genType == "GENCORR" ) {
-    generator = new AliGenCorrHF(1, VAR_GENCORR_QUARK, VAR_GENCORR_ENERGY);
-    generator->SetMomentumRange(0,9999);
-    generator->SetChildThetaRange(160.0,180.0);
-    generator->SetForceDecay(kSemiMuonic);
-  }
-  
-  generator->SetCutOnChild(1);
-  generator->SetChildPhiRange(0.,360.);
-  generator->SetTrackingFlag(1);
-  
-  return generator;
-}
index b9e4c62..a5399b7 100644 (file)
@@ -30,5 +30,5 @@ void rec(int dset=0)
   // MUON Tracker Residual Alignment
   reco.SetSpecificStorage("MUON/Align/Data","alien://folder=/alice/simulation/2008/v4-15-Release/Residual");
 
-     reco.Run();
+  reco.Run();
 }
index 7ebbd7f..8cf1234 100644 (file)
@@ -50,7 +50,7 @@
 #include "TParameter.h"
 #include "TPaveText.h"
 #include "TROOT.h"
-#include "TStyle.h"
+#include "TStopwatch.h"
 #include "TStyle.h"
 #include "TSystem.h"
 #include <cassert>
@@ -62,21 +62,49 @@ ClassImp(AliAnalysisMuMu)
 
 TString AliAnalysisMuMu::fgOCDBPath("raw://");
 
-TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-B-NOPF-MUON,CMUL7-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-MUON,CMUL8-S-NOPF-MUON,CMUL7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-MUON,CPBI1MUL-B-NOPF-MUON,CMULLO-B-NOPF-MUON");
+TString AliAnalysisMuMu::fgDefaultDimuonTriggers("CMUL7-B-NOPF-MUON");
+
+//,CMUL7-S-NOPF-ALLNOTRD,CMUL7-S-NOPF-MUON,CMUL8-S-NOPF-MUON,CMUL7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-ALLNOTRD,CMUU7-B-NOPF-MUON,CPBI1MUL-B-NOPF-MUON,CMULLO-B-NOPF-MUON");
+
+TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON");
 
-TString AliAnalysisMuMu::fgDefaultMuonTriggers("CMSL7-S-NOPF-MUON,CMSL7-S-NOPF-ALLNOTRD,CMSL8-S-NOPF-MUON,CMSL8-S-NOPF-ALLNOTRD,CMSL7-B-NOPF-MUON,CMUS1-B-NOPF-MUON,CMUS7-B-NOPF-MUON,CMSNGL-B-NOPF-MUON");
+//,CMSL7-S-NOPF-ALLNOTRD,CMSL8-S-NOPF-MUON,CMSL8-S-NOPF-ALLNOTRD,CMSL7-B-NOPF-MUON,CMUS1-B-NOPF-MUON,CMUS7-B-NOPF-MUON,CMSNGL-B-NOPF-MUON");
 
-TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
+TString AliAnalysisMuMu::fgDefaultMinbiasTriggers("CINT7-B-NOPF-ALLNOTRD");
 
-TString AliAnalysisMuMu::fgDefaultEventSelectionList("PSALL");
+//,CINT7-S-NOPF-ALLNOTRD,CINT8-B-NOPF-ALLNOTRD,CINT8-S-NOPF-ALLNOTRD,CINT1-B-NOPF-ALLNOTRD,CPBI2_B1-B-NOPF-ALLNOTRD");
+
+TString AliAnalysisMuMu::fgDefaultEventSelectionList("PSALL"); // for real data, for simulation see AliAnalysisMuMu ctor
 
 TString AliAnalysisMuMu::fgDefaultPairSelectionList("pMATCHLOWRABSBOTH");
 
 TString AliAnalysisMuMu::fgDefaultCentralitySelectionList("PP");
 
+//TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOW:2,PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
+TString AliAnalysisMuMu::fgDefaultFitTypeList("PSILOWalphaLow0.984nLow5.839alphaUp1.972nUp3.444:2,PSILOWMCTAILS:2");
+
+TString AliAnalysisMuMu::fgDefaultEventSelectionForSimulations("ALL");
+TString AliAnalysisMuMu::fgDefaultDimuonTriggerForSimulations("CMULLO-B-NOPF-MUON");
+
 Bool_t AliAnalysisMuMu::fgIsCompactGraphs(kFALSE);
 
 //_____________________________________________________________________________
+TString First(const TString s)
+{
+  TString rv;
+  
+  TObjArray* tokens = s.Tokenize(",");
+  
+  if (!tokens) return rv;
+  
+  rv = static_cast<TObjString*>(tokens->First())->String();
+  
+  delete tokens;
+  
+  return rv;
+}
+
+//_____________________________________________________________________________
 TString FindTrigger(const AliMergeableCollection& mc,
                     const char* base,
                     const char* selection,
@@ -125,7 +153,7 @@ TString FindTrigger(const AliMergeableCollection& mc,
 }
 
 //_____________________________________________________________________________
-AliAnalysisMuMu::AliAnalysisMuMu(const char* filename) : TObject(),
+AliAnalysisMuMu::AliAnalysisMuMu(const char* filename, const char* associatedSimFileName) : TObject(),
 fFilename(filename),
 fCounterCollection(0x0),
 fDimuonTriggers(fgDefaultDimuonTriggers),
@@ -134,14 +162,28 @@ fMinbiasTriggers(fgDefaultMinbiasTriggers),
 fEventSelectionList(fgDefaultEventSelectionList),
 fPairSelectionList(fgDefaultPairSelectionList),
 fCentralitySelectionList(fgDefaultCentralitySelectionList),
+fFitTypeList(fgDefaultFitTypeList),
 fBinning(0x0),
 fMergeableCollection(0x0),
 fRunNumbers(),
-fCorrectionPerRun(0x0)
+fCorrectionPerRun(0x0),
+fAssociatedSimulation(0x0)
 {
   // ctor
   
   GetCollections(fFilename,fMergeableCollection,fCounterCollection,fBinning,fRunNumbers);
+  
+  if (IsSimulation())
+  {
+    SetEventSelectionList("ALL");
+    SetDimuonTriggerList("CMULLO-B-NOPF-MUON");
+    SetFitTypeList("PSI1:1,COUNTJPSI:1");
+  }
+  
+  if ( strlen(associatedSimFileName) )
+  {
+    fAssociatedSimulation = new AliAnalysisMuMu(associatedSimFileName);
+  }
 }
 
 //_____________________________________________________________________________
@@ -152,6 +194,7 @@ AliAnalysisMuMu::~AliAnalysisMuMu()
   delete fBinning;
   delete fMergeableCollection;
   delete fCorrectionPerRun;
+  delete fAssociatedSimulation;
 }
 
 //_____________________________________________________________________________
@@ -360,6 +403,14 @@ void AliAnalysisMuMu::CentralityCheck(const char* filelist)
   
 }
 
+//_____________________________________________________________________________
+void AliAnalysisMuMu::CleanAllSpectra()
+{
+  /// Delete all the spectra we may have
+
+  MC()->RemoveByType("AliAnalysisMuMuSpectra");
+  Update();
+}
 
 //_____________________________________________________________________________
 void AliAnalysisMuMu::Compact(TGraph& g)
@@ -983,36 +1034,123 @@ void AliAnalysisMuMu::DrawMinv(const char* type,
                                const char* trigger,
                                const char* eventType,
                                const char* pairCut,
-                               const char* centrality) const
+                               const char* centrality,
+                               const char* subresultname,
+                               const char* flavour) const
 {
   /// Draw minv spectra for binning of given type
   
   if (!MC() || !BIN()) return;
   
-  TObjArray* bins = BIN()->CreateBinObjArray(particle,type);
+  TObjArray* bins = BIN()->CreateBinObjArray(particle,type,flavour);
   if (!bins)
   {
     AliError(Form("Could not get %s bins",type));
     return;
   }
   
+  Double_t xmin(-1);
+  Double_t xmax(-1);
+  
+  TString sparticle(particle);
+  if ( sparticle=="PSI" )
+  {
+    xmin = 2;
+    xmax = 6;
+  }
+  
+  Int_t nx(1);
+  Int_t ny(1);
+  
+  Int_t n = bins->GetEntries();
+  
+  if ( n == 2 )
+  {
+    nx = 2;
+  }
+  else if ( n > 2 )
+  {
+    ny = TMath::Nint(TMath::Sqrt(n));
+    nx = n/ny;
+  }
+  
+  TString stype(type);
+  stype.ToUpper();
+  
+  TString spectraName(Form("/%s/%s/%s/%s/%s-%s-%s",eventType,trigger,centrality,pairCut,particle,stype.Data(),flavour));
+  
+  AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(spectraName.Data()));
+  
+  AliDebug(1,Form("spectraName=%s spectra=%p",spectraName.Data(),spectra));
+
+  TObjArray* spectraBins(0x0);
+  if ( spectra )
+  {
+    spectraBins = spectra->Bins();
+  }
+  
+  TCanvas* c = new TCanvas;
+  c->Divide(nx,ny);
+  c->Draw();
+  gStyle->SetOptFit(1112);
+  
+  c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "AliAnalysisMuMu",
+              (void*)this, "ExecuteCanvasEvent(Int_t,Int_t,Int_t,TObject*)");
+
+  
   TIter next(bins);
   AliAnalysisMuMuBinning::Range* r;
+  Int_t ci(0);
   
   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
   {
     TString name(Form("/%s/%s/%s/%s/MinvUS%s",eventType,trigger,centrality,pairCut,r->AsString().Data()));
 
     AliDebug(1,name.Data());
-                 
+    
+    AliAnalysisMuMuResult* spectraBin(0x0);
+    
+    if ( spectraBins )
+    {
+      AliAnalysisMuMuResult* sr = static_cast<AliAnalysisMuMuResult*>(spectraBins->At(ci));
+      
+      spectraBin = sr->SubResult(subresultname);
+      
+      AliDebug(1,Form("spectraBin(%s)=%p",subresultname,spectraBin));
+    }
+    
     TH1* h = MC()->Histo(name.Data());
+    
+    if ( spectraBin )
+    {
+      h = spectraBin->Minv();
+    }
+    
     if (h)
     {
-      TString cname(name);
-      cname.ReplaceAll("/","_");
-      TCanvas* c = new TCanvas(cname.Data(),cname.Data());
-      c->SetLogy();
-      h->Draw("histe");
+      ++ci;
+      c->cd(ci);
+      gPad->SetLogy();
+      if (xmin>0)
+      {
+        h->GetXaxis()->SetRangeUser(xmin,xmax);
+      }
+      h->Draw("histes");
+      
+      TObject* f1 = h->GetListOfFunctions()->FindObject("fitTotal");
+      if (f1)
+      {
+        f1->Draw("same");
+      }
+      
+      gPad->Modified();
+      gPad->Update();
+      
+      TObject* stats = h->FindObject("stats");
+      if (stats)
+      {
+        stats->Draw("same");
+      }
     }
   }
   
@@ -1020,7 +1158,7 @@ void AliAnalysisMuMu::DrawMinv(const char* type,
 }
 
 //_____________________________________________________________________________
-void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle) const
+void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle, const char* flavour, const char* subresultname) const
 {
   /// Draw minv spectra for binning of given type
 
@@ -1028,7 +1166,63 @@ void AliAnalysisMuMu::DrawMinv(const char* type, const char* particle) const
            First(DimuonTriggerList()).Data(),
            First(EventSelectionList()).Data(),
            First(PairSelectionList()).Data(),
-           First(CentralitySelectionList()).Data());
+           First(CentralitySelectionList()).Data(),
+           subresultname,
+           flavour);
+}
+
+//___________________________________________________________________
+void AliAnalysisMuMu::ExecuteCanvasEvent(Int_t event, Int_t /*px*/, Int_t /*py*/, TObject *sel)
+{
+  // Actions in reponse to mouse button events.
+  
+  TCanvas* c = static_cast<TCanvas*>(gTQSender);
+  TPad* pad = static_cast<TPad*>(c->GetSelectedPad());
+  if (!pad) return;
+  
+//  if ((event == kButton1Down) ||
+  if (event == kButton1Double) 
+  {
+    
+//    Float_t x = pad->AbsPixeltoX(px);
+//    Float_t y = pad->AbsPixeltoY(py);
+//    x = pad->PadtoX(x);
+//    y = pad->PadtoY(y);
+
+//    std::cout << "event=" << event << " px=" << px << " py=" << py << " ";
+    
+    if ( sel && sel->InheritsFrom("TH1") )
+    {
+      TCanvas* clocal = new TCanvas;
+      clocal->SetLogy();
+      clocal->Draw();
+      sel->Draw();
+    }
+    else
+    {
+      TList* list = pad->GetListOfPrimitives();
+      TIter next(list);
+      TObject* h;
+      
+      while ( ( h = next() ) )
+      {
+        if ( h->InheritsFrom("TH1") )
+        {
+          TCanvas* clocal = new TCanvas;
+          clocal->SetLogy();
+          clocal->Draw();
+          h->Draw();
+          break;
+        }
+      }
+      
+    }
+
+//      std::cout  << std::endl;
+
+      pad->Modified();
+  }
+  
 }
 
 //_____________________________________________________________________________
@@ -1123,32 +1317,22 @@ AliAnalysisMuMu::FitParticle(const char* particle,
     return 0x0;
   }
   
-//  AliDebug(1,Form("will project fMergeableCollection=%p...",fMergeableCollection));
-//  
-//  AliMergeableCollection* hp = fMergeableCollection->Project(Form("/%s/%s/%s/%s/",eventType,trigger,centrality,pairCut));
-//  
-//  AliDebug(1,Form("projection=hp=%p",hp));
-  
-//  if (!hp)
-//  {
-//    AliError(Form("projection %s/%s/%s/%s failed",eventType,trigger,centrality,pairCut));
-//    delete bins;
-//    return 0x0;
-//  }
-  
-  binning.Print();
+//  binning.Print();
   
   AliAnalysisMuMuSpectra* spectra(0x0);
   
   AliAnalysisMuMuBinning::Range* bin;
   TIter next(bins);
-  Int_t nrebin=1;
+  
+  TObjArray* fitTypeArray = fFitTypeList.Tokenize(",");
+  TIter nextFitType(fitTypeArray);
+  TObjString* fitType;
+  TString flavour;
   
   while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next())) )
   {
     TString hname(Form("MinvUS%s",bin->AsString().Data()));
     
-//    TH1* hminv = hp->Histo(hname.Data());
     TH1* hminv = fMergeableCollection->Histo(Form("/%s/%s/%s/%s",eventType,trigger,centrality,pairCut),hname.Data());
     
     if (!hminv)
@@ -1178,24 +1362,38 @@ AliAnalysisMuMu::FitParticle(const char* particle,
     
     r->SetNofTriggers(ntrigger);
     
-//    r->AddFit("PSILOW",nrebin);
-    if (IsSimulation())
-    {
-      r->AddFit("PSI1",1);
-      r->AddFit("COUNTPSI",1);
-    }
-    else
+    nextFitType.Reset();
+    
+    while ( ( fitType = static_cast<TObjString*>(nextFitType())) )
     {
-      r->AddFit("PSILOW",nrebin+1);
+      AliDebug(1,Form("<<<<<< fitType=%s bin=%s",fitType->String().Data(),bin->Flavour().Data()));
+      
+      if ( fitType->String().BeginsWith("PSILOWMCTAILS") )
+      {
+        std::vector<Double_t> par;
+        par = GetMCCB2Tails(*bin);
+        if (!par.empty())
+        {
+          r->AddFit(fitType->String().Data(),par.size(),&par[0]);
+        }
+      }
+      else
+      {
+        r->AddFit(fitType->String());
+      }
     }
-//    r->AddFit("PSILOW",nrebin+2);
-//    r->AddFit("PSILOW",nrebin+3);
-    
-//    r->Print();
   
+    flavour = bin->Flavour();
+    
     if (!spectra)
     {
-      spectra = new AliAnalysisMuMuSpectra(binning.GetName());
+      TString spectraName(binning.GetName());
+      if ( flavour.Length() > 0 )
+      {
+        spectraName += "-";
+        spectraName += flavour;
+      }
+      spectra = new AliAnalysisMuMuSpectra(spectraName.Data());
     }
     
     spectra->AdoptResult(*bin,r);
@@ -1208,12 +1406,62 @@ AliAnalysisMuMu::FitParticle(const char* particle,
     
   } // loop on bins
   
+  delete fitTypeArray;
   delete bins;
   
   return spectra;
 }
 
 //_____________________________________________________________________________
+std::vector<Double_t>
+AliAnalysisMuMu::GetMCCB2Tails(const AliAnalysisMuMuBinning::Range& bin) const
+{
+  /// Get the tails from the associated simulation
+  
+  std::vector<Double_t> par;
+  
+  if (!SIM())
+  {
+    AliError("Cannot get MC tails without an associated simulation file !");
+    return par;
+  }
+
+
+  AliAnalysisMuMuSpectra* s = static_cast<AliAnalysisMuMuSpectra*>(SIM()->GetSpectra(bin.Type().Data(),bin.Flavour().Data()));
+  
+  if (!s)
+  {
+    AliError(Form("Could not find spectra %s,%s for associated simulation",bin.Type().Data(),bin.Flavour().Data()));
+    fAssociatedSimulation->MC()->Print("*:Ali*");
+    return par;
+  }
+  else
+  {
+    AliDebug(1,Form("AliAnalysisMuMuSpectra* s = reinterpret_cast<AliAnalysisMuMuSpectra*>(%p)",s));
+  }
+  
+  AliAnalysisMuMuResult* r = s->GetResultForBin(bin);
+  
+  if ( r )
+  {
+    AliAnalysisMuMuResult* r1 = r->SubResult("JPSI:1");
+    if  (r1)
+    {
+      TF1* func = static_cast<TF1*>(r1->Minv()->GetListOfFunctions()->FindObject("fitTotal"));
+      if (func)
+      {
+        par.push_back(func->GetParameter("alphaLow"));
+        par.push_back(func->GetParameter("nLow"));
+        par.push_back(func->GetParameter("alphaUp"));
+        par.push_back(func->GetParameter("nUp"));
+      }
+    }
+  }
+  
+  return par;
+}
+
+//_____________________________________________________________________________
 ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t runNumber)
 {
   // Get the expected (from OCDB scalers) trigger count
@@ -1240,6 +1488,36 @@ ULong64_t AliAnalysisMuMu::GetTriggerScalerCount(const char* triggerList, Int_t
 }
 
 //_____________________________________________________________________________
+AliAnalysisMuMuSpectra* AliAnalysisMuMu::GetSpectra(const char* what, const char* flavour) const
+{
+  /// get a given spectra
+  
+  TString swhat(what);
+  TString sflavour(flavour);
+  swhat.ToUpper();
+  sflavour.ToUpper();
+  
+  TString spectraName(Form("/PSALL/%s/PP/%s/PSI-%s",
+                           First(fgDefaultDimuonTriggers).Data(),
+                           First(fgDefaultPairSelectionList).Data(),
+                           swhat.Data()));
+
+  if (sflavour.Length()>0)
+  {
+    spectraName += "-";
+    spectraName += sflavour.Data();
+  }
+
+  if (IsSimulation())
+  {
+    spectraName.ReplaceAll("PSALL",fgDefaultEventSelectionForSimulations.Data());
+    spectraName.ReplaceAll(First(fgDefaultDimuonTriggers).Data(),fgDefaultDimuonTriggerForSimulations.Data());
+  }
+
+  return SPECTRA(spectraName.Data());
+}
+
+//_____________________________________________________________________________
 UInt_t AliAnalysisMuMu::GetSum(AliCounterCollection& cc, const char* triggerList,
                                const char* eventSelection, Int_t runNumber)
 {
@@ -1345,19 +1623,21 @@ Bool_t AliAnalysisMuMu::IsSimulation() const
 {
   // whether or not we have MC information
   
-  return ( fMergeableCollection->Histo("/INPUT/ALL/MinvUS") != 0x0 );
+  return ( fMergeableCollection->Histo(Form("/INPUT/%s/MinvUS",fgDefaultEventSelectionForSimulations.Data())) != 0x0 );
 }
 
 //_____________________________________________________________________________
 Int_t
-AliAnalysisMuMu::Jpsi(const char* what)
+AliAnalysisMuMu::Jpsi(const char* what, const char* binningFlavour)
 {
   // Fit the J/psi (and psiprime) peaks for the triggers in fDimuonTriggers list
-  // what="" => fit only fully integrated MinvUS
+  // what="integrated" => fit only fully integrated MinvUS
   // what="pt" => fit MinvUS in pt bins
   // what="y" => fit MinvUS in y bins
   // what="pt,y" => fit MinvUS in (pt,y) bins
   
+  TStopwatch timer;
+  
   if (!fMergeableCollection)
   {
     AliError("No mergeable collection. Consider Upgrade()");
@@ -1371,11 +1651,6 @@ AliAnalysisMuMu::Jpsi(const char* what)
   TObjArray* pairCutArray = fPairSelectionList.Tokenize(",");
   TObjArray* whatArray = TString(what).Tokenize(",");
   
-  if ( whatArray->GetEntries() <= 0 )
-  {
-    whatArray->Add(new TObjString(""));
-  }
-  
   TIter nextTrigger(triggerArray);
   TIter nextEventType(eventTypeArray);
   TIter nextPairCut(pairCutArray);
@@ -1392,7 +1667,7 @@ AliAnalysisMuMu::Jpsi(const char* what)
     
     if ( fBinning && swhat->String().Length() > 0 )
     {
-      binning = fBinning->Project("psi",swhat->String().Data());
+      binning = fBinning->Project("psi",swhat->String().Data(),binningFlavour);
     }
     else
     {
@@ -1469,13 +1744,18 @@ AliAnalysisMuMu::Jpsi(const char* what)
   delete triggerArray;
   delete eventTypeArray;
   delete pairCutArray;
-  
+
+  timer.Print();
+
   if (nfits)
   {
-    ReOpen(fFilename,"UPDATE");
-    fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
-    ReOpen(fFilename,"READ");
+    Update();
+//    ReOpen(fFilename,"UPDATE");
+//    fMergeableCollection->Write("MC",TObjArray::kOverwrite);// | TObjArray::kSingleKey);
+//    ReOpen(fFilename,"READ");
   }
+  
+  
   return nfits;
   
 }
@@ -1853,14 +2133,8 @@ TGraph* AliAnalysisMuMu::PlotEventSelectionEvolution(const char* trigger1, const
       {
         Double_t xcorr,ycorr;
         fCorrectionPerRun->GetPoint(i,xcorr,ycorr); // note that the fact that xcorr==runNumber has been checked by the SetCorrectionPerRun method
-        if ( TMath::Abs(ycorr)  > 1E-12 )
-        {
-          y /= ycorr;
-        }
-        else
-        {
-          y = 0.0;
-        }
+        y *= ycorr;
+        // FIXME: should get the correction error here
       }
       
       if ( asRejection ) y = 100*(1.0 - y);
@@ -2078,11 +2352,10 @@ TFile* AliAnalysisMuMu::ReOpen(const char* filename, const char* mode) const
   
   if (f)
   {
-    f->Close();
     delete f;
   }
   
-  f = TFile::Open(filename,"update");
+  f = TFile::Open(filename,mode);
   
   if ( !f || !f->IsOpen() )
   {
@@ -2126,7 +2399,7 @@ void AliAnalysisMuMu::SetColorScheme()
 }
 
 //_____________________________________________________________________________
-Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr)
+Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr, const char* formula)
 {
     /// Sets the graph used to correct values per run
   delete fCorrectionPerRun;
@@ -2148,8 +2421,30 @@ Bool_t AliAnalysisMuMu::SetCorrectionPerRun(const TGraph& corr)
     ++i;
   }
   
-  fCorrectionPerRun = static_cast<TGraph*>(corr.Clone());
+  fCorrectionPerRun = new TGraphErrors(corr.GetN());
+
+  TFormula* tformula(0x0);
+  if ( strlen(formula) > 0 )
+  {
+    tformula = new TFormula("SetCorrectionPerRunFormula",formula);
+  }
+
+  i = 0;
+  
+  for ( std::set<int>::const_iterator it = RunNumbers().begin(); it != RunNumbers().end(); ++it )
+  {
+    Double_t y = corr.GetY()[i];
+    
+    if ( tformula )
+    {
+      y = tformula->Eval(y);
+    }
+    fCorrectionPerRun->SetPoint(i,corr.GetX()[i],y);
+    ++i;
+  }
 
+  delete formula;
+  
   return kTRUE;
 }
 
@@ -2165,7 +2460,7 @@ void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuResult& r)
   if (!hinput)
   {
     AliError(Form("Got a simulation file where I did not find histogram /INPUT/INYRANGE/%s",hname.Data()));
-  
+
   }
   else
   {
@@ -2174,6 +2469,15 @@ void AliAnalysisMuMu::SetNofInputParticles(AliAnalysisMuMuResult& r)
 }
 
 //_____________________________________________________________________________
+AliAnalysisMuMuSpectra* AliAnalysisMuMu::SPECTRA(const char* fullpath) const
+{
+  /// Shortcut method to get to a spectra
+  if (!MC()) return 0x0;
+  
+  return static_cast<AliAnalysisMuMuSpectra*>(MC()->GetObject(fullpath));
+}
+
+//_____________________________________________________________________________
 void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
                                            Bool_t compact,
                                            Bool_t orderByTriggerCount)
@@ -2235,11 +2539,13 @@ void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
       total += n;
       totalExpected += expected;
       
-      msg += TString::Format("%30s %9lld expected %9lld ",strigger->String().Data(),n,expected);
+      msg += TString::Format("%30s %9lld expected %9lld [%s] ",strigger->String().Data(),n,expected,
+                             (n>expected ? "!" : " "));
       
       if ( expected > 0 ) {
         msg += TString::Format("fraction %5.1f %%",n*100.0/expected);
       }
+
       if (!compact)
       {
         msg += "\n";
@@ -2267,7 +2573,8 @@ void AliAnalysisMuMu::TriggerCountCoverage(const char* triggerList,
   {
     ++n;
     current += it->first;
-    std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",current*100.0/total,n ) << std::endl;
+    Double_t percent = ( total > 0.0 ? current*100.0/total : 0.0);
+    std::cout << Form("%10lld",it->first) << " " << it->second << " percentage of total = " << Form("%7.2f %% %3d",percent,n ) << std::endl;
   }
 
   std::cout << Form("--- TOTAL %lld expected %lld fraction %5.1f %%",
@@ -2288,6 +2595,18 @@ void AliAnalysisMuMu::UnsetCorrectionPerRun()
 }
 
 //_____________________________________________________________________________
+void AliAnalysisMuMu::Update()
+{
+  /// update the current file with memory
+  ReOpen(fFilename,"UPDATE");
+
+  MC()->Write("MC",TObject::kSingleKey);
+
+  ReOpen(fFilename,"READ");
+}
+
+//_____________________________________________________________________________
 Bool_t AliAnalysisMuMu::Upgrade(const char* filename)
 {
   /// Upgrade a file
@@ -2352,43 +2671,194 @@ Bool_t AliAnalysisMuMu::Upgrade()
 }
 
 //_____________________________________________________________________________
-AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* realFile, const char* simFile, const char* particle, const char* type)
+AliAnalysisMuMuSpectra* AliAnalysisMuMu::CorrectSpectra(const char* type, const char* flavour)
 {
+  /// Correct one spectra
   
-  AliAnalysisMuMu real(realFile);
-    AliAnalysisMuMu sim(simFile);
+  if (!SIM())
+  {
+    AliError("Cannot compute corrected yield without associated MC file !");
+    return 0x0;
+  }
+
+  const char* accEffSubResultName="COUNTJPSI:1";
   
+  AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
+  AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
   
-  AliAnalysisMuMuSpectra* realpt = static_cast<AliAnalysisMuMuSpectra*>(real.MC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/%s",type)));
-  AliAnalysisMuMuSpectra* simpt = static_cast<AliAnalysisMuMuSpectra*>(sim.MC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/%s",type)));
+  if ( !realSpectra )
+  {
+    AliError("could not get real spectra");
+    return 0x0;
+  }
   
-  if ( !realpt )
+  if ( !simSpectra)
   {
-    AliErrorClass("could not get real spectra");
+    AliError("could not get sim spectra");
     return 0x0;
   }
   
-  if ( !simpt )
+  realSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
+
+  Update();
+  
+  return realSpectra;
+}
+
+//_____________________________________________________________________________
+AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* type, const char* flavour)
+{
+  if (!SIM())
+  {
+    AliError("Cannot compute corrected yield without associated MC file !");
+    return 0x0;
+  }
+  
+  const char* accEffSubResultName="COUNTJPSI:1";
+  
+  AliAnalysisMuMuSpectra* realSpectra = GetSpectra(type,flavour);
+  
+  if ( !realSpectra )
+  {
+    AliError("could not get real spectra");
+    return 0x0;
+  }
+  
+  if (!realSpectra->HasValue("CoffNofJpsi"))
+  {
+    if (!CorrectSpectra(type,flavour))
+    {
+      AliError("Could not get corrected spectra");
+      return 0x0;
+    }
+  }
+  
+  AliAnalysisMuMuSpectra* simSpectra = SIM()->GetSpectra(type,flavour);
+  
+  if ( !simSpectra)
   {
     AliErrorClass("could not get sim spectra");
     return 0x0;
   }
   
-  AliInfoClass("REAL >>>");
-  realpt->Print("4");
-  AliInfoClass("SIM >>>");
-  simpt->Print("1");
+  Double_t nofCMUL7 = CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
+  Double_t nofCINT7 = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
+  Double_t nofCINT7w0MUL = CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
   
-  AliInfoClass("CORRECTED >>>");
-
-  AliAnalysisMuMuSpectra* spectra = static_cast<AliAnalysisMuMuSpectra*>(realpt->Clone());
-  spectra->Correct(*simpt,particle);
+  AliAnalysisMuMuBinning* binning = realSpectra->Binning();
+  TObjArray* bins = binning->CreateBinObjArray();
+  TIter nextBin(bins);
+  AliAnalysisMuMuBinning::Range* bin;
+  Int_t i(0);
+  AliAnalysisMuMuResult* r;
+  
+  while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
+  {
+    r = static_cast<AliAnalysisMuMuResult*>(realSpectra->Bins()->At(i));
+   
+    StdoutToAliDebug(1,std::cout << "bin=";r->Print(););
+    
+    AliAnalysisMuMuResult* rsim = static_cast<AliAnalysisMuMuResult*>(simSpectra->Bins()->At(i));
+    
+    Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
+    Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
+                                                                nofCINT7,TMath::Sqrt(nofCINT7),
+                                                                nofCMUL7,TMath::Sqrt(nofCMUL7));
+    
+    r->Set("MBR",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
+                                                                                                nofCINT7,TMath::Sqrt(nofCINT7)));
+    
+    Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
+    
+    Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
+                                                                 mbeq,mbeqError);
+    
+    r->Set("YJpsi",yield,yieldError);
+    
+    r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
+    r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
+    
+    ++i;
+  }
   
-  spectra->Print("1");
+  delete bins;
   
-  return spectra;
+  Update();
+  
+  return realSpectra;
 }
 
+////_____________________________________________________________________________
+//AliAnalysisMuMuSpectra* AliAnalysisMuMu::ComputeYield(const char* realFile, const char* simFile,
+//                                                      const  char* type)
+//{
+//  const char* accEffSubResultName="COUNTJPSI-1";
+//
+//  AliAnalysisMuMu real(realFile);
+//  AliAnalysisMuMu sim(simFile);
+//  
+//  
+//  AliAnalysisMuMuSpectra* realSpectra = static_cast<AliAnalysisMuMuSpectra*>(real.MC()->GetObject(Form("/PSALL/CMUL7-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
+//  AliAnalysisMuMuSpectra* simSpectra = static_cast<AliAnalysisMuMuSpectra*>(sim.MC()->GetObject(Form("/ALL/CMULLO-B-NOPF-MUON/PP/pMATCHLOWRABSBOTH/PSI-%s",type)));
+//  
+//  if ( !realSpectra )
+//  {
+//    AliErrorClass("could not get real spectra");
+//    return 0x0;
+//  }
+//  
+//  if ( !simSpectra)
+//  {
+//    AliErrorClass("could not get sim spectra");
+//    return 0x0;
+//  }
+//  
+//  AliAnalysisMuMuSpectra* corrSpectra = static_cast<AliAnalysisMuMuSpectra*>(realSpectra->Clone());
+//  corrSpectra->Correct(*simSpectra,"Jpsi",accEffSubResultName);
+//  
+//  Double_t nofCMUL7 = real.CC()->GetSum(Form("trigger:CMUL7-B-NOPF-MUON/event:PSALL"));
+//  Double_t nofCINT7 = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD/event:PSALL"));
+//  Double_t nofCINT7w0MUL = real.CC()->GetSum(Form("trigger:CINT7-B-NOPF-ALLNOTRD&0MUL/event:PSALL"));
+//  
+//  AliAnalysisMuMuBinning* binning = corrSpectra->Binning();
+//  TObjArray* bins = binning->CreateBinObjArray();
+//  TIter nextBin(bins);
+//  AliAnalysisMuMuBinning::Range* bin;
+//  Int_t i(0);
+//  AliAnalysisMuMuResult* r;
+//  
+//  while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
+//  {
+//    r = static_cast<AliAnalysisMuMuResult*>(corrSpectra->Bins()->At(i));
+//
+//    AliAnalysisMuMuResult* rsim = static_cast<AliAnalysisMuMuResult*>(simSpectra->Bins()->At(i));
+//    
+//    Double_t mbeq = nofCINT7w0MUL / ( nofCINT7 * nofCMUL7);
+//    Double_t mbeqError = mbeq * AliAnalysisMuMuResult::ErrorABC( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
+//                                                                nofCINT7,TMath::Sqrt(nofCINT7),
+//                                                                nofCMUL7,TMath::Sqrt(nofCMUL7));
+//    
+//    r->Set("MBR",nofCINT7/nofCINT7w0MUL,(nofCINT7/nofCINT7w0MUL)*AliAnalysisMuMuResult::ErrorAB( nofCINT7w0MUL, TMath::Sqrt(nofCINT7w0MUL),
+//                                                                                                nofCINT7,TMath::Sqrt(nofCINT7)));
+//    
+//    Double_t yield =  r->GetValue("CorrNofJpsi") * mbeq;
+//    
+//    Double_t yieldError = yield * AliAnalysisMuMuResult::ErrorAB( r->GetValue("CorrNofJpsi"), r->GetErrorStat("CorrNofJpsi"),
+//                                                                 mbeq,mbeqError);
+//    
+//    r->Set("YJpsi",yield,yieldError);
+//        
+//    r->Set("NofInputJpsi",rsim->GetValue("NofInputJpsi",accEffSubResultName),rsim->GetErrorStat("NofInputJpsi",accEffSubResultName));
+//    r->Set("AccEffJpsi",rsim->GetValue("AccEffJpsi",accEffSubResultName),rsim->GetErrorStat("AccEffJpsi",accEffSubResultName));
+//    
+//    ++i;
+//  }
+//
+//  delete bins;
+//  
+//  return corrSpectra;
+//}
+
 //_____________________________________________________________________________
 AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* realFile, const char* simFile, const char* type,
                                                const char* direction)
@@ -2399,7 +2869,7 @@ AliAnalysisMuMuSpectra* AliAnalysisMuMu::RABy(const char* realFile, const char*
   const Double_t ymax=TMath::Log(sqrts*1000.0/3.096916);
   const Double_t tab = 0.093e-6; // nb^-1
   const Double_t tabError = 0.0035E-6; // nb^-1
-  const char* accEffSubResultName="COUNTJPSI-1";
+  const char* accEffSubResultName="COUNTJPSI:1";
   
   TF1 ydist("ydist","[0]*TMath::Exp(-(x*x)/(2.0*0.39*0.39))",0.,0.5);
   ydist.SetParameter(0,1.);
index fbd0ddd..b0a92ca 100644 (file)
@@ -19,6 +19,7 @@
 #include <string>
 #include <TString.h>
 #include <vector>
+#include "RQ_OBJECT.h"
 
 class AliAnalysisMuMuResult;
 class AliAnalysisMuMuSpectra;
@@ -32,16 +33,22 @@ class TMap;
 
 class AliAnalysisMuMu : public TObject
 {
+
+  RQ_OBJECT("AliAnalysisMuMu")
   
+
 public:
   
-  enum EColor {
+  enum EColor
+  {
     kBlue=1500,
     kOrange=1501,
     kGreen=1502
   };
   
-  AliAnalysisMuMu(const char* filename="LHC12c_muons_AOD000_000179687.saf.root");
+  AliAnalysisMuMu(const char* filename="LHC12c_muons_AOD000_000179687.saf.root",
+                  const char* associatedSimFileName="");
+  
   virtual ~AliAnalysisMuMu();
   
   /* Basic checks */
@@ -79,9 +86,17 @@ public:
                                       const char* centrality,
                                       const AliAnalysisMuMuBinning& binning);
 
+  AliAnalysisMuMuSpectra* CorrectSpectra(const char* type, const char* flavour="");
+
+  AliAnalysisMuMuSpectra* ComputeYield(const char* type, const char* flavour="");
+  
+  void CleanAllSpectra();
+  
   ///------
   
-  static AliAnalysisMuMuSpectra* CorrectSpectra(const char* realFile="ds.list.saf.root", const char* simFile="ds.sim.list.saf.root", const char* particle="Jpsi", const char* type="PSI-PT");
+//  static AliAnalysisMuMuSpectra* ComputeYield(const char* realFile="ds.list.saf.root",
+//                                              const char* simFile="ds.sim.list.saf.root",
+//                                              const  char* type="PSI-Y VS PT");
 
   static AliAnalysisMuMuSpectra* RABy(const char* realFile="ds.list.saf.root", const char* simFile="ds.sim.list.saf.root", const char* type="", const char* direction="pPb");
 
@@ -119,11 +134,13 @@ public:
                                AliAnalysisMuMuBinning*& bin,
                                std::set<int>& runnumbers);
   
+  AliAnalysisMuMuSpectra* GetSpectra(const char* what, const char* flavour="") const;
+
   static UInt_t GetSum(AliCounterCollection& cc, const char* triggerList, const char* eventSelection, Int_t runNumber=-1);
   
   static ULong64_t GetTriggerScalerCount(const char* triggerList, Int_t runNumber);
   
-  Int_t Jpsi(const char* what="");
+  Int_t Jpsi(const char* what="integrated", const char* binningFlavour="");
   
   Bool_t IsSimulation() const;
   
@@ -149,14 +166,20 @@ public:
   TString CentralitySelectionList() const { return fCentralitySelectionList; }
   void SetCentralitySelectionList(const char* centralitySelectionList) { fCentralitySelectionList = centralitySelectionList; }
 
+  TString FitTypeList() const { return fFitTypeList; }
+  void SetFitTypeList(const char* fitTypelist) { fFitTypeList = fitTypelist; }
+  
   static void SetDefaultDimuonTriggerList(const char* dimuonTriggerList) { fgDefaultDimuonTriggers = dimuonTriggerList; }
   static void SetDefaultMuonTriggerList(const char* muonTriggerList) { fgDefaultMuonTriggers = muonTriggerList; }
   static void SetDefaultMinbiasTriggerList(const char* minbiasTriggerList) { fgDefaultMinbiasTriggers = minbiasTriggerList; }
   static void SetDefaultEventSelectionList(const char* eventSelectionList) { fgDefaultEventSelectionList = eventSelectionList; }
   static void SetDefaultPairSelectionList(const char* pairSelectionList) { fgDefaultPairSelectionList = pairSelectionList; }
   static void SetDefaultCentralitySelectionList(const char* centralitySelectionList) { fgDefaultCentralitySelectionList = centralitySelectionList; }
+  static void SetDefaultFitTypes(const char* fitTypes) { fgDefaultFitTypeList = fitTypes; }
   
-  
+  static void SetDefaultEventSelectionForSimulations(const char* value) { fgDefaultEventSelectionForSimulations = value; }
+  static void SetDefaultDimuonTriggerForSimulations(const char* value) { fgDefaultDimuonTriggerForSimulations = value; }
+
   AliMergeableCollection* MC() const { return fMergeableCollection; }
   AliCounterCollection* CC() const { return fCounterCollection; }
   AliAnalysisMuMuBinning* BIN() const { return fBinning; }
@@ -184,14 +207,26 @@ public:
                 const char* trigger,
                 const char* eventType,
                 const char* pairCut,
-                const char* centrality) const;
+                const char* centrality,
+                const char* subresultname="",
+                const char* flavour="") const;
 
-  void DrawMinv(const char* type="PT", const char* particle="PSI") const;
+  void DrawMinv(const char* type="PT", const char* particle="PSI", const char* flavour="", const char* subresultname="") const;
   
-  Bool_t SetCorrectionPerRun(const TGraph& corr);
+  Bool_t SetCorrectionPerRun(const TGraph& corr, const char* formula="");
   
   void UnsetCorrectionPerRun();
   
+  void ExecuteCanvasEvent(Int_t event, Int_t px, Int_t py, TObject *sel);
+
+  std::vector<Double_t> GetMCCB2Tails(const AliAnalysisMuMuBinning::Range& bin) const;
+  
+  AliAnalysisMuMu* SIM() const { return fAssociatedSimulation; }
+  
+  AliAnalysisMuMuSpectra* SPECTRA(const char* fullpath) const;
+  
+  void Update();
+
 private:
   AliAnalysisMuMu(const AliAnalysisMuMu& rhs); // not implemented on purpose
   AliAnalysisMuMu& operator=(const AliAnalysisMuMu& rhs); // not implemented on purpose
@@ -206,6 +241,7 @@ private:
 
   static TString ExpandPathName(const char* file);
   
+  
   TString fFilename; // file containing the result collections (of objects and counters) from AliAnalysisTaskMuMu
   AliCounterCollection* fCounterCollection; // collection of counters in file
   TString fDimuonTriggers; // list of dimuon triggers to consider
@@ -214,7 +250,8 @@ private:
   TString fEventSelectionList; // list of event types to consider
   TString fPairSelectionList; // list of pair cuts to consider
   TString fCentralitySelectionList; // list of centrality cuts to consider
-  
+  TString fFitTypeList; // list of fit types to perform
+
   AliAnalysisMuMuBinning* fBinning; // binning
   
   static TString fgOCDBPath; // OCDB to be used (raw:// by default)
@@ -225,6 +262,9 @@ private:
   static TString fgDefaultEventSelectionList; // default list of event selections
   static TString fgDefaultPairSelectionList; // default list of pair selections
   static TString fgDefaultCentralitySelectionList; // default list of centrality selections
+  static TString fgDefaultFitTypeList; // default list of fit types
+  static TString fgDefaultEventSelectionForSimulations; // default event selection (simulations)
+  static TString fgDefaultDimuonTriggerForSimulations; // default dimuon trigger (simulations)
   
   static Bool_t fgIsCompactGraphs; // whether the graph produced should be compact
   
@@ -234,7 +274,9 @@ private:
   
   TGraph* fCorrectionPerRun; // correction factor per run
   
-  ClassDef(AliAnalysisMuMu,10) // class to analysis results from AliAnalysisTaskMuMuXXX tasks
+  AliAnalysisMuMu* fAssociatedSimulation; // associated simulations (if any)
+
+  ClassDef(AliAnalysisMuMu,11) // class to analysis results from AliAnalysisTaskMuMuXXX tasks
 };
 
 #endif
index 78ccc1b..fd17027 100644 (file)
@@ -38,6 +38,7 @@ ClassImp(AliAnalysisMuMuResult)
 
 const std::map<std::string,Double_t>& MassMap()
 {
+  /// a simple map of masses...
   static std::map<std::string,Double_t> massMap;
   // not super elegant, but that way we do not depend on TDatabasePDG and thus
   // can decide on our particle naming
@@ -52,26 +53,11 @@ const std::map<std::string,Double_t>& MassMap()
 }
 
 
-Double_t ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr)
+Double_t funcCB(Double_t* xx, Double_t* par)
 {
-  Double_t e(0.0);
-  
-  if ( TMath::Abs(a) > 1E-12 )
-  {
-    e += (aerr*aerr)/(a*a);
-  }
+  /// Crystal ball
   
-  if ( TMath::Abs(b) > 1E-12 )
-  {
-    e += (berr*berr)/(b*b);
-  }
-  
-  return TMath::Sqrt(e);
-}
-
-Double_t funcCB(Double_t* xx, Double_t* par)
-{ 
-  Double_t N = par[0];
+  Double_t norm = par[0];
   Double_t alpha = par[1];
   Double_t n = par[2];
   Double_t mean = par[3];
@@ -79,23 +65,24 @@ Double_t funcCB(Double_t* xx, Double_t* par)
   
   Double_t x = xx[0];
   
-  Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
-  Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
+  Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
+  Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
   
   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
   
   if ( y > alpha*-1.0 ) 
   {
-    return N*TMath::Exp(-0.5*y*y);
+    return norm*TMath::Exp(-0.5*y*y);
   }
   else 
   {
-    return N*A*TMath::Power(B-y,-n);
+    return norm*a*TMath::Power(b-y,-n);
   }
 }
 
 Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
 {
+  /// crystal ball + expo + gaussian
   Double_t x = xx[0];
   
   Double_t g = par[0]*TMath::Gaus(x,par[1],par[2]);
@@ -108,8 +95,10 @@ Double_t funcJpsiGCBE(Double_t* xx, Double_t* par)
 }
 
 Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
-{ 
-  Double_t N = par[0];
+{
+  // custom fit for jpsi + psi prime
+  
+  Double_t norm = par[0];
   Double_t alpha = par[1];
   Double_t n = par[2];
   Double_t mean = par[3];
@@ -119,10 +108,10 @@ Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
   
   Double_t x = xx[0];
   
-  Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
-  Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
-  Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
-  Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
+  Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
+  Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
+  Double_t c = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
+  Double_t d = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
   
   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
   
@@ -130,15 +119,15 @@ Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
   
   if ( y > alphaprime )
   {
-    cb = N*C*TMath::Power(D+y,-nprime);
+    cb = norm*c*TMath::Power(d+y,-nprime);
   }
   else if ( y > alpha*-1.0 ) 
   {
-    cb = N*TMath::Exp(-0.5*y*y);
+    cb = norm*TMath::Exp(-0.5*y*y);
   }
   else 
   {
-    cb = N*A*TMath::Power(B-y,-n);
+    cb = norm*a*TMath::Power(b-y,-n);
   }
   
   if ( x < mean )
@@ -155,8 +144,10 @@ Double_t funcJpsiPsiPrimeCustom(Double_t* xx, Double_t* par)
 
 
 Double_t funcCB2(Double_t* xx, Double_t* par)
-{ 
-  Double_t N = par[0];
+{
+  /// CB2 = extended crystal ball
+  
+  Double_t norm = par[0];
   Double_t alpha = par[1];
   Double_t n = par[2];
   Double_t mean = par[3];
@@ -166,29 +157,31 @@ Double_t funcCB2(Double_t* xx, Double_t* par)
   
   Double_t x = xx[0];
   
-  Double_t A = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
-  Double_t B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
-  Double_t C = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
-  Double_t D = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
+  Double_t a = TMath::Power(n/TMath::Abs(alpha),n)*TMath::Exp(-0.5*alpha*alpha);
+  Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
+  Double_t c = TMath::Power(nprime/TMath::Abs(alphaprime),nprime)*TMath::Exp(-0.5*alphaprime*alphaprime);
+  Double_t d = nprime/TMath::Abs(alphaprime) - TMath::Abs(alphaprime);
   
   Double_t y = ( TMath::Abs(sigma) > 1E-12 ? (x-mean)/sigma : 0 );
   
   if ( y > alphaprime )
   {
-    return N*C*TMath::Power(D+y,-nprime);
+    return norm*c*TMath::Power(d+y,-nprime);
   }
   else if ( y > alpha*-1.0 ) 
   {
-    return N*TMath::Exp(-0.5*y*y);
+    return norm*TMath::Exp(-0.5*y*y);
   }
   else 
   {
-    return N*A*TMath::Power(B-y,-n);
+    return norm*a*TMath::Power(b-y,-n);
   }
 }
 
 Double_t funcJpsiPsiPrime(Double_t* xx, Double_t* par)
 {
+  /// CB + CB2
+  
   Double_t jpsi = funcCB(xx,par);
   Double_t psiprime = funcCB2(xx,par+5);
   
@@ -221,6 +214,8 @@ Double_t funcJpsiCBE(Double_t* xx, Double_t* par)
 
 Double_t funcJpsiPCBE(Double_t* xx, Double_t* par)
 {
+  // CB + expo + pol2
+  
   Double_t x = xx[0];
 
   Double_t pol2 = par[0] + par[1]*x + par[2]*x*x;
@@ -247,8 +242,54 @@ Double_t funcJpsiECBE(Double_t* xx, Double_t* par)
   return e1+e2+jpsi;
 }
 
+Double_t funcJpsiNA48(Double_t*x, Double_t* par)
+{
+  /// Fit function from e.q. 4.8 of Ruben's PhD.
+  Double_t c1 = par[0];
+  Double_t c2 = par[1];
+  Double_t d1 = par[2];
+  Double_t d2 = par[3];
+  Double_t g1 = par[4];
+  Double_t g2 = par[5];
+  Double_t m0 = par[6];
+  Double_t sigma1 = par[7];
+  Double_t sigma2 = par[8];
+  Double_t b1 = par[9];
+  Double_t b2 = par[10];
+  Double_t norm = par[11];
+  
+  Double_t m = x[0];
+  
+  Double_t rv(0);
+  
+  if ( m <= c1*m0 )
+  {
+    Double_t e = d1-g1*TMath::Sqrt(c1*m0-m);
+    rv = TMath::Power(b1*(c1*m0-m),e);
+    rv += sigma1;
+  }
+  else if( m >= c1*m0 && m <= m0 )
+  {
+    rv = sigma1;
+  }
+  else if ( m >= m0 && m < c2*m0 )
+  {
+    rv = sigma2;
+  }
+  else if( m >= c2*m0 )
+  {
+    Double_t e = d2-g2*TMath::Sqrt(m-c2*m0);
+    rv = TMath::Power(b2*(m-c2*m0),e);
+    rv += sigma2;
+  }
+  
+  return norm*TMath::Exp(-(m-m0)*(m-m0)/(2.0*rv*rv));
+}
+
 const char* NormalizeName(const char* name, const char* suffix)
 {
+  /// Remove - and / from the name, and adds _suffix
+  
   TString str(Form("%s_%s",name,suffix));
   
   str.ReplaceAll("-","_");
@@ -346,6 +387,9 @@ Double_t fitFunctionCB2VWG(Double_t *x, Double_t *par)
 //------------------------------------------------------------------------------
 Double_t func2CB2VWG(Double_t *x, Double_t *par)
 {
+  /// 2 extended crystal balls + variable width gaussian
+  /// width of the second CB related to the first (free) one.
+  
   Double_t par2[7] = {
     par[11],
     3.68609+(par[5]-3.096916)/3.096916*3.68609,
@@ -375,7 +419,13 @@ fSubResults(0x0),
 fMap(0x0),
 fMother(0x0),
 fKeys(0x0),
-fWeight(0.0)
+fWeight(0.0),
+fRebin(0),
+fTriggerClass(),
+fEventSelection(),
+fPairSelection(),
+fCentralitySelection(),
+fAlias()
 {
 }
 
@@ -391,8 +441,13 @@ AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv)
   fMap(0x0),
   fMother(0x0),
 fKeys(0x0),
-fWeight(0.0)
-
+fWeight(0.0),
+fRebin(0),
+fTriggerClass(),
+fEventSelection(),
+fPairSelection(),
+fCentralitySelection(),
+fAlias()
 {
   SetMinv(hminv);
 }
@@ -402,7 +457,7 @@ AliAnalysisMuMuResult::AliAnalysisMuMuResult(const TH1& hminv,
                                              const char* fitType,
                                              Int_t nrebin)
 :
-TNamed(Form("%s-%d",fitType,nrebin),""),
+TNamed(Form("%s:%d",fitType,nrebin),""),
 fNofRuns(1),
 fNofTriggers(-1),
 fMinv(0x0),
@@ -411,8 +466,13 @@ fSubResults(0x0),
 fMap(0x0),
 fMother(0x0),
 fKeys(0x0),
-fWeight(0.0)
-
+fWeight(0.0),
+fRebin(nrebin),
+fTriggerClass(),
+fEventSelection(),
+fPairSelection(),
+fCentralitySelection(),
+fAlias()
 {
   SetMinv(hminv);
 }
@@ -434,8 +494,13 @@ fSubResults(0x0),
 fMap(0x0),
 fMother(0x0),
 fKeys(0x0),
-fWeight(0.0)
-
+fWeight(0.0),
+fRebin(1),
+fTriggerClass(triggerName),
+fEventSelection(eventSelection),
+fPairSelection(pairSelection),
+fCentralitySelection(centSelection),
+fAlias()
 {
   SetMinv(hminv);
 }
@@ -452,7 +517,9 @@ fSubResults(0x0),
 fMap(0x0),
 fMother(0x0),
 fKeys(0x0),
-fWeight(rhs.fWeight)
+fWeight(rhs.fWeight),
+fRebin(rhs.fRebin),
+fAlias()
 {
   /// copy ctor
   /// Note that the mother is lost
@@ -473,11 +540,17 @@ fWeight(rhs.fWeight)
     fMap = static_cast<TMap*>(rhs.fMap->Clone());
   }
   
+  if ( rhs.fAlias.Length() > 0 )
+  {
+    fAlias = rhs.fAlias;
+  }
 }
 
 //_____________________________________________________________________________
 AliAnalysisMuMuResult& AliAnalysisMuMuResult::operator=(const AliAnalysisMuMuResult& rhs)
 {
+  /// Assignment operator
+  
   if (this!=&rhs)
   {
     delete fMinv;
@@ -510,6 +583,14 @@ AliAnalysisMuMuResult& AliAnalysisMuMuResult::operator=(const AliAnalysisMuMuRes
     fNofTriggers = rhs.NofTriggers();
     fBin = rhs.Bin();
     fWeight = rhs.fWeight;
+    fRebin = rhs.fRebin;
+    
+    fAlias="";
+    
+    if ( rhs.fAlias.Length() > 0 )
+    {
+      fAlias = rhs.fAlias;
+    }
   }
   
   return *this;
@@ -518,6 +599,7 @@ AliAnalysisMuMuResult& AliAnalysisMuMuResult::operator=(const AliAnalysisMuMuRes
 //_____________________________________________________________________________
 AliAnalysisMuMuResult::~AliAnalysisMuMuResult()
 {
+  // dtor
   delete fMap;
   delete fMinv;
   delete fSubResults;
@@ -527,6 +609,7 @@ AliAnalysisMuMuResult::~AliAnalysisMuMuResult()
 //_____________________________________________________________________________
 const AliAnalysisMuMuBinning::Range& AliAnalysisMuMuResult::Bin() const
 {
+  /// Get the bin of this result
   if ( !Mother() ) return fBin;
   else return Mother()->Bin();
 }
@@ -534,6 +617,7 @@ const AliAnalysisMuMuBinning::Range& AliAnalysisMuMuResult::Bin() const
 //_____________________________________________________________________________
 TObject* AliAnalysisMuMuResult::Clone(const char* /*newname*/) const
 {
+  /// Clone this result
   return new AliAnalysisMuMuResult(*this);
 }
 
@@ -564,7 +648,10 @@ Bool_t AliAnalysisMuMuResult::Correct(const AliAnalysisMuMuResult& other, const
     
     return kTRUE;
   }
-  
+  else
+  {
+    AliError(Form("Result does not have Nof%s : cannot correct it !",particle));
+  }
   return kFALSE;
 }
 
@@ -880,6 +967,8 @@ AliAnalysisMuMuResult::SubResult AliAnalysisMuMuResult::FitJpsiPsiPrimeCB(TH1& h
 //_____________________________________________________________________________
 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiGCBE(TH1& h)
 {
+  /// Fit Jpsi spectra with crystal ball + gaussian + exponential
+  
   std::cout << "Fit with jpsi alone (gaus + CB + expo)" << std::endl;
   
   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
@@ -1221,6 +1310,8 @@ void AliAnalysisMuMuResult::FitJpsiECBE(TH1& h)
 //_____________________________________________________________________________
 AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi(TH1& h)
 {
+  /// Fit Jpsi spectra using extended crystall ball (CB2) with free tails
+  
   std::cout << "Fit with jpsi alone" << std::endl;
 
   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
@@ -1233,15 +1324,15 @@ AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi(TH1& h)
   const Double_t xmax(8.0);
 
   TF1* fitTotal = new TF1("fitTotal",funcCB2,xmin,xmax,7);
-  fitTotal->SetParNames("N","alpha","n","mean","sigma","alphaprime","nprime");
-  fitTotal->SetParameters(h.GetMaximum(),1,5,3.0,0.07,1.5,3);
+  fitTotal->SetParNames("N","alphaLow","nLow","mean","sigma","alphaUp","nUp");
+  fitTotal->SetParameters(h.GetMaximum(),1,5,3.1,0.07,1.5,3);
   fitTotal->SetParLimits(0,0,h.GetMaximum()*2); // N
   fitTotal->SetParLimits(1,0,10); // alpha
-  fitTotal->SetParLimits(2,1,10); // n
-  fitTotal->SetParLimits(3,1,4); // mean
+  fitTotal->SetParLimits(2,0.1,10); // n
+  fitTotal->SetParLimits(3,3,3.1); // mean
   fitTotal->SetParLimits(4,0.01,1); // sigma
   fitTotal->SetParLimits(5,0,10); // alpha
-  fitTotal->SetParLimits(6,1,10); // n
+  fitTotal->SetParLimits(6,0.1,10); // n
   
   hfit->Fit(fitTotal,"+","",2,5);
   
@@ -1259,17 +1350,157 @@ AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi(TH1& h)
 }
 
 //_____________________________________________________________________________
-AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h)
+AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiCB2VWG(const TH1& h)
 {
-  std::cout << "Fit with jpsi + psiprime VWG" << std::endl;
+  /// Fit Jpsi spectra using extended crystal ball (CB2) + variable width gaussian (VWG)
+  
+  std::cout << "Fit with jpsi VWG" << std::endl;
   
   Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
   
-  AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSI2CB2VWG",nrebin);
+  AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSICB2VWG",nrebin);
   
   
   TH1* hfit = r->Minv();
   
+  const Double_t xmin(2.0);
+  const Double_t xmax(5.0);
+  
+//  // gaussian variable width
+//  Double_t sigma = par[2]+par[3]*((x[0]-par[1])/par[1]);
+//  return par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*sigma*sigma));
+//  Double_t CrystalBallExtended(Double_t *x,Double_t *par)
+//  //par[0] = Normalization
+//  //par[1] = mean
+//  //par[2] = sigma
+//  //par[3] = alpha
+//  //par[4] = n
+//  //par[5] = alpha'
+//  //par[6] = n'
+
+  TF1* fitTotal = new TF1("fitTotal",fitFunctionCB2VWG,xmin,xmax,11);
+  fitTotal->SetParNames("kVWG","mVWG","sVWG1","sVWG2","norm","mean","sigma","alpha","n","alpha'","n'");
+  
+  fitTotal->SetParameter(0, 10000.); // kVWG
+  fitTotal->SetParameter(1, 1.9); // mVWG
+  
+  fitTotal->SetParameter(2, 0.5); // sVWG1
+  fitTotal->SetParLimits(2, 0., 100.);
+  
+  fitTotal->SetParameter(3, 0.3); // sVWG2
+  fitTotal->SetParLimits(3, 0., 100.);
+  
+  fitTotal->SetParameter(4, h.GetMaximum()); // norm
+  
+  fitTotal->SetParameter(5, 3.1); // mean
+  fitTotal->SetParLimits(5, 3.0, 3.2);
+  
+  fitTotal->SetParameter(6, 0.08); // sigma
+  fitTotal->SetParLimits(6, 0.04, 0.20);
+  
+  fitTotal->SetParameter(7,1.0); // alpha
+  fitTotal->SetParameter(8,5); // n
+  fitTotal->SetParameter(9,2.0); // alpha'
+  fitTotal->SetParameter(10,4); // n'
+  
+//  fitTotal->FixParameter(7, 0.93);
+//  fitTotal->FixParameter(8, 5.59);
+//  fitTotal->FixParameter(9, 2.32);
+//  fitTotal->FixParameter(10, 3.39);
+//  fitTotal->SetParameter(11, 10.);
+  
+  const char* fitOption = "SIER"; //+";
+  
+  TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
+  
+  r->Set("MeanJpsi",fitTotal->GetParameter(5),fitTotal->GetParError(5));
+  r->Set("SigmaJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
+  
+  double m = r->GetValue("MeanJpsi");
+  double s = r->GetValue("SigmaJpsi");
+  double n = 3.0;
+  
+  TF1* fcb = new TF1("fcb",CrystalBallExtended,xmin,xmax,7);
+  fcb->SetParameters(fitTotal->GetParameter(4),
+                     fitTotal->GetParameter(5),
+                     fitTotal->GetParameter(6),
+                     fitTotal->GetParameter(7),
+                     fitTotal->GetParameter(8),
+                     fitTotal->GetParameter(9),
+                     fitTotal->GetParameter(10));
+  
+  
+  fcb->SetLineColor(1);
+  fcb->SetNpx(1000);
+  TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(4));
+  TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(4));
+  l1->SetLineColor(6);
+  l2->SetLineColor(6);
+  hfit->GetListOfFunctions()->Add(fcb);
+  hfit->GetListOfFunctions()->Add(l1);
+  hfit->GetListOfFunctions()->Add(l2);
+  
+  Double_t cbParameters[7];
+  Double_t covarianceMatrix[7][7];
+  
+  for ( int ix = 0; ix < 7; ++ix )
+  {
+    cbParameters[ix] = fitTotal->GetParameter(ix+4);
+  }
+  
+  for ( int iy = 0; iy < 5; ++iy )
+  {
+    for ( int ix = 0; ix < 5; ++ix )
+    {
+      covarianceMatrix[ix][iy] = (fitResult->GetCovarianceMatrix())(ix+4,iy+4);
+    }
+  }
+  
+  double njpsi = fcb->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
+  double nerr = fcb->IntegralError(m-n*s,m+n*s,&cbParameters[0],&covarianceMatrix[0][0])/h.GetBinWidth(1);
+  
+  r->Set("NofJpsi",njpsi,nerr);
+  
+  return r;
+}
+
+//_____________________________________________________________________________
+AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h,
+                                                             Double_t alphaLow,
+                                                             Double_t nLow,
+                                                             Double_t alphaUp,
+                                                             Double_t nUp)
+{
+  /// Fit using extended crystal ball + variable width gaussian
+  
+  std::cout << Form("Fit with jpsi + psiprime VWG alphaLow=%5.2f nLow=%5.2f alphaUp=%5.2f nUp=%5.2f",
+                    alphaLow,nLow,alphaUp,nUp) << std::endl;
+  
+  Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
+  
+  TString resultName("JPSI2CB2VWG");
+  if ( alphaLow > 0 )
+  {
+    resultName += TString::Format("alphaLow=%5.2f",alphaLow);
+  }
+  if ( nLow > 0 )
+  {
+    resultName += TString::Format("nLow=%5.2f",nLow);
+  }
+  if ( alphaUp > 0 )
+  {
+    resultName += TString::Format("alphaUp=%5.2f",alphaUp);
+  }
+  if ( nUp > 0 )
+  {
+    resultName += TString::Format("nUp=%5.2f",nUp);
+  }
+  resultName.ReplaceAll(" ","");
+  
+  AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,resultName.Data(),nrebin);
+  
+  TH1* hfit = r->Minv();
+  
   const Double_t xmin(2.2);
   const Double_t xmax(5.0);
   
@@ -1284,17 +1515,56 @@ AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h)
   fitTotal->SetParameter(3, 0.3);
   fitTotal->SetParLimits(3, 0., 100.);
   fitTotal->SetParameter(4, 100.);
-  fitTotal->SetParameter(5, 3.13);
+  fitTotal->SetParameter(5, 3.1);
   fitTotal->SetParLimits(5, 3.08, 3.2);
   fitTotal->SetParameter(6, 0.08);
   fitTotal->SetParLimits(6, 0.05, 0.15);
-  fitTotal->FixParameter(7, 0.93);
-  fitTotal->FixParameter(8, 5.59);
-  fitTotal->FixParameter(9, 2.32);
-  fitTotal->FixParameter(10, 3.39);
+  
+//  r = FitJpsi2CB2VWG(*hminv,0.93,5.59,2.32,3.39);
+
+  if ( alphaLow > 0 )
+  {
+    fitTotal->FixParameter(7, alphaLow);
+  }
+  else
+  {
+    fitTotal->SetParameter(7,0.9);
+    fitTotal->SetParLimits(7,0.1,10.0);
+  }
+  
+  if ( nLow > 0 )
+  {
+    fitTotal->FixParameter(8, nLow);
+  }
+  else
+  {
+    fitTotal->SetParameter(8,5.0);
+    fitTotal->SetParLimits(8,0.0,10.0);
+  }
+  
+  if ( alphaUp > 0 )
+  {
+    fitTotal->FixParameter(9, alphaUp);
+  }
+  else
+  {
+    fitTotal->SetParameter(9, 2.0);
+    fitTotal->SetParLimits(9,0.1,10.0);
+  }
+  
+  if ( nUp > 0 )
+  {
+    fitTotal->FixParameter(10, nUp);    
+  }
+  else
+  {
+    fitTotal->SetParameter(10,3.0);
+    fitTotal->SetParLimits(10,0.0,10.0);
+  }
+  
   fitTotal->SetParameter(11, 10.);
 
-  const char* fitOption = "SIER"; //+";
+  const char* fitOption = "SER"; //+";
   
   TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
   
@@ -1349,6 +1619,97 @@ AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsi2CB2VWG(const TH1& h)
   return r;
 }
 
+//_____________________________________________________________________________
+AliAnalysisMuMuResult* AliAnalysisMuMuResult::FitJpsiNA48(const TH1& h)
+{
+  /// fit using functional form from Ruben Shahoyan's thesis (2001) (eq. 4.8.)
+  
+  std::cout << "Fit with jpsi NA50 Ruben eq. 4.8" << std::endl;
+  
+  Int_t nrebin = fMinv->GetXaxis()->GetNbins() / h.GetXaxis()->GetNbins();
+  
+  AliAnalysisMuMuResult* r = new AliAnalysisMuMuResult(h,"JPSINA",nrebin);
+  
+  TH1* hfit = r->Minv();
+  
+  const Double_t xmin(2.0);
+  const Double_t xmax(5.0);
+  
+  TF1* fitTotal = new TF1("fitTotal",funcJpsiNA48,xmin,xmax,12);
+  
+  fitTotal->SetParName( 0, "c1");
+  fitTotal->FixParameter(0,0.97);
+  
+  fitTotal->SetParName( 1, "c2");
+  fitTotal->FixParameter(1,1.05);
+  
+  fitTotal->SetParName( 2, "d1");
+  fitTotal->SetParameter(2,0.0);
+  fitTotal->SetParLimits(2,0,1);
+  
+  fitTotal->SetParName( 3, "d2");
+  fitTotal->SetParameter(3,0.0);
+  fitTotal->SetParLimits(3,0,1);
+  
+  fitTotal->SetParName( 4, "g1");
+  fitTotal->SetParameter(4,0.0);
+  fitTotal->SetParLimits(4,0.0,1);
+  
+  fitTotal->SetParName( 5, "g2");
+  fitTotal->SetParameter(5,0.0);
+  fitTotal->SetParLimits(5,0.0,1);
+  
+  fitTotal->SetParName( 6, "m0");
+  fitTotal->SetParameter(6,3.1);
+  fitTotal->SetParLimits(6,2.8,3.4);
+
+  fitTotal->SetParName( 7, "sigma1");
+  fitTotal->SetParameter(7,0.05);
+  fitTotal->SetParLimits(7,0.01,0.2);
+  
+  fitTotal->SetParName( 8, "sigma2");
+  fitTotal->SetParameter(8,0.05);
+  fitTotal->SetParLimits(8,0.01,0.2);
+
+  fitTotal->SetParName( 9, "b1");
+  fitTotal->SetParameter(9,1.0);
+  fitTotal->SetParLimits(9,0,1);
+  
+  fitTotal->SetParName(10, "b2");
+  fitTotal->SetParameter(10,1.0);
+  fitTotal->SetParLimits(10,0,1);
+  
+  fitTotal->SetParName(11, "norm");
+  fitTotal->SetParameter(11,h.GetMaximum());
+  
+  const char* fitOption = "SIER"; //+";
+  
+  TFitResultPtr fitResult = hfit->Fit(fitTotal,fitOption,"");
+  
+  r->Set("MeanJpsi",fitTotal->GetParameter(6),fitTotal->GetParError(6));
+  r->Set("SigmaJpsi",
+         fitTotal->GetParameter(7)+fitTotal->GetParameter(8),
+         0.0);
+
+  double m = r->GetValue("MeanJpsi");
+  double s = r->GetValue("SigmaJpsi");
+  double n = 3.0;
+  
+  TLine* l1 = new TLine(m-n*s,0,m-n*s,fitTotal->GetParameter(11));
+  TLine* l2 = new TLine(m+n*s,0,m+n*s,fitTotal->GetParameter(11));
+  l1->SetLineColor(6);
+  l2->SetLineColor(6);
+  hfit->GetListOfFunctions()->Add(l1);
+  hfit->GetListOfFunctions()->Add(l2);
+  
+  double njpsi = fitTotal->Integral(m-n*s,m+n*s)/h.GetBinWidth(1);
+  double nerr = fitTotal->IntegralError(m-n*s,m+n*s)/h.GetBinWidth(1);
+  
+  r->Set("NofJpsi",njpsi,nerr);
+  
+  return r;
+}
+
 /*
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::FitUpsilon(TH1& h)
@@ -1386,6 +1747,8 @@ void AliAnalysisMuMuResult::FitUpsilon(TH1& h)
 //_____________________________________________________________________________
 Double_t AliAnalysisMuMuResult::ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr)
 {
+  /// Compute the quadratic sum of 2 errors
+
   Double_t e(0.0);
   
   if ( TMath::Abs(a) > 1E-12 )
@@ -1404,6 +1767,8 @@ Double_t AliAnalysisMuMuResult::ErrorAB(Double_t a, Double_t aerr, Double_t b, D
 //_____________________________________________________________________________
 Double_t AliAnalysisMuMuResult::ErrorABC(Double_t a, Double_t aerr, Double_t b, Double_t berr, Double_t c, Double_t cerror)
 {
+  /// Compute the quadratic sum of 3 errors
+  
   Double_t e(0.0);
   
   if ( TMath::Abs(a) > 1E-12 )
@@ -1426,15 +1791,37 @@ Double_t AliAnalysisMuMuResult::ErrorABC(Double_t a, Double_t aerr, Double_t b,
 
 
 //_____________________________________________________________________________
-Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t nrebin)
+Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t npar, Double_t* par)
 {
-  //  new TCanvas;
+  // Add a fit to this result
+  
+  TString msg(Form("fitType=%s npar=%d par[]=",fitType,npar));
+  
+  for ( Int_t i = 0; i < npar; ++i )
+  {
+    msg += TString::Format("%e,",par[i]);
+  }
+  
+  msg += TString::Format(" minv=%p %d",fMinv,fMinv?TMath::Nint(fMinv->GetEntries()):0);
   
   if ( !fMinv ) return kFALSE;
   
   TString sFitType(fitType);
   sFitType.ToUpper();
+  Int_t nrebin(1);
+  
+  if (sFitType.CountChar(':'))
+  {
+    Int_t index = sFitType.Index(":");
+    nrebin = TString(sFitType(index+1,sFitType.Length()-index-1)).Atoi();
+    sFitType = sFitType(0,index);
+  }
+  
+  msg += TString::Format(" nrebin=%d",nrebin);
   
+  AliDebug(1,msg.Data());
+  
+
   if ( fMinv->GetEntries()<100 && !sFitType.Contains("COUNT")) return kFALSE;
   
   TH1* hminv = static_cast<TH1*>(fMinv->Clone());
@@ -1444,39 +1831,55 @@ Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t nrebin)
 
   AliAnalysisMuMuResult* r(0x0);
   
-//  if ( sFitType == "PSI12" )
-//  {
-//    r = FitJpsiPsiPrimeCustom(*hminv);
-//  }
-//
   if ( sFitType=="PSI1")
-  {//
+  {
     r = FitJpsi(*hminv);
-//    r = FitJpsiGCBE(*hminv);
   }
-  
-  if ( sFitType == "PSILOW")
+  else if ( sFitType == "PSILOW")
   {
-    r = FitJpsi2CB2VWG(*hminv);
+    r = FitJpsi2CB2VWG(*hminv,-1,-1,-1,-1); // free tails
   }
-
-  if ( sFitType == "COUNTPSI" )
+  else if ( sFitType == "PSILOWMCTAILS" )
+  {
+    if ( npar!= 4 )
+    {
+      AliError("Cannot use PSILOWMCTAILS without being given the MC tails !");
+      delete hminv;
+      return kFALSE;
+    }
+    r = FitJpsi2CB2VWG(*hminv,par[0],par[1],par[2],par[3]);
+    if (r)
+    {
+      r->SetAlias("MCTAILS");
+    }
+  }
+  else if ( sFitType.BeginsWith("PSILOWALPHA") )
+  {
+    Float_t lpar[] = { 0.0, 0.0, 0.0, 0.0 };
+    
+    AliDebug(1,Form("sFitType=%s",sFitType.Data()));
+    
+    sscanf(sFitType.Data(),"PSILOWALPHALOW%fNLOW%fALPHAUP%fNUP%f",
+           &lpar[0],&lpar[1],&lpar[2],&lpar[3]);
+    
+    AliDebug(1,Form("PSILOW ALPHALOW=%f NLOW=%f ALPHAUP=%f NUP=%f",lpar[0],lpar[1],lpar[2],lpar[3]));
+    
+    if ( lpar[0] == 0.0 && lpar[1] == 0.0 && lpar[0] == 0.0 && lpar[1] == 0.0 )
+    {
+      AliError("Cannot work with zero tails !");
+    }
+    else
+    {
+      r = FitJpsi2CB2VWG(*hminv,lpar[0],lpar[1],lpar[2],lpar[3]);      
+    }
+  }
+  else if ( sFitType == "COUNTJPSI" )
   {
     r = new AliAnalysisMuMuResult(*hminv,"COUNTJPSI",1);
     Double_t n = CountParticle(*hminv,"Jpsi");
     r->Set("NofJpsi",n,TMath::Sqrt(n));
   }
   
-//  if (r)
-//  {
-//    r->SetNofInputParticles("Jpsi",GetValue("NofInputJpsi"));
-//  }
-
-  if (!fSubResults)
-  {
-    fSubResults = new TObjArray;
-    fSubResults->SetOwner(kTRUE);
-  }
   
   if ( r )
   {
@@ -1484,6 +1887,13 @@ Bool_t AliAnalysisMuMuResult::AddFit(const char* fitType, Int_t nrebin)
     r->SetBin(Bin());
     r->SetNofTriggers(NofTriggers());
     r->SetNofRuns(NofRuns());
+
+    if (!fSubResults)
+    {
+      fSubResults = new TObjArray;
+      fSubResults->SetOwner(kTRUE);
+    }
+
     fSubResults->Add(r);
   }
   
@@ -1590,6 +2000,9 @@ Double_t AliAnalysisMuMuResult::GetValue(const char* name, const char* subResult
 //_____________________________________________________________________________
 Bool_t AliAnalysisMuMuResult::HasValue(const char* name, const char* subResultName) const
 {
+  /// Whether this result (or subresult if subResultName is provided) has a property
+  /// named "name"
+  
   if ( strlen(subResultName) > 0 )
   {
     if ( !fSubResults)
@@ -1732,6 +2145,11 @@ Long64_t AliAnalysisMuMuResult::Merge(TCollection* list)
     {
       AliAnalysisMuMuResult* result = dynamic_cast<AliAnalysisMuMuResult*>(currObj);
     
+      if (!result)
+      { 
+        continue;
+      }
+      
       if (!result->HasValue(key->String()))
       {
         AliError(Form("Did not find key %s in of the result to merge",key->String().Data()));
@@ -1814,6 +2232,7 @@ Long64_t AliAnalysisMuMuResult::Merge(TCollection* list)
 //_____________________________________________________________________________
 Int_t AliAnalysisMuMuResult::NofRuns() const
 {
+  /// Get the number of runs
   if ( !Mother() ) return fNofRuns;
   else return Mother()->NofRuns();
 }
@@ -1821,6 +2240,8 @@ Int_t AliAnalysisMuMuResult::NofRuns() const
 //_____________________________________________________________________________
 Int_t AliAnalysisMuMuResult::NofTriggers() const
 {
+  /// Get the number of triggers
+  
   if ( !Mother() ) return fNofTriggers;
   else return Mother()->NofTriggers();
 }
@@ -1828,6 +2249,8 @@ Int_t AliAnalysisMuMuResult::NofTriggers() const
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::Print(Option_t* opt) const
 {
+  /// printout
+  
   TString sopt(opt);
   sopt.ToUpper();
   
@@ -1840,7 +2263,14 @@ void AliAnalysisMuMuResult::Print(Option_t* opt) const
   pot.ReplaceAll("ALL","");
   pot.ReplaceAll("FULL","");
 
-  std::cout << pot.Data() << Form("%50s- NRUNS %d - NTRIGGER %10d - %s%s",
+  std::cout << pot.Data();
+  
+  if ( fAlias.Length() > 0 )
+  {
+    std::cout << Form("%s - ",fAlias.Data());
+  }
+  
+  std::cout << Form("%50s - NRUNS %d - NTRIGGER %10d - %s%s",
                GetName(),
                     NofRuns(),
                     NofTriggers(),
@@ -1852,6 +2282,11 @@ void AliAnalysisMuMuResult::Print(Option_t* opt) const
       std::cout << " (" << fSubResults->GetEntries()-1 << " subresults)";
     }
   
+  if (!fSubResults )
+  {
+    std::cout << Form(" - REBIN %d",fRebin) << std::endl;
+  }
+  
   std::cout << std::endl;
 
   if ( sopt.Contains("DUMP") )
@@ -1881,7 +2316,7 @@ void AliAnalysisMuMuResult::Print(Option_t* opt) const
     std::cout << Form("\t\tMBR %e +- %e",GetValue("MBR"),GetErrorStat("MBR")) << std::endl;
   }
   
-  if ( fSubResults && fSubResults->GetEntries() > 1 && ( sopt.Contains("ALL") || sopt.Contains("FULL") ) )
+  if ( fSubResults /* && fSubResults->GetEntries() > 1 */ && ( sopt.Contains("ALL") || sopt.Contains("FULL") ) )
   {
     std::cout << pot.Data() << "\t===== sub results ===== " << std::endl;
     
@@ -1900,6 +2335,8 @@ void AliAnalysisMuMuResult::Print(Option_t* opt) const
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::PrintParticle(const char* particle, const char* opt) const
 {
+  /// Print all information about one particule type
+  
   Double_t npart = GetValue(Form("Nof%s",particle));
   if (npart<=0) return;
   
@@ -1921,7 +2358,7 @@ void AliAnalysisMuMuResult::PrintParticle(const char* particle, const char* opt)
 }
 
 //_____________________________________________________________________________
-void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat) const
+void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat)
 {
   // print one value and its associated error
   
@@ -1933,7 +2370,7 @@ void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_
   {
     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f MeV/c^2",key,value*1E3,1E3*errorStat) << std::endl;
   }
-  else if ( TString(key).BeginsWith("Nof") )
+  else if ( TString(key).Contains("Nof") )
   {
     std::cout << opt << Form("\t\t%20s %9.2f +- %5.2f",key,value,errorStat) << std::endl;
   }
@@ -1950,6 +2387,8 @@ void AliAnalysisMuMuResult::PrintValue(const char* key, const char* opt, Double_
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::Set(const char* name, Double_t value, Double_t errorStat)
 {
+  /// Set a (value,error) pair with a given name
+  
   if ( !fMap )
   {
     fMap = new TMap;
@@ -1988,6 +2427,8 @@ void AliAnalysisMuMuResult::Set(const char* name, Double_t value, Double_t error
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::SetBin(const AliAnalysisMuMuBinning::Range& bin)
 {
+  /// Set the bin
+  
   if (!Mother()) fBin = bin;
   else Mother()->SetBin(bin);
 }
@@ -1995,6 +2436,8 @@ void AliAnalysisMuMuResult::SetBin(const AliAnalysisMuMuBinning::Range& bin)
 //_____________________________________________________________________________
 void AliAnalysisMuMuResult::SetNofInputParticles(const TH1& hminv)
 {
+  /// Set the number of input particle from the invariant mass spectra
+  
   static const char* particleNames[] = { "Jpsi" , "PsiPrime", "Upsilon","UpsilonPrime" };
 
   for ( Int_t i = 0; i < 4; ++i )
@@ -2075,3 +2518,24 @@ void AliAnalysisMuMuResult::SetMinv(const TH1& hminv)
   fMinv->SetDirectory(0);
 }
 
+//_____________________________________________________________________________
+AliAnalysisMuMuResult*
+AliAnalysisMuMuResult::SubResult(const char* subResultName) const
+{
+  /// get a given subresult
+  if (!fSubResults)
+  {
+    return 0x0;
+  }
+  TIter next(fSubResults);
+  AliAnalysisMuMuResult* r;
+  while ( ( r = static_cast<AliAnalysisMuMuResult*>(next())) )
+  {
+    if ( r->Alias() == subResultName )
+    {
+      return r;
+    }
+  }
+  return 0x0;
+}
+
index a2d2d97..0b009a7 100644 (file)
@@ -42,7 +42,7 @@ public:
                         const AliAnalysisMuMuBinning::Range& bin);
   
   AliAnalysisMuMuResult(const AliAnalysisMuMuResult& rhs);
-  AliAnalysisMuMuResult& operator=(const AliAnalysisMuMuResult&);
+  AliAnalysisMuMuResult& operator=(const AliAnalysisMuMuResult& rhs);
   
   virtual ~AliAnalysisMuMuResult();
 
@@ -66,21 +66,17 @@ public:
   
   void Print(Option_t* opt="") const;
   
-  Bool_t AddFit(const char* fitType, Int_t rebin=1);
+  Bool_t AddFit(const char* fitType, Int_t npar=0, Double_t* par=0x0);
 
   AliAnalysisMuMuResult* CountJpsi(TH1& h);
 
-//  void FitJpsiPsiPrimeCB(TH1& h);
   AliAnalysisMuMuResult*  FitJpsi(TH1& h);
-//  void FitJpsiCBE(TH1& h);
-//  void FitJpsiECBE(TH1& h);
-//  void FitJpsiPCBE(TH1& h);
-//  void FitUpsilon(TH1& h);
 
-  AliAnalysisMuMuResult* FitJpsi2CB2VWG(const TH1& h);
-  AliAnalysisMuMuResult* FitJpsiGCBE(TH1& h);
+  AliAnalysisMuMuResult* FitJpsiNA48(const TH1& h);
+  AliAnalysisMuMuResult* FitJpsiCB2VWG(const TH1& h);
+  AliAnalysisMuMuResult* FitJpsi2CB2VWG(const TH1& h, Double_t alphaLow=-1.0, Double_t nLow=-1.0, Double_t alphaUp=-1.0, Double_t nUp=-1.0);
   
-//  SubResult* FitJpsiPsiPrimeCustom(TH1& h);
+  AliAnalysisMuMuResult* FitJpsiGCBE(TH1& h);
   
   Int_t NofRuns() const;
   
@@ -96,14 +92,10 @@ public:
 
   void SetMinv(const TH1& hminv);
 
+  AliAnalysisMuMuResult* SubResult(const char* subResultName) const;
+  
   TObjArray* SubResults() const { return fSubResults; }
   
-  static Double_t CountParticle(const TH1& hminv, const char* particle, Double_t sigma=-1);
-
-  static Double_t ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr);
-
-  static Double_t ErrorABC(Double_t a, Double_t aerr, Double_t b, Double_t berr, Double_t c, Double_t cerror);
-
   Long64_t Merge(TCollection* list);
 
   AliAnalysisMuMuResult* Mother() const { return fMother; }
@@ -113,6 +105,18 @@ public:
   Double_t Weight() const { return fWeight > 0  ? fWeight : fNofTriggers; }
   
   void SetWeight(Double_t w) { fWeight=w; }
+
+  static Double_t CountParticle(const TH1& hminv, const char* particle, Double_t sigma=-1.0);
+  
+  static Double_t ErrorAB(Double_t a, Double_t aerr, Double_t b, Double_t berr);
+  
+  static Double_t ErrorABC(Double_t a, Double_t aerr, Double_t b, Double_t berr, Double_t c, Double_t cerror);
+
+  static void PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat);
+
+  void SetAlias(const char* alias) { fAlias = alias; }
+  
+  TString Alias() const { if ( fAlias.Length()>0) return fAlias; else return GetName(); }
   
 private:
   
@@ -122,9 +126,7 @@ private:
     kErrorStat=1
   };
   
-
   void PrintParticle(const char* particle, const char* opt) const;
-  void PrintValue(const char* key, const char* opt, Double_t value, Double_t errorStat) const;
 
 private:
   Int_t fNofRuns; // number of runs used to get this result
@@ -136,8 +138,16 @@ private:
   AliAnalysisMuMuResult* fMother; // mother result
   mutable THashList* fKeys; //! keys we have in our internal map (or the one of our subresults)
   Double_t fWeight; // weight of this result (default 1.0)
+  Int_t fRebin; // rebin level of minv spectra
+  
+  TString fTriggerClass; // trigger class for this result
+  TString fEventSelection; // event selection for this result
+  TString fPairSelection; // pair selection for this result
+  TString fCentralitySelection; // centrality selection for this result
+
+  TString fAlias; // alias name
   
-  ClassDef(AliAnalysisMuMuResult,5) // a class to hold invariant mass analysis results (counts, yields, AccxEff, R_AB, etc...)
+  ClassDef(AliAnalysisMuMuResult,8) // a class to hold invariant mass analysis results (counts, yields, AccxEff, R_AB, etc...)
 };
 
 #endif
index 6646c03..a5aa8b8 100644 (file)
@@ -1,3 +1,26 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+// AliAnalysisMuMuSpectra : a class to encapsulate results from MuMu analysis
+//
+// Spectra can be merged and converted into histograms
+//
+// author: L. Aphecetche (Subatech)
+//
+
 #include "AliAnalysisMuMuSpectra.h"
 
 #include "AliLog.h"
@@ -16,12 +39,74 @@ TNamed(name,title),
 fBinning(0x0),
 fBins(0x0)
 {
+  // default ctor
+}
+
+//______________________________________________________________________________
+AliAnalysisMuMuSpectra::AliAnalysisMuMuSpectra(const AliAnalysisMuMuSpectra& rhs)
+: TNamed(rhs.GetName(),rhs.GetTitle()),
+fBinning(0x0),
+fBins(0x0)
+{
+  // copy ctor
+
+  if ( rhs.fBinning )
+  {
+    fBinning = new AliAnalysisMuMuBinning(*rhs.fBinning);
+  }
+
+  TIter next(rhs.fBins);
+  AliAnalysisMuMuBinning::Range* bin;
+  
+  while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
+  {
+    if (!fBins)
+    {
+      fBins = new TObjArray;
+      fBins->SetOwner(kTRUE);
+    }
+    fBins->Add(bin);
+  }
+}
+
+//______________________________________________________________________________
+AliAnalysisMuMuSpectra&
+AliAnalysisMuMuSpectra::operator=(const AliAnalysisMuMuSpectra& rhs)
+{
+  // assignment operator
+  
+  if (this==&rhs) return *this;
+  
+  delete fBinning;
+  fBinning = 0x0;
+  delete fBins;
+  fBins = 0x0;
+  
+  if ( rhs.fBinning )
+  {
+    fBinning = new AliAnalysisMuMuBinning(*rhs.fBinning);
+  }
+  
+  TIter next(rhs.fBins);
+  AliAnalysisMuMuBinning::Range* bin;
   
+  while ( ( bin = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
+  {
+    if (!fBins)
+    {
+      fBins = new TObjArray;
+      fBins->SetOwner(kTRUE);
+    }
+    fBins->Add(bin);
+  }
+  
+  return *this;
 }
 
 //______________________________________________________________________________
 AliAnalysisMuMuSpectra::~AliAnalysisMuMuSpectra()
 {
+  // dtor
   delete fBinning;
   delete fBins;
 }
@@ -30,6 +115,7 @@ AliAnalysisMuMuSpectra::~AliAnalysisMuMuSpectra()
 void AliAnalysisMuMuSpectra::AdoptResult(const AliAnalysisMuMuBinning::Range& bin,
                                          AliAnalysisMuMuResult* result)
 {
+  // adopt (i.e. we are becoming the owner) a result for a given bin
   if (!fBinning)
   {
     fBinning = new AliAnalysisMuMuBinning;
@@ -76,10 +162,10 @@ Bool_t AliAnalysisMuMuSpectra::Correct(const AliAnalysisMuMuSpectra& accEff, con
   {
     AliAnalysisMuMuResult* thisResult = static_cast<AliAnalysisMuMuResult*>(fBins->At(i));
     AliAnalysisMuMuResult* accResult = static_cast<AliAnalysisMuMuResult*>(bins->At(i));
-    AliInfoClass(Form("i=%d",i ));
-    StdoutToAliInfoClass(thisResult->Print("full");
-                         std::cout << "----" << std::endl;
-                         accResult->Print("full"));
+//    AliInfoClass(Form("i=%d",i ));
+//    StdoutToAliInfoClass(thisResult->Print("full");
+//                         std::cout << "----" << std::endl;
+//                         accResult->Print("full"));
     
     thisResult->Correct(*accResult,particle,subResultName);
   }
@@ -90,8 +176,60 @@ Bool_t AliAnalysisMuMuSpectra::Correct(const AliAnalysisMuMuSpectra& accEff, con
 }
 
 //______________________________________________________________________________
+AliAnalysisMuMuResult*
+AliAnalysisMuMuSpectra::GetResultForBin(const AliAnalysisMuMuBinning::Range& bin) const
+{
+  /// Get result for a given bin
+  /// Warning: this requires a loop on bins
+  
+  if ( IsEmpty() ) return 0x0;
+  
+  TObjArray* bins = fBinning->CreateBinObjArray();
+  
+  Int_t j(-1);
+  
+  StdoutToAliDebug(1,std::cout << "searching for "; bin.Print());
+  
+  for ( Int_t i = 0; i <= bins->GetLast() && j < 0 ; ++i )
+  {
+    AliAnalysisMuMuBinning::Range* b = static_cast<AliAnalysisMuMuBinning::Range*>(bins->At(i));
+
+    StdoutToAliDebug(1,b->Print(););
+    
+    if ( bin == *b )
+    {
+      j = i;
+    }
+  }
+  
+  delete bins;
+  
+  if (j>=0)
+  {
+    return static_cast<AliAnalysisMuMuResult*>(fBins->At(j));
+  }
+  else
+  {
+    StdoutToAliDebug(1,std::cout << "Could not find result for bin:" << std::endl; bin.Print(););
+  }
+  return 0x0;
+}
+
+//______________________________________________________________________________
+Bool_t AliAnalysisMuMuSpectra::HasValue(const char* what) const
+{
+    // whether or not our result(s) has a given property
+  if ( IsEmpty() ) return kFALSE;
+  
+  AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(fBins->First());
+  
+  return r->HasValue(what);
+}
+
+//______________________________________________________________________________
 Bool_t AliAnalysisMuMuSpectra::IsEmpty() const
 {
+  // whether this spectra is empty or not
   return ( fBins==0x0 || fBins->GetEntries()<=0 );
 }
 
@@ -155,8 +293,10 @@ Long64_t AliAnalysisMuMuSpectra::Merge(TCollection* list)
 }
 
 //_____________________________________________________________________________
-TH1* AliAnalysisMuMuSpectra::Plot(const char* what, const char* subresult) const
+TH1* AliAnalysisMuMuSpectra::Plot(const char* what, const char* subresult, Bool_t divideByBinWidth) const
 {
+  // Convert the spectra into an histogram
+  
   TString swhat(what);
   swhat.ToUpper();
   
@@ -166,27 +306,22 @@ TH1* AliAnalysisMuMuSpectra::Plot(const char* what, const char* subresult) const
   
   TH1* h(0x0);
   
-  AliDebugClass(1,Form("nbins=%d nresults=%d",binArray->GetEntries(),fBins->GetEntries()));
+  AliDebug(1,Form("nbins=%d nresults=%d",binArray->GetEntries(),fBins->GetEntries()));
   
   for ( Int_t j = 0; j < TMath::Min(binArray->GetEntries(),fBins->GetEntries()); ++j )
   {
     AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(fBins->At(j));
     
-    if ( strlen(subresult) > 0 )
+    if ( strlen(subresult) > 0 && r->SubResults() )
     {
-      TObjArray* sr = r->SubResults();
-      if (!sr) continue;
       TString sub(subresult);
       sub.ToUpper();
-      r = static_cast<AliAnalysisMuMuResult*>(sr->FindObject(sub.Data()));
+      r = r->SubResult(sub.Data());
       if (!r) continue;
     }
     
     const AliAnalysisMuMuBinning::Range& b = r->Bin();
     
-    r->Print();
-    b.Print();
-    
     if (!h)
     {
       h = new TH1F(r->GetName(),r->GetName(),binArray->GetEntries(),bins);
@@ -195,18 +330,19 @@ TH1* AliAnalysisMuMuSpectra::Plot(const char* what, const char* subresult) const
     
     Double_t y = r->GetValue(what);
     Double_t yerr = r->GetErrorStat(what);
-    
-//    if ( !swhat.BeginsWith("ACC") )
-      if ( swhat.Contains("NOF") )
+  
+    if ( divideByBinWidth && b.WidthX()>0 )
     {
       y /= (b.WidthX());
       yerr /= (b.WidthX());
     }
     
+    std::cout << b.AsString();
+    AliAnalysisMuMuResult::PrintValue(swhat.Data(),"",y,yerr);
+    
     h->SetBinContent(j+1,y);
     h->SetBinError(j+1,yerr);
     
-    AliInfoClass(Form("%e +- %e",y,yerr));
   }
   
   delete binArray;
@@ -219,21 +355,17 @@ TH1* AliAnalysisMuMuSpectra::Plot(const char* what, const char* subresult) const
 //______________________________________________________________________________
 void AliAnalysisMuMuSpectra::Print(Option_t* opt) const
 {
+  // printout
+  
   if (!IsEmpty())
   {
     TString sopt(opt);
-    sopt.ToUpper();
-    if ( sopt.Contains("BINNING") )
-    {
-      fBinning->Print(opt);
-    }
     Int_t nmax = sopt.Atoi();
     if ( nmax <= 0 ) nmax = fBins->GetEntries();
     for ( Int_t i = 0; i < nmax; ++i )
     {
-      AliAnalysisMuMuBinning::Range* r = static_cast<AliAnalysisMuMuBinning::Range*>(fBins->At(i));
+      AliAnalysisMuMuResult* r = static_cast<AliAnalysisMuMuResult*>(fBins->At(i));
       if (r) r->Print(opt);
     }
-    
   }
 }
index 361749e..3080d82 100644 (file)
@@ -24,6 +24,9 @@ class AliAnalysisMuMuSpectra : public TNamed
 {
 public:
   AliAnalysisMuMuSpectra(const char* name="", const char* title="");
+  AliAnalysisMuMuSpectra(const AliAnalysisMuMuSpectra& rhs);
+  AliAnalysisMuMuSpectra& operator=(const AliAnalysisMuMuSpectra& rhs);
+  
   virtual ~AliAnalysisMuMuSpectra();
 
   void AdoptResult(const AliAnalysisMuMuBinning::Range& bin, AliAnalysisMuMuResult* result);
@@ -32,7 +35,7 @@ public:
   
   Long64_t Merge(TCollection* list);
 
-  TH1* Plot(const char* what="NofJpsi", const char* subresult="") const;
+  TH1* Plot(const char* what="NofJpsi", const char* subresult="", Bool_t divideByBinWidth=kTRUE) const;
 
   void Print(Option_t* opt="") const;
 
@@ -42,6 +45,10 @@ public:
   
   Bool_t Correct(const AliAnalysisMuMuSpectra& accEff, const char* particle, const char* subResultName="");
   
+  AliAnalysisMuMuResult* GetResultForBin(const AliAnalysisMuMuBinning::Range& bin) const;
+  
+  Bool_t HasValue(const char* what="NofJpsi") const;
+  
 private:
   AliAnalysisMuMuBinning* fBinning; // internal binning
   TObjArray* fBins; // the results (bin by bin)
index 98b7234..461e667 100644 (file)
@@ -537,6 +537,8 @@ AliAnalysisTriggerScalers::GetTriggerScaler(Int_t runNumber, const char* level,
   AliTriggerConfiguration* tc = static_cast<AliTriggerConfiguration*>(GetOCDBObject("GRP/CTP/Config",runNumber));
   AliTriggerRunScalers* trs = static_cast<AliTriggerRunScalers*>(GetOCDBObject("GRP/CTP/Scalers",runNumber));
   AliGRPObject* grp = static_cast<AliGRPObject*>(GetOCDBObject("GRP/GRP/Data",runNumber));
+
+  if (!tc || !trs || !grp) return 0x0;
   
   TString diCurrent(Form("L3:%5.0f;DIP:%5.0f [L3:%d;DIP:%d]",
                          grp->GetL3Current((AliGRPObject::Stats)0),
@@ -544,8 +546,6 @@ AliAnalysisTriggerScalers::GetTriggerScaler(Int_t runNumber, const char* level,
                          grp->GetL3Polarity(),
                          grp->GetDipolePolarity()));
   
-  if (!tc || !trs || !grp) return 0x0;
-  
   const TObjArray& trClasses = tc->GetClasses();
   
   const TObjArray* scalers = trs->GetScalersRecords();  
@@ -1002,6 +1002,7 @@ void AliAnalysisTriggerScalers::IntegratedLuminosity(const char* triggerList,
     if (!lumiB)
     {
       AliError(Form("Did not find lumiTrigger %s for run %09d",lumiTriggerClassName.Data(),runNumber));
+      continue;
     }
         
     Float_t pacCorrection(1.0);
@@ -1260,6 +1261,7 @@ TGraph* AliAnalysisTriggerScalers::PlotTriggerEvolution(const char* triggerClass
   /// - mu ( = -TMath::Log( 1 - P(0) ) where P(0) is the proba to have zero collisions in a bunch crossing)
   /// - pileupfactor = mu/(1-exp(-mu)) : the factor to apply to correct the cint1b count rate
   /// - vsnb = L0B/(NumberOfInteractingBunches*11245)
+  /// - NBCX = NumberOfInteractingBunches
   
   TString swhat(what);
   swhat.ToUpper();
@@ -1389,7 +1391,12 @@ TGraph* AliAnalysisTriggerScalers::PlotTriggerEvolution(const char* triggerClass
         UInt_t l2a = scaler->GetL2CA() - refl2a;
         UInt_t timelapse = seconds - reft;
         
-        if ( l0b <= 2 || ( l0a <= 2 && l0a != 0 ) || timelapse <= 9 ) continue;
+        if ( removeZeros && ( l0b <= 2 || ( l0a <= 2 && l0a != 0 ) || timelapse <= 9 ) )
+        {
+          AliInfo(Form("Skipping point for %s l0b %d l0a %d timelapse %d ts %s",
+                       triggerClassName,l0b,l0a,timelapse,ts.AsString()));
+          continue;
+        }
         
         reft = seconds;
         refl0b = scaler->GetLOCB();
@@ -1460,6 +1467,11 @@ TGraph* AliAnalysisTriggerScalers::PlotTriggerEvolution(const char* triggerClass
           value = l0b/(11245.0*numberOfInteractingBunches);
           error = -1.0; // FIXME
         }
+        else if ( swhat.Contains("NBCX"))
+        {
+          value = numberOfInteractingBunches;
+          error = 1.0; // FIXME          
+        }
         else
         {
           value = timelapse;
@@ -1468,7 +1480,7 @@ TGraph* AliAnalysisTriggerScalers::PlotTriggerEvolution(const char* triggerClass
         
         if ( ! swhat.Contains("OVER") && ! swhat.Contains("RATIO") &&
             ! swhat.Contains("MU") && ! swhat.Contains("PILEUPFACTOR") &&
-            ! swhat.Contains("RAW") )
+            ! swhat.Contains("RAW") & ! swhat.Contains("NBCX") )
         {
           value /= timelapse;
         }
@@ -1584,8 +1596,11 @@ TGraph* AliAnalysisTriggerScalers::PlotTriggerRatioEvolution(const char* trigger
 {
   /// Plots the evolution of one trigger ratio
   
-  TGraph* g1 = PlotTriggerEvolution(triggerClassName1,what1,kFALSE);
-  TGraph* g2 = PlotTriggerEvolution(triggerClassName2,what2,kFALSE);
+  Bool_t draw(kFALSE);
+  Bool_t removeZeros(kFALSE);
+  
+  TGraph* g1 = PlotTriggerEvolution(triggerClassName1,what1,draw,0x0,removeZeros);
+  TGraph* g2 = PlotTriggerEvolution(triggerClassName2,what2,draw,0x0,removeZeros);
   
   if (!g1 || !g2) return 0x0;
   
index e9121b4..1f84786 100644 (file)
@@ -1,3 +1,58 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//
+// AliMuonAccEffSubmitter : a class to help submit Acc x Eff simulations
+// anchored to real runs for J/psi, upsilon, single muons, etc...
+//
+// This class is dealing with 3 different directories :
+//
+// - template directory ($ALICE_ROOT/PWG/muondep/AccEffTemplates) containing the
+//   basic template files to be used for a simuation. A template can contain
+//   some variables that will be replaced during during the copy from template
+//   to local dir
+//
+// - local directory, where the files from the template directory, are copied
+//   once the class has been configured properly (i.e. using the various Set, Use,
+//   etc... methods). Some other files (e.g. JDL ones) are generated from
+//   scratch and also copied into this directory.
+//   At this point one could(should) check the files, as they are the ones
+//   to be copied to the remote directory for the production
+//
+// - remote directory, the alien directory where the files will be copied
+//   (from the local directory) before the actual submission
+//
+// ==========================================================
+//
+// Basic usage
+//
+// AliMuonAccEffSubmitter a;
+// a.UseOCDBSnapshots(kFALSE);
+// a.SetRemoteDir("/alice/cern.ch/user/l/laphecet/Analysis/LHC13d/simjpsi/pp503z0");
+// a.ShouldOverwriteFiles(true);
+// a.MakeNofEventsPropToTriggerCount("CMUL7-B-NOPF-MUON");
+// a.SetVar("VAR_GENLIB_PARNAME","\"pp 5.03\"");
+// a.SetRunList(195682);
+// a.Print();
+// a.Run("test"); // will do everything but the submit
+// a.Submit(false); // actual submission
+//
+//
+// author: Laurent Aphecetche (Subatech
+//
+
 #include "AliMuonAccEffSubmitter.h"
 
 #include "AliAnalysisTriggerScalers.h"
@@ -8,6 +63,7 @@
 #include "TMap.h"
 #include "TMath.h"
 #include "TObjString.h"
+#include "TROOT.h"
 #include "TString.h"
 #include "TSystem.h"
 #include <vector>
@@ -18,12 +74,13 @@ namespace
 }
 
 //______________________________________________________________________________
-AliMuonAccEffSubmitter::AliMuonAccEffSubmitter()
+AliMuonAccEffSubmitter::AliMuonAccEffSubmitter(const char* generator)
 : TObject(),
 fScalers(0x0),
 fRemoteDir(""),
 fReferenceTrigger(""),
 fRatio(1.0),
+fFixedNofEvents(10000),
 fMaxEventsPerChunk(5000),
 fLocalDir(gSystem->pwd()),
 fOCDBPath("raw://"),
@@ -52,13 +109,30 @@ fSnapshotDir(fLocalDir)
     fIsValid = kFALSE;
   }
 
-  SetPackages("VO_ALICE@AliRoot::v5-03-Rev-09","VO_ALICE@GEANT3::v1-14-6","VO_ALICE@ROOT::v5-34-02-1");
+  SetPackages("VO_ALICE@AliRoot::v5-03-Rev-18","VO_ALICE@GEANT3::v1-14-8","VO_ALICE@ROOT::v5-34-05-1");
   
-//  SetVar("VAR_ENERGY","5.03");
-  SetVar("VAR_GENLIB_TYPE","AliGenMUONlib::kJpsi");
-  SetVar("VAR_GENLIB_PARNAME","\"pPb 5.03\"");
   SetVar("VAR_OCDB_PATH","\"raw://\"");
+
+  SetVar("VAR_GENPARAM_GENLIB_TYPE","AliGenMUONlib::kJpsi");
+  SetVar("VAR_GENPARAM_GENLIB_PARNAME","\"pPb 5.03\"");
+
+  SetVar("VAR_GENCORRHF_QUARK","5");
+  SetVar("VAR_GENCORRHF_ENERGY","5");
+
+  SetVar("VAR_GENPARAMCUSTOM_PDGPARTICLECODE","443");
+
+  // default values below are from J/psi p+Pb (from muon_calo pass)
+  SetVar("VAR_GENPARAMCUSTOM_Y_P0","4.08E5");
+  SetVar("VAR_GENPARAMCUSTOM_Y_P1","7.1E4");
+  
+  SetVar("VAR_GENPARAMCUSTOM_PT_P0","1.13E9");
+  SetVar("VAR_GENPARAMCUSTOM_PT_P1","18.05");
+  SetVar("VAR_GENPARAMCUSTOM_PT_P2","2.05");
+  SetVar("VAR_GENPARAMCUSTOM_PT_P3","3.34");
+
   UseOCDBSnapshots(kTRUE);
+  
+  SetGenerator(generator);
 }
 
 //______________________________________________________________________________
@@ -72,6 +146,36 @@ AliMuonAccEffSubmitter::~AliMuonAccEffSubmitter()
 }
 
 //______________________________________________________________________________
+Bool_t AliMuonAccEffSubmitter::CheckCompilation(const char* file) const
+{
+  /// Check whether file can be compiled or not
+  /// FIXME: use gSystem->TempFileName for tmpfile !
+  
+  Bool_t rv(kTRUE);
+  
+  TString sfile(gSystem->BaseName(file));
+  TString tmpfile(Form("tmpfile_%s",sfile.Data()));
+  
+  gSystem->Exec(Form("cp %s %s",file,tmpfile.Data()));
+  
+  ReplaceVars(tmpfile.Data());
+  
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+  gSystem->AddIncludePath("-I$ALICE_ROOT/EVGEN");
+
+  if (gROOT->LoadMacro(Form("%s++",tmpfile.Data())))
+  {
+    AliError(Form("macro %s can not be compiled. Please check.",file));
+    rv = kFALSE;
+  }
+  
+  gSystem->Exec(Form("rm %s",tmpfile.Data()));
+  
+  return rv;
+}
+
+
+//______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::CheckLocal() const
 {
   /// Check whether all required local files are there
@@ -124,6 +228,7 @@ void AliMuonAccEffSubmitter::CleanRemote() const
 //______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::CopyFile(const char* localFile)
 {
+  /// copy a local file to remote destination
   TString local;
   
   if ( gSystem->IsAbsoluteFileName(localFile) )
@@ -168,6 +273,7 @@ Bool_t AliMuonAccEffSubmitter::CopyFile(const char* localFile)
   
   if ( ok )
   {
+    AliDebug(1,Form("cp %s alien://%s",local.Data(),remote.Data()));
     return TFile::Cp(local.Data(),Form("alien://%s",remote.Data()));
   }
   else
@@ -213,6 +319,8 @@ Bool_t AliMuonAccEffSubmitter::CopyLocalFilesToRemote()
   
   if (!IsValid()) return kFALSE;
   
+  AliDebug(1,"");
+  
   if ( CheckRemoteDir() )
   {
     TString sdir(gSystem->ExpandPathName(LocalDir()));
@@ -239,6 +347,8 @@ Bool_t AliMuonAccEffSubmitter::CopyTemplateFilesToLocal()
   
   if (!IsValid()) return kFALSE;
 
+  AliDebug(1,"");
+
   TIter next(TemplateFileList());
   TObjString* file;
   
@@ -313,6 +423,9 @@ Bool_t AliMuonAccEffSubmitter::CopyTemplateFilesToLocal()
 //______________________________________________________________________________
 std::ostream* AliMuonAccEffSubmitter::CreateJDLFile(const char* name) const
 {
+  /// Create a JDL file
+  AliDebug(1,"");
+
   TString jdl(Form("%s/%s",fLocalDir.Data(),name));
   
   if ( !ShouldOverwriteFiles() && !gSystem->AccessPathName(jdl.Data()) )
@@ -336,6 +449,11 @@ std::ostream* AliMuonAccEffSubmitter::CreateJDLFile(const char* name) const
 ///______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::GenerateMergeJDL(const char* name)
 {
+  /// Create the JDL for merging jobs
+  /// FIXME: not checked !
+  
+  AliDebug(1,"");
+
   std::ostream* os = CreateJDLFile(name);
   
   if (!os)
@@ -405,6 +523,8 @@ Bool_t AliMuonAccEffSubmitter::GenerateRunJDL(const char* name)
   /// Generate (locally) the JDL to perform the simulation+reco+aod filtering
   /// (to be then copied to the grid and finally submitted)
   
+  AliDebug(1,"");
+
   std::ostream* os = CreateJDLFile(name);
   
   if (!os)
@@ -482,6 +602,9 @@ Bool_t AliMuonAccEffSubmitter::GenerateRunJDL(const char* name)
 //______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::GetLastStage(const char* remoteDir) const
 {
+  /// Get the last staging phase already performed
+  /// FIXME : not checked !
+  
   Int_t n = 0, lastStage = 0;
   gSystem->Exec(Form("alien_ls -F %s | grep Stage_.*/ > __stage__", remoteDir));
   ifstream f("__stage__");
@@ -494,6 +617,45 @@ Bool_t AliMuonAccEffSubmitter::GetLastStage(const char* remoteDir) const
 }
 
 //______________________________________________________________________________
+TObjArray* AliMuonAccEffSubmitter::GetVariables(const char* file) const
+{
+  /// Find the variables in the file
+  
+  std::ifstream in(file);
+  char line[1024];
+  TObjArray* variables(0x0);
+  
+  while ( in.getline(line,1023,'\n') )
+  {
+    TString sline(line);
+    while (sline.Contains("VAR_") && !sline.BeginsWith("//") )
+    {
+      Int_t i1 = sline.Index("VAR_");
+      Int_t i2(i1);
+      
+      while ( ( i2 < sline.Length() ) && ( isalnum(sline[i2]) || sline[i2]=='_' ) ) ++i2;
+      
+      if (!variables)
+      {
+        variables = new TObjArray;
+        variables->SetOwner(kTRUE);
+      }
+      
+      TString var = sline(i1,i2-i1);
+      if ( !variables->FindObject(var) )
+      {
+        variables->Add(new TObjString(var));
+      }
+      sline.ReplaceAll(var,"");
+    }
+  }
+  
+  in.close();
+  
+  return variables;
+}
+
+//______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::HasVars(const char* file) const
 {
   /// Whether or not the file contains variables that have to
@@ -538,6 +700,8 @@ Bool_t AliMuonAccEffSubmitter::MakeOCDBSnapshots()
   
   if (!fScalers) return kFALSE;
   
+  AliDebug(1,"");
+
   const std::vector<int>& runs = fScalers->GetRunList();
 
   Bool_t ok(kTRUE);
@@ -549,10 +713,11 @@ Bool_t AliMuonAccEffSubmitter::MakeOCDBSnapshots()
     TString ocdbSim(Form("%s/OCDB/%d/OCDB_sim.root",SnapshotDir().Data(),runNumber));
     TString ocdbRec(Form("%s/OCDB/%d/OCDB_rec.root",SnapshotDir().Data(),runNumber));
 
-    if ( !gSystem->AccessPathName(ocdbSim.Data()) && 
+    if ( !gSystem->AccessPathName(ocdbSim.Data()) &&
          !gSystem->AccessPathName(ocdbRec.Data()) )
     {
       AliWarning(Form("Local OCDB snapshots already there for run %d. Will not redo them. If you want to force them, delete them by hand !",runNumber));
+      continue;
     }
     else
     {
@@ -754,6 +919,8 @@ UInt_t AliMuonAccEffSubmitter::NofRuns() const
 void AliMuonAccEffSubmitter::Output(std::ostream& out, const char* key,
                                     const TObjArray& values) const
 {
+  /// output to ostream of key,{values} pair
+  
   out << key << " = ";
   
   Int_t n = values.GetEntries();
@@ -795,6 +962,8 @@ void AliMuonAccEffSubmitter::Output(std::ostream& out, const char* key, const ch
                                     const char* v5, const char* v6, const char* v7,
                                     const char* v8, const char* v9) const
 {
+  /// output to ostream
+  
   TObjArray values;
   values.SetOwner(kTRUE);
   
@@ -815,6 +984,8 @@ void AliMuonAccEffSubmitter::Output(std::ostream& out, const char* key, const ch
 //______________________________________________________________________________
 void AliMuonAccEffSubmitter::Print(Option_t* /*opt*/) const
 {
+  /// Printout
+  
   if (!IsValid())
   {
     std::cout << std::string(80,'*') << std::endl;
@@ -822,9 +993,14 @@ void AliMuonAccEffSubmitter::Print(Option_t* /*opt*/) const
     std::cout << std::string(80,'*') << std::endl;
   }
     
-  std::cout << "Template directory = " << fTemplateDir.Data() << std::endl;
-  std::cout << "Local    directory = " << fLocalDir.Data() << std::endl;
-  std::cout << "Remote   directory = " << fRemoteDir.Data() << std::endl;
+  std::cout << "Template  directory = " << fTemplateDir.Data() << std::endl;
+  std::cout << "Local     directory = " << fLocalDir.Data() << std::endl;
+  std::cout << "Remote    directory = " << fRemoteDir.Data() << std::endl;
+  
+  if ( fSnapshotDir != fLocalDir )
+  {
+    std::cout << "Snapshots directory = " << fSnapshotDir.Data() << std::endl;
+  }
   
   std::cout << "OCDB path = " << fOCDBPath.Data() << std::endl;
   
@@ -920,8 +1096,10 @@ Bool_t AliMuonAccEffSubmitter::RemoteFileExists(const char *lfn) const
 }
 
 //______________________________________________________________________________
-Bool_t AliMuonAccEffSubmitter::ReplaceVars(const char* file)
+Bool_t AliMuonAccEffSubmitter::ReplaceVars(const char* file) const
 {
+  /// Replace the variables (i.e. things starting by VAR_) found in file
+  
   std::ifstream in(file);
   char line[1024];
   TObjArray lines;
@@ -1053,6 +1231,9 @@ void AliMuonAccEffSubmitter::SetPackages(const char* aliroot,
                                          const char* geant3,
                                          const char* api)
 {
+  /// Set the packages to be used by the jobs
+  /// Must be a valid combination, see http://alimonitor.cern.ch/packages/
+  ///
   fPackageAliroot = aliroot;
   fPackageRoot = root;
   fPackageGeant3 = geant3;
@@ -1085,8 +1266,63 @@ TString AliMuonAccEffSubmitter::GetRemoteDir(const char* dir, Bool_t create)
 }
 
 //______________________________________________________________________________
+Bool_t AliMuonAccEffSubmitter::SetGenerator(const char* generator)
+{
+  // set the variable to select the generator macro in Config.C
+  
+  gSystem->Load("libEVGEN");
+    
+  fIsValid = kFALSE;
+  
+  TString generatorFile(Form("%s/%s.C",fTemplateDir.Data(),generator));
+  
+  Int_t nofMissingVariables(0);
+  
+  // first check we indeed have such a macro
+  if (!gSystem->AccessPathName(generatorFile.Data()))
+  {
+    TObjArray* variables = GetVariables(generatorFile.Data());
+    
+    TIter next(variables);
+    TObjString* var;
+    
+    while ( ( var = static_cast<TObjString*>(next())) )
+    {
+      if ( !fVars->GetValue(var->String()) )
+      {
+        ++nofMissingVariables;
+        AliError(Form("file %s expect the variable %s to be defined, but we've not defined it !",generatorFile.Data(),var->String().Data()));
+      }
+    }
+    
+    delete variables;
+    
+    if ( !nofMissingVariables )
+    {
+      if (CheckCompilation(generatorFile.Data()))
+      {
+        fIsValid = kTRUE;
+        SetVar("VAR_GENERATOR",Form("%s",generator));        
+        TemplateFileList()->Add(new TObjString(Form("%s.C",generator)));
+        return kTRUE;
+      }
+    }
+    else
+    {
+      return kFALSE;
+    }
+  }
+  else
+  {
+    AliError(Form("Can not work with the macro %s",generatorFile.Data()));
+  }
+  return kFALSE;
+}
+
+//______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::SetMergedDir(const char* dir, Bool_t create)
 {
+  // Set the merged directory to be used
   fMergedDir = GetRemoteDir(dir,create);
   return (fMergedDir.Length()>0);
 }
@@ -1094,6 +1330,7 @@ Bool_t AliMuonAccEffSubmitter::SetMergedDir(const char* dir, Bool_t create)
 //______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::SetRemoteDir(const char* dir, Bool_t create)
 {
+  // Set the remote directory to be used
   fRemoteDir = GetRemoteDir(dir,create);
   return (fIsValid = (fRemoteDir.Length()>0));
 }
@@ -1130,8 +1367,44 @@ void AliMuonAccEffSubmitter::SetRunList(int runNumber)
 }
 
 //______________________________________________________________________________
+void AliMuonAccEffSubmitter::SetOCDBPath(const char* ocdbPath)
+{
+  /// Sets the OCDB path to be used
+  
+  fOCDBPath = ocdbPath;
+  
+  if (fScalers)
+  {
+    // redefine trigger scalers to use the new ocdb path
+    AliAnalysisTriggerScalers* ts = new AliAnalysisTriggerScalers(fScalers->GetRunList(),
+                                                                  fOCDBPath.Data());
+    
+    delete fScalers;
+    fScalers = ts;
+  }
+}
+
+
+//______________________________________________________________________________
+void AliMuonAccEffSubmitter::SetOCDBSnapshotDir(const char* dir)
+{
+  // change the directory used for snapshot
+  
+  if (gSystem->AccessPathName(Form("%s/OCDB",dir)))
+  {
+    AliError(Form("Snapshot top directory (%s) should contain an OCDB subdir with runnumbers in there",dir));
+  }
+  else
+  {
+    fSnapshotDir = dir;
+  }
+}
+
+//______________________________________________________________________________
 Bool_t AliMuonAccEffSubmitter::SetVar(const char* varname, const char* value)
 {
+  /// Set a variable
+  
   TString s(varname);
   s.ToUpper();
   if (!s.BeginsWith("VAR_"))
@@ -1169,6 +1442,8 @@ Int_t AliMuonAccEffSubmitter::Submit(Bool_t dryRun)
   
   if (!IsValid()) return 0;
   
+  AliDebug(1,"");
+
   gGrid->Cd(RemoteDir());
   
   if (!RemoteFileExists(RunJDLName()))
@@ -1316,6 +1591,8 @@ TObjArray* AliMuonAccEffSubmitter::TemplateFileList() const
 //______________________________________________________________________________
 void AliMuonAccEffSubmitter::UpdateLocalFileList(Bool_t clearSnapshots)
 {
+  /// Update the list of local files
+  
   if (!fScalers) return;
   
   if ( clearSnapshots )
@@ -1360,6 +1637,11 @@ void AliMuonAccEffSubmitter::UpdateLocalFileList(Bool_t clearSnapshots)
 //______________________________________________________________________________
 void AliMuonAccEffSubmitter::UseOCDBSnapshots(Bool_t flag)
 {
+  /// Whether or not to use OCDB snapshots
+  /// Using OCDB snapshots will speed-up both the sim and reco initialization
+  /// phases on each worker node, but takes time to produce...
+  /// So using them is not always a win-win...
+  
   fUseOCDBSnapshots = flag;
   if ( flag )
   {
index 3e75349..de37a75 100644 (file)
@@ -1,6 +1,13 @@
 #ifndef ALIMUONACCEFFSUBMITTER_H
 #define ALIMUONACCEFFSUBMITTER_H
 
+//
+// AliMuonAccEffSubmitter : a class to help submit Acc x Eff simulations
+// anchored to real runs for J/psi, upsilon, single muons, etc...
+//
+// author: Laurent Aphecetche (Subatech)
+//
+
 #include "TObject.h"
 #include "TString.h"
 #include "Riostream.h"
@@ -11,12 +18,12 @@ class TMap;
 class AliMuonAccEffSubmitter : public TObject
 {
 public:
-  AliMuonAccEffSubmitter();
+  AliMuonAccEffSubmitter(const char* generator="GenParamCustom");
 
   virtual ~AliMuonAccEffSubmitter();
 
   Bool_t SetRemoteDir(const char* dir, Bool_t create = kTRUE);
-  void SetLocalDir(const char* LocalDir) { fLocalDir = LocalDir; }
+  void SetLocalDir(const char* localdir) { fLocalDir = localdir; }
   Bool_t SetMergedDir(const char* dir, Bool_t create = kTRUE);
 
   void UseOCDBSnapshots(Bool_t flag);
@@ -41,8 +48,8 @@ public:
   UInt_t NofRuns() const;
   
   void SetRatio(Float_t ratio) { fRatio = ratio; }
-  void SetOCDBPath(const char* ocdbPath) { fOCDBPath = ocdbPath; }
-                   
+  void SetOCDBPath(const char* ocdbPath);
+  
   void MakeNofEventsPropToTriggerCount(const char* trigger="CMUL8-S-NOPF-MUON", Float_t ratio=1.0) { fReferenceTrigger = trigger; fRatio = ratio; }
   void MakeNofEventsFixed(Int_t nevents) { fFixedNofEvents = nevents; fReferenceTrigger=""; fRatio=0.0; }
   
@@ -92,8 +99,16 @@ public:
   
   void SetOCDBSnapshotDir(const char* dir);
 
+  Bool_t SetGenerator(const char* generator);
+  
+  TObjArray* GetVariables(const char* file) const;
+  
+  Bool_t IsValid() const { return fIsValid; }
+  
 private:
 
+  TString SnapshotDir() const { return fSnapshotDir; }
+  
   TString GetRemoteDir(const char* dir, Bool_t create);
 
   std::ostream* CreateJDLFile(const char* name) const;
@@ -113,45 +128,47 @@ private:
   
   void Output(std::ostream& out, const char* key, const TObjArray& values) const;
 
-  Bool_t ReplaceVars(const char* file);
+  Bool_t ReplaceVars(const char* file) const;
 
   TObjArray* TemplateFileList() const;
 
   TObjArray* LocalFileList() const;
   
-  Bool_t IsValid() const { return fIsValid; }
-  
   Bool_t HasVars(const char* localFile) const;
 
   void UpdateLocalFileList(Bool_t clearSnapshot=kFALSE);
   
-  TString SnapshotDir() const { return fSnapshotDir; }
+  Bool_t CheckCompilation(const char* file) const;
+  
+private:
+  AliMuonAccEffSubmitter(const AliMuonAccEffSubmitter& rhs);
+  AliMuonAccEffSubmitter& operator=(const AliMuonAccEffSubmitter& rhs);
   
 private:
-  AliAnalysisTriggerScalers* fScalers;
-  TString fRemoteDir;
-  TString fReferenceTrigger;
-  Float_t fRatio;
-  Int_t fFixedNofEvents;
-  Int_t fMaxEventsPerChunk;
-  TString fLocalDir;
-  TString fOCDBPath;
-  TString fTemplateDir;
-  TString fPackageAliroot;
-  TString fPackageGeant3;
-  TString fPackageRoot;
-  TString fPackageApi;
-  TString fMergedDir;
-  Int_t fSplitMaxInputFileNumber;
-  Int_t fCompactMode;
-  Bool_t fShouldOverwriteFiles;
-  TMap* fVars;
-  TString fExternalConfig;
-  Bool_t fUseOCDBSnapshots;
-  Bool_t fIsValid;
-  mutable TObjArray* fTemplateFileList;
-  mutable TObjArray* fLocalFileList;
-  TString fSnapshotDir;
+  AliAnalysisTriggerScalers* fScalers; // helper class used to handle the runlist and the scalers
+  TString fRemoteDir; // remote directory to used in alien
+  TString fReferenceTrigger; // reference trigger (if any) to be used to get the number of events to be used per run
+  Float_t fRatio; // ratio simulated events vs real events
+  Int_t fFixedNofEvents; // fixed number of events to be used per run
+  Int_t fMaxEventsPerChunk; // max events to generate per subjob
+  TString fLocalDir; // local directory
+  TString fOCDBPath; // OCDB path
+  TString fTemplateDir; // template directory
+  TString fPackageAliroot; // which aliroot package to use
+  TString fPackageGeant3; // which geant3 package to use
+  TString fPackageRoot; // which root package to use (for valid root,geant3,aliroot combinations see http://alimonitor.cern.ch/packages/)
+  TString fPackageApi; // which API package to use
+  TString fMergedDir; // merge directory
+  Int_t fSplitMaxInputFileNumber; // used for merging jdl
+  Int_t fCompactMode; // controls which outputs are kept (0=everything, 1=only aods)
+  Bool_t fShouldOverwriteFiles; // whether any copy (of template to local) is allowed to overwrite existing files
+  TMap* fVars; // map of the variables we can replace in template files
+  TString fExternalConfig; // path to an (optional) external config file
+  Bool_t fUseOCDBSnapshots; // whether to use OCDB snapshots or not
+  Bool_t fIsValid; // whether this object is valid (i.e. properly configured)
+  mutable TObjArray* fTemplateFileList; // list of template files
+  mutable TObjArray* fLocalFileList; // list of local files
+  TString fSnapshotDir; // directory for OCDB snapshots
   
   ClassDef(AliMuonAccEffSubmitter,1) // Helper class to submit AccxEff single particle simulations
 };