]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
- Added template class for fitting (work in progress).
authorpulvir <pulvir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Aug 2010 20:26:13 +0000 (20:26 +0000)
committerpulvir <pulvir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Aug 2010 20:26:13 +0000 (20:26 +0000)
- Updated the 2010 cuts after some comparison tests for bug fixing.
- Added and example macro to run with AliEn plugin
- Corrected some cuts
- Removed a bug in analysis task

21 files changed:
PWG2/PWG2resonancesLinkDef.h
PWG2/RESONANCES/AliRsnAnalysisPhi7TeV.cxx
PWG2/RESONANCES/AliRsnAnalysisPhi7TeV.h
PWG2/RESONANCES/AliRsnAnalysisSE.cxx
PWG2/RESONANCES/AliRsnAnalysisSE.h
PWG2/RESONANCES/AliRsnCutESD2010.cxx
PWG2/RESONANCES/AliRsnCutPrimaryVertex.cxx
PWG2/RESONANCES/AliRsnCutPrimaryVertex.h
PWG2/RESONANCES/AliRsnFitResult.cxx [new file with mode: 0644]
PWG2/RESONANCES/AliRsnFitResult.h [new file with mode: 0644]
PWG2/RESONANCES/AliRsnFunction.cxx
PWG2/RESONANCES/AliRsnFunction.h
PWG2/RESONANCES/AliRsnValue.cxx
PWG2/RESONANCES/AliRsnValue.h
PWG2/RESONANCES/macros/test/AddAnalysisTaskRsnTest.C
PWG2/RESONANCES/macros/test/ConfigTaskRsn2010.C
PWG2/RESONANCES/macros/test/ConfigTaskRsnTest.C
PWG2/RESONANCES/macros/test/PluginDataByRun.C [new file with mode: 0644]
PWG2/RESONANCES/macros/test/runAlienPlugin.C [new file with mode: 0644]
PWG2/RESONANCES/macros/test/runLocal.C
PWG2/libPWG2resonances.pkg

index 417ebdf29a1cfd90627032283a7008b99201d736..3248077911374910e237f4953c03d3300a55cd34 100644 (file)
@@ -43,5 +43,6 @@
 
 #pragma link C++ class AliRsnMonitorTrack+;
 #pragma link C++ class AliRsnAnalysisMonitorTask+;
+#pragma link C++ class AliRsnFitResult+;
 
 #endif
index 02630176d7b06e01d0d17105fa1590dcf43bbdcc..1ca7fd11d82ffcb33e9483228530ce24dc0f15a7 100644 (file)
@@ -31,6 +31,9 @@
 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name) :
   AliAnalysisTaskSE(name),
   fUseMC(kFALSE),
+  fCheckITS(kTRUE),
+  fCheckTPC(kTRUE),
+  fCheckTOF(kTRUE),
   fPDG(0),
   fCh(0),
   fIM(0.0),
@@ -71,6 +74,9 @@ AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name) :
 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) :
   AliAnalysisTaskSE(copy),
   fUseMC(copy.fUseMC),
+  fCheckITS(copy.fCheckITS),
+  fCheckTPC(copy.fCheckTPC),
+  fCheckTOF(copy.fCheckTOF),
   fPDG(0),
   fCh(0),
   fIM(0.0),
@@ -110,6 +116,9 @@ AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7
 //
 
   fUseMC = copy.fUseMC;
+  fCheckITS = copy.fCheckITS;
+  fCheckTPC = copy.fCheckTPC;
+  fCheckTOF = copy.fCheckTOF;
 
   fMaxVz   = copy.fMaxVz;
   fMaxITSband = copy.fMaxITSband;
@@ -259,9 +268,6 @@ Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
 
   static Int_t evNum = 0;
   evNum++;
-
-  // debug message
-  AliDebug(AliLog::kDebug + 1, Form("Event %d -- number of tracks = %d", evNum, esd->GetNumberOfTracks()));
   
   // get the best primary vertex:
   // first try the one with tracks
@@ -269,11 +275,13 @@ Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
   Double_t            vzTrk = 1000.0;
   Double_t            vzSPD = 1000.0;
+  Int_t               ncTrk = -1;
+  Int_t               ncSPD = -1;
+  if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
+  if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
   if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
   if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
-  AliDebug(AliLog::kDebug + 1, Form("Event %d -- vertex with tracks: contributors = %d, abs(vz) = %f", evNum, vTrk->GetNContributors(), vzTrk));
-  AliDebug(AliLog::kDebug + 1, Form("Event %d -- vertex with SPD,    contributors = %d, abs(vz) = %f", evNum, vSPD->GetNContributors(), vzSPD));
-  if(vTrk->GetNContributors() > 0)
+  if(vTrk && ncTrk > 0)
   {
     // fill the histograms
     fVertexX[0]->Fill(vTrk->GetXv());
@@ -286,7 +294,7 @@ Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
     else
       return kFarTracksPrimaryVertex;
   }
-  else if (vSPD->GetNContributors() > 0)
+  else if (vSPD && ncSPD > 0)
   {
     // fill the histograms
     fVertexX[1]->Fill(vSPD->GetXv());
@@ -311,18 +319,30 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
 // This function works with the ESD object
 //
 
+  static Int_t lastRun = -1;
+
   // ITS stuff #1 create the response function
   Bool_t isMC = (stack != 0x0);
   AliITSPIDResponse itsrsp(isMC);
 
   // TOF stuff #1: init OCDB
   Int_t run = esd->GetRunNumber();
-  AliCDBManager *cdb = AliCDBManager::Instance();
-  cdb->SetDefaultStorage("raw://");
-  cdb->SetRun(run);
+  if (run != lastRun)
+  {
+    cout << "Run = " << run << " -- LAST = " << lastRun << endl;
+    lastRun = run;
+    AliCDBManager *cdb = AliCDBManager::Instance();
+    cdb->SetDefaultStorage("raw://");
+    cdb->SetRun(run);
+    fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
+    fTOFcalib->Init();
+  }
+  //AliCDBManager *cdb = AliCDBManager::Instance();
+  //cdb->SetDefaultStorage("raw://");
+  //cdb->SetRun(run);
   // TOF stuff #2: init calibration
-  fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
-  fTOFcalib->Init();
+  //fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
+  //fTOFcalib->Init();
   // TOF stuff #3: calibrate
   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
@@ -349,8 +369,8 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
   // loop on all tracks
   ULong_t  status;
   Int_t    i, k, charge, npos = 0, nneg = 0, nITS;
-  Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofRel;
-  Bool_t   okTOF, okTrack, isTPC, isITSSA;
+  Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofDiff, tofRel;
+  Bool_t   okQuality, okTOF, okTPC, okITS, okTrack, isTPC, isITSSA;
   UChar_t  itsCluMap;
   for (i = 0; i < ntracks; i++)
   {
@@ -362,6 +382,7 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
     mom     = track->P();
     isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
     isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
+    //cout << isTPC << ' ' << isITSSA << endl;
     
     // accept only tracks which are TPC+ITS or ITS standalone
     if (!isTPC && !isITSSA) continue;
@@ -374,44 +395,95 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
     if (isTPC)
     {
       // check standard ESD cuts
-      if (!fESDtrackCutsTPC.IsSelected(track)) okTrack = kFALSE;
+      okQuality = fESDtrackCutsTPC.IsSelected(track);
+      //cout << "GLOBAL -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
+      if (!okQuality) continue;
       
       // check TPC dE/dx
-      tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kKaon));
-      if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
-      if (tpcNSigma > tpcMaxNSigma) okTrack = kFALSE;
+      if (fCheckTPC)
+      {
+        tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kKaon));
+        if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
+        okTPC = (tpcNSigma <= tpcMaxNSigma);
+        //cout << "ALT -- TPC    -- nsigma = " << tpcNSigma << ", max = " << tpcMaxNSigma << " --> " << (okTPC ? "OK" : "FAILED") << endl;
+        //cout << "ALTTPC -- " << fTPCpar[0] << ' ' << fTPCpar[1] << ' ' << fTPCpar[2] << ' ' << fTPCpar[3] << ' ' << fTPCpar[4] << endl;
+      }
+      else
+      {
+        okTPC = kTRUE;
+      }
       
       // check TOF (only if momentum is large than function asymptote and flags are OK)
-      okTOF = kTRUE;
-      if (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) && mom > TMath::Max(b1, b2))
+      if (fCheckTOF)
       {
-        track->GetIntegratedTimes(times);
-        tofTime  = (Double_t)track->GetTOFsignal();
-        tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kKaon], AliPID::ParticleMass(AliPID::kKaon));
-        tofRef   = times[AliPID::kKaon];
-        if (tofRef > 0.0)
+        if (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) && mom > TMath::Max(b1, b2))
         {
-          tofRel   = (tofTime - tofRef) / tofRef;
-          ymax     = a1 / (mom - b1) + c1;
-          ymin     = a2 / (mom - b2) + c2;
-          okTOF    = (tofRel >= ymin && tofRel <= ymax);
+          track->GetIntegratedTimes(times);
+          tofTime  = (Double_t)track->GetTOFsignal();
+          tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kKaon], AliPID::ParticleMass(AliPID::kKaon));
+          tofRef   = times[AliPID::kKaon];
+          if (tofRef > 0.0)
+          {
+            tofDiff  = (tofTime - tofRef);
+            tofRel   = (tofTime - tofRef) / tofRef;
+            ymax     = a1 / (mom - b1) + c1;
+            ymin     = a2 / (mom - b2) + c2;
+            okTOF    = (tofRel >= ymin && tofRel <= ymax);
+            //cout << "TOF    -- diff = " << tofDiff << ", rel diff = " << tofRel << ", range = " << ymin << " to " << ymax << ", sigma = " << tofSigma << " --> " << (okTOF ? "OK" : "FAILED") << endl;
+          }
+          else
+          {
+            okTOF = kTRUE;
+            //cout << "TOF    -- not checked due to ZERO reference time" << endl;
+          }
         }
+        else
+        {
+          okTOF = kTRUE;
+          //cout << "TOF    -- not checked because TOF pid absent" << endl;
+        }
+      }
+      else
+      {
+        okTOF = kTRUE;
       }
-      if (!okTOF) okTrack = kFALSE;
+      
+      okTrack = okQuality && okTPC && okTOF;
+      //cout << "GLOBAL -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
     }
     else
     {
       // check standard ESD cuts
-      if (!fESDtrackCutsITS.IsSelected(track)) okTrack = kFALSE;
+      okQuality = fESDtrackCutsITS.IsSelected(track);
+      //cout << "ITSSA  -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
+      if (!okQuality) continue;
       
       // check dE/dx
-      itsSignal = track->GetITSsignal();
-      itsCluMap = track->GetITSClusterMap();
-      nITS      = 0;
-      for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
-      if (nITS < 3) okTrack = kFALSE;; // track not good for PID
-      itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kKaon, nITS, kTRUE);
-      if (TMath::Abs(itsNSigma) > fMaxITSband) okTrack = kFALSE;
+      if (fCheckITS)
+      {
+        itsSignal = track->GetITSsignal();
+        itsCluMap = track->GetITSClusterMap();
+        nITS      = 0;
+        for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
+        if (nITS < 3) 
+        {
+          okITS = kFALSE;
+          //cout << "ITS    -- not checked due to too few PID clusters" << endl;
+        }
+        else
+        {
+          itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kKaon, nITS, kTRUE);
+          okITS = (TMath::Abs(itsNSigma) <= fMaxITSband);
+          //cout << "ITS    -- nsigma = " << itsNSigma << ", max = " << fMaxITSband << " --> " << (okITS ? "OK" : "FAILED") << endl;
+        }
+      }
+      else
+      {
+        okITS = kTRUE;
+      }
+      
+      okTrack = okQuality && okITS;
+      //cout << "ITSSA  -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
     }
     
     // skip tracks not passing cuts
@@ -439,7 +511,25 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
   neg.Set(nneg);
   itspos.Set(npos);
   itsneg.Set(nneg);
-
+  
+  /*
+  // fill unlike-sign
+  Int_t i1, i2;
+  for (i1 = 0; i1 < npos; i1++)
+    for (i2 = 0; i2 < nneg; i2++)
+      AddEntryFromESD(esd, i1, i2, itspos[i1], itsneg[i2], 0, stack);
+      
+  // fill ++ like-sign
+  for (i1 = 0; i1 < npos; i1++)
+    for (i2 = i1 + 1; i2 < npos; i2++)
+      AddEntryFromESD(esd, i1, i2, itspos[i1], itspos[i2], 1, 0x0);
+      
+  // fill -- like-sign
+  for (i1 = 0; i1 < nneg; i1++)
+    for (i2 = i1 + 1; i2 < nneg; i2++)
+      AddEntryFromESD(esd, i1, i2, itsneg[i1], itsneg[i2], -1, 0x0);
+      */
+      
   // loop on unlike-sign pairs to compute invariant mass signal
   Int_t           ip, in, lp, ln;
   AliPID          pid;
@@ -492,8 +582,8 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
       fITS[0] = itspos[ip];
       fITS[1] = itsneg[in];
 
-      if (fIM < 0.9 || fIM >  1.4) continue;
-      if (fPt < 0.0 || fPt > 20.0) continue;
+      //if (fIM < 0.9 || fIM >  1.4) continue;
+      //if (fPt < 0.0 || fPt > 20.0) continue;
       
       fRsnTreeComp->Fill();
     }
@@ -527,8 +617,8 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
       fITS[0] = itspos[i1];
       fITS[1] = itspos[i2];
 
-      if (fIM < 0.9 || fIM >  1.4) continue;
-      if (fPt < 0.0 || fPt > 20.0) continue;
+      //if (fIM < 0.9 || fIM >  1.4) continue;
+      //if (fPt < 0.0 || fPt > 20.0) continue;
       
       fRsnTreeComp->Fill();
     }
@@ -556,8 +646,8 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
       fITS[0] = itsneg[i1];
       fITS[1] = itsneg[i2];
 
-      if (fIM < 0.9 || fIM >  1.4) continue;
-      if (fPt < 0.0 || fPt > 20.0) continue;
+      //if (fIM < 0.9 || fIM >  1.4) continue;
+      //if (fPt < 0.0 || fPt > 20.0) continue;
       
       fRsnTreeComp->Fill();
     }
@@ -566,6 +656,63 @@ void AliRsnAnalysisPhi7TeV::ProcessESD
   PostData(1, fRsnTreeComp);
 }
 
+//__________________________________________________________________________________________________
+void AliRsnAnalysisPhi7TeV::AddEntryFromESD(AliESDEvent *esd, Int_t i1, Int_t i2, Int_t its1, Int_t its2, Short_t charge, AliStack *stack)
+{
+//
+// Add an entry to the output TTree computed from two tracks in the ESD
+//
+
+  // loop on unlike-sign pairs to compute invariant mass signal
+  AliPID          pid;
+  Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
+  Double_t        phimass = 1.019455;
+  TLorentzVector  v1, v2, vsum, vref;
+  
+  AliESDtrack *t1 = esd->GetTrack(i1);
+  AliESDtrack *t2 = esd->GetTrack(i2);
+
+  v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
+  v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
+  vsum = v1 + v2;
+  vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
+  
+  // if stack is present, search for true pair
+  fPDG = 0;
+  Int_t l1 = TMath::Abs(t1->GetLabel());
+  Int_t l2 = TMath::Abs(t2->GetLabel());
+  if (stack) 
+  {
+    TParticle *part1 = stack->Particle(l1);
+    TParticle *part2 = stack->Particle(l2);
+    if (part1 && part2)
+    {
+      if (part1->GetFirstMother() == part2->GetFirstMother())
+      {
+        if (part1->GetFirstMother() > 0)
+        {
+          TParticle *mum = stack->Particle(part1->GetFirstMother());
+          fPDG = mum->GetPdgCode();
+        }
+      }
+    }
+  }
+  
+  fPDG = TMath::Abs(fPDG);
+  fCh     = charge;
+  fIM     = (Float_t)vsum.M();
+  fPt     = (Float_t)vsum.Perp();
+  fEta    = (Float_t)vsum.Eta();
+  fY      = (Float_t)vref.Rapidity();
+  fITS[0] = its1;
+  fITS[1] = its2;
+
+  //if (fIM < 0.9 || fIM >  1.4) return;
+  //if (fPt < 0.0 || fPt > 20.0) return;
+      
+  fRsnTreeComp->Fill();
+}
+  
 //__________________________________________________________________________________________________
 void AliRsnAnalysisPhi7TeV::ProcessMC(AliStack *stack)
 {
index 1945586b767933238b61f258ab8941081436c83b..ff0739457ee9bf8a3d2789f52824d70f6c2ec1d6 100644 (file)
@@ -42,6 +42,9 @@ class AliRsnAnalysisPhi7TeV : public AliAnalysisTaskSE
     virtual ~AliRsnAnalysisPhi7TeV();
 
     void             SetUseMC(Bool_t yn = kTRUE) {fUseMC = yn;}
+    void             SetCheckITS(Bool_t yn = kTRUE) {fCheckITS = yn;}
+    void             SetCheckTPC(Bool_t yn = kTRUE) {fCheckTPC = yn;}
+    void             SetCheckTOF(Bool_t yn = kTRUE) {fCheckTOF = yn;}
     
     void             SetMaxVz(Double_t v)   {fMaxVz = v;}
     
@@ -70,8 +73,12 @@ class AliRsnAnalysisPhi7TeV : public AliAnalysisTaskSE
 
     void     ProcessESD(AliESDEvent *esd, AliStack *stack);
     void     ProcessMC(AliStack *stack);
+    void     AddEntryFromESD(AliESDEvent *esd, Int_t i1, Int_t i2, Int_t its1, Int_t its2, Short_t charge, AliStack *stack = 0x0);
 
     Bool_t   fUseMC;      // use MC or data?
+    Bool_t   fCheckITS;   // chec ITS PID?
+    Bool_t   fCheckTPC;   // chec TPC PID?
+    Bool_t   fCheckTOF;   // chec TOF PID?
     
     Short_t  fPDG;        // PDG code
     Short_t  fCh;         // control flag for like/unlike sign
index 6483face6d02256aed1faeb4b3b04b3b738beeb5..ad757f8a265bf4bd91413ce61a5e45dc0ace0a3c 100644 (file)
@@ -54,6 +54,16 @@ AliRsnAnalysisSE::AliRsnAnalysisSE(const AliRsnAnalysisSE& copy) :
   AliDebug(AliLog::kDebug+2,"->");
 }
 
+//_____________________________________________________________________________
+AliRsnAnalysisSE::~AliRsnAnalysisSE()
+{
+//
+// Destructor
+//
+
+  fOutList->Clear();
+}
+
 //_____________________________________________________________________________
 void AliRsnAnalysisSE::RsnUserCreateOutputObjects()
 {
@@ -116,7 +126,7 @@ void AliRsnAnalysisSE::RsnUserExec(Option_t*)
       if ((zeroEventPercent>fZeroEventPercentWarning)&&(fEntry>100))
         AliWarning(Form("%3.2f%% Events are with zero tracks (CurrentEvent=%d)!!!",zeroEventPercent,fEntry));
     }
-    return;
+    //return;
   }
 
   // if general event cuts are added to the task (recommended)
index a914092fa838a9017ace85bcca9fb11259fb7d65..d4112d4744609ccc370cf9a083550daefa471d02 100644 (file)
@@ -24,7 +24,7 @@ class AliRsnAnalysisSE : public AliRsnVAnalysisTaskSE
   public:
     AliRsnAnalysisSE(const char *name = "AliRsnAnalysisSE", Bool_t useKine = kFALSE);
     AliRsnAnalysisSE(const AliRsnAnalysisSE& copy);
-    virtual ~AliRsnAnalysisSE() {;};
+    virtual ~AliRsnAnalysisSE();
 
     // Implement this
     virtual void    RsnUserCreateOutputObjects();
index 096c96068d5c0790a7b2e2c2b1d549b3958f2c12..910bc42c4ff619b85a027469a95bc726478d323f 100644 (file)
@@ -207,7 +207,7 @@ Bool_t AliRsnCutESD2010::IsSelected(TObject *obj1, TObject* /*obj2*/)
   ULong_t  status;
   Int_t    k, nITS;
   Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofRel;
-  Bool_t   okTOF, isTPC, isITSSA;
+  Bool_t   okQuality, okTOF, okTPC, okITS, okTrack, isTPC, isITSSA;
   UChar_t  itsCluMap;
   
   // get commonly used variables
@@ -222,20 +222,27 @@ Bool_t AliRsnCutESD2010::IsSelected(TObject *obj1, TObject* /*obj2*/)
   if (isTPC)
   {
     // check standard ESD cuts
-    if (!fESDtrackCutsTPC.IsSelected(track)) return kFALSE;
-    
+    okQuality = fESDtrackCutsTPC.IsSelected(track);
+    //cout << "GLOBAL -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
+    if (!okQuality) return kFALSE;
+  
     // check TPC dE/dx
     if (fCheckTPC)
     {
       tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kKaon));
       if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
-      if (tpcNSigma > tpcMaxNSigma) return kFALSE;
+      okTPC = (tpcNSigma <= tpcMaxNSigma);
+      //cout << "RSN -- TPC    -- nsigma = " << tpcNSigma << ", max = " << tpcMaxNSigma << " --> " << (okTPC ? "OK" : "FAILED") << endl;
+      //cout << "RSNTPC -- " << fTPCpar[0] << ' ' << fTPCpar[1] << ' ' << fTPCpar[2] << ' ' << fTPCpar[3] << ' ' << fTPCpar[4] << endl;
+    }
+    else
+    {
+      okTPC = kTRUE;
     }
     
     // check TOF (only if momentum is large than function asymptote and flags are OK)
     if (fCheckTOF)
     {
-      okTOF = kTRUE;
       if (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) && mom > TMath::Max(b1, b2))
       {
         track->GetIntegratedTimes(times);
@@ -248,29 +255,63 @@ Bool_t AliRsnCutESD2010::IsSelected(TObject *obj1, TObject* /*obj2*/)
           ymax     = a1 / (mom - b1) + c1;
           ymin     = a2 / (mom - b2) + c2;
           okTOF    = (tofRel >= ymin && tofRel <= ymax);
+          //cout << "TOF    -- diff = " << tofDiff << ", rel diff = " << tofRel << ", range = " << ymin << " to " << ymax << ", sigma = " << tofSigma << " --> " << (okTOF ? "OK" : "FAILED") << endl;
+        }
+        else
+        {
+          okTOF = kTRUE;
+          //cout << "TOF    -- not checked due to ZERO reference time" << endl;
         }
       }
-      if (!okTOF) return kFALSE;
+      else
+      {
+        okTOF = kTRUE;
+        //cout << "TOF    -- not checked because TOF pid absent" << endl;
+      }
+    }
+    else
+    {
+      okTOF = kTRUE;
     }
+    
+    // combine checks
+    okTrack = okQuality && okTPC && okTOF;
+    //cout << "GLOBAL -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
   }
   else
   {
     // check standard ESD cuts
-    if (!fESDtrackCutsITS.IsSelected(track)) return kFALSE;
+    okQuality = fESDtrackCutsITS.IsSelected(track);
+    //cout << "ITSSA  -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
+    if (!okQuality) return kFALSE;
     
     // check dE/dx
     if (fCheckITS)
     {
-      if ((status & AliESDtrack::kITSpid)  != 0) return kFALSE;
       itsSignal = track->GetITSsignal();
       itsCluMap = track->GetITSClusterMap();
       nITS      = 0;
       for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
-      if (nITS < 3) return kFALSE; // track not good for PID
-      itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kKaon, nITS, kTRUE);
-      if (TMath::Abs(itsNSigma) > fMaxITSband) return kFALSE;
+      if (nITS < 3) 
+      {
+        okITS = kFALSE;
+        //cout << "ITS    -- not checked due to too few PID clusters" << endl;
+      }
+      else
+      {
+        itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kKaon, nITS, kTRUE);
+        okITS = (TMath::Abs(itsNSigma) <= fMaxITSband);
+        //cout << "ITS    -- nsigma = " << itsNSigma << ", max = " << fMaxITSband << " --> " << (okITS ? "OK" : "FAILED") << endl;
+      }
     }
+    else
+    {
+      okITS = kTRUE;
+    }
+    
+    okTrack = okQuality && okITS;
+    //cout << "ITSSA  -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
   }
   
-  return kTRUE;
+  return okTrack;
 }
index 8b8dc5393fd2de5d36031601c79e168a9854c304..7adfe357acccd4b9714855f768a6b5f433eb5c7d 100644 (file)
@@ -44,7 +44,7 @@ AliRsnCutPrimaryVertex::AliRsnCutPrimaryVertex
 //
 
   fMinD = 0.0;
-  fMaxD = maxVz;
+  fMaxD = maxVz + 1E-6;
 }
 
 //_________________________________________________________________________________________________
@@ -54,6 +54,9 @@ Bool_t AliRsnCutPrimaryVertex::IsSelected(TObject *obj1, TObject* /*obj2*/)
 // Cut checker
 //
 
+  static Int_t evNum = 0;
+  evNum++;
+  
   // retrieve ESD event
   AliRsnEvent *rsn = dynamic_cast<AliRsnEvent*>(obj1);
   if (!rsn) return kFALSE;
@@ -69,33 +72,42 @@ Bool_t AliRsnCutPrimaryVertex::IsSelected(TObject *obj1, TObject* /*obj2*/)
   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
   const AliESDVertex *vTPC  = esd->GetPrimaryVertexTPC();
+  Int_t               ncTrk = -1;
+  Int_t               ncSPD = -1;
+  Int_t               ncTPC = -1;
   Double_t            vzTrk = 2.0 * fMaxD;
   Double_t            vzSPD = 2.0 * fMaxD;
   Double_t            vzTPC = 2.0 * fMaxD;
   if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
   if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
   if (vTPC) vzTPC = TMath::Abs(vTPC->GetZv());
-  AliDebug(AliLog::kDebug + 1, Form("Vertex with tracks: contributors = %d, abs(vz) = %f", vTrk->GetNContributors(), vzTrk));
-  AliDebug(AliLog::kDebug + 1, Form("Vertex with SPD,    contributors = %d, abs(vz) = %f", vSPD->GetNContributors(), vzSPD));
-  AliDebug(AliLog::kDebug + 1, Form("Vertex with TPC,    contributors = %d, abs(vz) = %f", vTPC->GetNContributors(), vzTPC));
-  if(vTrk->GetStatus())
+  if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
+  if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
+  if (vTPC) ncTPC = (Int_t)vTPC->GetNContributors();
+  if(vTrk && ncTrk > 0)
   {
-    fCutValueI = vTrk->GetNContributors();
+    fCutValueI = ncTrk;
     fCutValueD = vzTrk;
   }
-  else if (vSPD->GetStatus())
+  else if (vSPD && ncSPD > 0)
   {
-    fCutValueI = vSPD->GetNContributors();
+    fCutValueI = ncSPD;
     fCutValueD = vzSPD;
   }
-  else if (vTPC->GetStatus() && fAcceptTPC)
+  else if (vTPC && fAcceptTPC && ncTPC > 0)
   {
-    fCutValueI = vTPC->GetNContributors();
+    fCutValueI = ncTPC;
     fCutValueD = vzTPC;
   }
   else
-    return kFALSE;
+  {
+    fCutValueI = -1;
+    fCutValueD = 2.0 * fMaxD;
+  }
+  
+  // output
+  Bool_t result = ((!OkRangeI()) && OkRangeD());
   
-  return OkRangeI() && OkRangeD();
+  return result;
 }
 
index 7606d6f7be39d99cffbb86491d78a8e96d6d1c5d..14599be88138a60662f03953a160de1bb553aea9 100644 (file)
@@ -31,7 +31,7 @@ class AliRsnCutPrimaryVertex : public AliRsnCut
   public:
 
     AliRsnCutPrimaryVertex();
-    AliRsnCutPrimaryVertex(const char *name, Double_t maxVz, Int_t minContributors = 0, Bool_t acceptTPC = kFALSE);
+    AliRsnCutPrimaryVertex(const char *name, Double_t maxVz, Int_t minContributors = 1, Bool_t acceptTPC = kFALSE);
     virtual ~AliRsnCutPrimaryVertex() {;};
 
     virtual Bool_t IsSelected(TObject *obj1, TObject *obj2 = 0x0);
diff --git a/PWG2/RESONANCES/AliRsnFitResult.cxx b/PWG2/RESONANCES/AliRsnFitResult.cxx
new file mode 100644 (file)
index 0000000..6bcffbe
--- /dev/null
@@ -0,0 +1,276 @@
+//
+//
+//
+
+#include <TF1.h>
+#include <TH1.h>
+#include <TMath.h>
+
+#include "AliLog.h"
+#include <AliRsnFitResult.h>
+
+ClassImp(AliRsnFitResult)
+
+//__________________________________________________________________________________________________
+AliRsnFitResult::AliRsnFitResult(const char *name, ESignalType sType, EBackgroundType bType) :
+  TNamed(name, ""),
+  fSignalType(sType),
+  fBackgroundType(bType)
+{
+//
+// Default constructor
+//
+
+  InitFunctions();
+}
+  
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::NormSquareRoot(Double_t *x, Double_t *par)
+{
+//
+// Computes a square-root like function normalized 
+// to the integral in the background range specified
+// for this object.
+// This is obtained by dividing the function value by
+// its integral analytically computed in that range.
+// Returns 0 in all cases where a floating point exception
+// could be raised by computing square root of a negative number.
+//
+
+  if (x[0] < par[1]) return 0.0;
+  
+  Double_t fcnVal = TMath::Sqrt(x[0] - par[1]);
+
+  Double_t fcnMax = fPeakRange[1] - par[1];
+  Double_t fcnMin = fPeakRange[0] - par[1];
+  Double_t fcnInt = (2.0 / 3.0) * ((fcnMax*TMath::Sqrt(fcnMax) - fcnMin*TMath::Sqrt(fcnMin)));
+  
+  return par[0] * fcnVal / fcnInt;
+}
+
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::NormLinear(Double_t *x, Double_t *par)
+{
+//
+// Computes a 1st order polynomial function normalized 
+// to the integral in the background range specified
+// for this object.
+// This is obtained by dividing the function value by
+// its integral analytically computed in that range.
+// Returns 0 in all cases where a floating point exception
+// could be raised by computing square root of a negative number.
+//
+  
+  Double_t fcnVal = 1.0 + x[0] * par[1];
+  Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + 0.5*par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]);
+  
+  return par[0] * fcnVal / fcnInt;
+}
+  
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::NormPoly2(Double_t *x, Double_t *par)
+{
+//
+// Computes a 2nd order polynomial function normalized 
+// to the integral in the background range specified
+// for this object.
+// This is obtained by dividing the function value by
+// its integral analytically computed in that range.
+// Returns 0 in all cases where a floating point exception
+// could be raised by computing square root of a negative number.
+//
+  
+  Double_t fcnVal = 1.0 + x[0] * par[1] + x[0]*x[0] * par[2];
+  Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0])/2.0 + par[2]*(fPeakRange[1]*fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]*fPeakRange[0])/3.0;
+  
+  return par[0] * fcnVal / fcnInt;
+}
+
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::NormBreitWigner(Double_t *x, Double_t *par)
+{
+//
+// Computes a Breit-Wigner function normalized 
+// to the integral in the background range specified
+// for this object.
+// This is obtained by dividing the function value by
+// its integral analytically computed in that range.
+// Returns 0 in all cases where a floating point exception
+// could be raised by computing square root of a negative number.
+//
+  
+  Double_t fcnVal = 1.0 / ((x[0] - par[1])*(x[0] - par[1]) + 0.25*par[2]*par[2]);
+  Double_t fcnInt = 2.0 / par[2] * (TMath::ATan((fPeakRange[1] - par[1]) / (0.5 * par[2])) - TMath::ATan((fPeakRange[0] - par[1]) / (0.5 * par[2])));
+  
+  return par[0] * fcnVal / fcnInt;
+}
+
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::NormGaus(Double_t *x, Double_t *par)
+{
+//
+// Computes a Gaussian function normalized 
+// to the integral in the background range specified
+// for this object.
+// This is obtained by dividing the function value by
+// its integral analytically computed in that range.
+// Returns 0 in all cases where a floating point exception
+// could be raised by computing square root of a negative number.
+//
+  
+  Double_t fcnVal      = TMath::Gaus(x[0], par[1], par[2], kTRUE);
+  Double_t fcnIntLeft  = 0.5 * TMath::Erf(((par[1] - fPeakRange[0]) / par[2]) / TMath::Sqrt(2.0));
+  Double_t fcnIntRight = 0.5 * TMath::Erf(((fPeakRange[1] - par[1]) / par[2]) / TMath::Sqrt(2.0));
+  Double_t fcnInt      = fcnIntLeft + fcnIntRight;
+  
+  return par[0] * fcnVal / fcnInt;
+}
+
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::Signal(Double_t *x, Double_t *par)
+{
+//
+// Computes the signal function according to chosen option
+//
+
+  switch (fSignalType)
+  {
+    case kBreitWigner:
+      return NormBreitWigner(x, par);
+    case kGaus:
+      return NormGaus(x, par);
+    default:
+      AliError("Invalid signal function");
+      return 0.0;
+  }
+}
+  
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::Background(Double_t *x, Double_t *par)
+{
+//
+// Computes the background function according to chosen option
+//
+
+  switch (fBackgroundType)
+  {
+    case kSquareRoot:
+      return NormSquareRoot(x, par);
+    case kLinear:
+      return NormLinear(x, par);
+    case kPoly2:
+      return NormPoly2(x, par);
+    default:
+      AliError("Invalid background function");
+      return 0.0;
+  }
+}
+
+//__________________________________________________________________________________________________
+Double_t AliRsnFitResult::Sum(Double_t *x, Double_t *par)
+{
+//
+// Computes the sum of the signal and the background chosen.
+// First 3 parameters are for the signal (they are always 3),
+// and other parameters, up to the necessary amount, are for
+// the background.
+//
+
+  Double_t signal     = Signal(x, par);
+  Double_t background = Background(x, &par[3]);
+
+  return signal + background;
+}
+
+//__________________________________________________________________________________________________
+Int_t AliRsnFitResult::GetNParBackground()
+{
+//
+// Tells how many parameters has the background function
+//
+
+  switch (fBackgroundType)
+  {
+    case kSquareRoot: return 2;
+    case kLinear    : return 2;
+    case kPoly2     : return 3;
+    default         : return 0;
+  }
+}
+
+//__________________________________________________________________________________________________
+Bool_t AliRsnFitResult::InitFunctions()
+{
+//
+// Initialize functions
+//
+
+  //if (!fSignalFcn) fSignalFcn = new TF1(Form("sg%s", GetName()), Signal, fFullRange[0], fFullRange[1], 3);
+  //if (!fBackgroundFcn) fBackgroundFcn = new TF1(Form("bg%s", GetName()), Background, fFullRange[0], fFullRange[1], GetNParBackground());
+  //if (!fSumFcn) fSumFcn = new TF1(Form("sum%s", GetName()), Sum, fFullRange[0], fFullRange[1], 3+GetNParBackground());
+}
+
+//__________________________________________________________________________________________________
+void AliRsnFitResult::SetHistogram(TH1F *histogram)
+{
+//
+// Clones the passed histogram to proceed with fits
+//
+
+  if (fHistogram) delete fHistogram;
+  
+  fHistogram = (TH1F*)histogram->Clone();
+  fHistogram->SetName(Form("h_%s_%s", GetName(), histogram->GetName()));
+  fHistogram->SetTitle(histogram->GetTitle());
+}
+
+//__________________________________________________________________________________________________
+Bool_t AliRsnFitResult::SingleFit(const char *opt, Double_t mass, Double_t width)
+{
+//
+// Executes a single fit of the function
+//
+/*
+  // require histogram
+  if (!fHistogram)
+  {
+    AliError("Require an initialized histogram!");
+    return kFALSE;
+  }
+  
+  // if necessary, initialize functions
+  if (!fSignalFcn || !fBackgroundFcn || !fSumFcn) InitFunctions();
+  
+  // step #0: know how many parameter we have
+  Int_t npar = GetNParBackground();
+  
+  // step #1:
+  // fit outside peak to roughly initialize the background
+  for (Int_t i = 0; i < npar; i++) fBackgroundFcn->SetParameter(i, 0.5);
+  status = (Int_t)h->Fit(fcnOut, opt);
+  if (status) return kFALSE;
+
+  // step #2:
+  // estimate signal/background in the peak
+  Int_t    imax = h->GetMaximumBin();
+  Double_t xmax = h->GetBinCenter(imax);
+  Double_t ymax = h->GetMaximum() - fBackgroundFcn->Eval(xmax);
+
+  // step #3:
+  // fit whole range with signal + background
+  // which is initialized to the reference values for mass and width
+  fSignalFcn->SetParameter(0, ymax);
+  fSignalFcn->SetParameter(1, mass);
+  fSignalFcn->SetParameter(2, width);
+  for (Int_t i = 0; i < npar; i++) fcnSum->SetParameter(i + 3, fBackgroundFcn->GetParameter(i));
+  status = (Int_t)h->Fit(fcnSum, opt);
+  if (status) return kFALSE;
+  
+  // get integral and its error here
+  fValue[kSumIntegralError] = fSumFcn->IntegralError(fitPeakRangeMin, fitPeakRangeMax);
+  fValue[kSumIntegral]      = fSumFcn->Integral     (fitPeakRangeMin, fitPeakRangeMax);
+  fValue[kChi2]             = fSumFcn->GetChisquare();
+  fValue[kNDF]              = fSumFcn->GetNDF();
+*/  
+  return kTRUE;
+}
diff --git a/PWG2/RESONANCES/AliRsnFitResult.h b/PWG2/RESONANCES/AliRsnFitResult.h
new file mode 100644 (file)
index 0000000..bcf7f87
--- /dev/null
@@ -0,0 +1,84 @@
+//
+// Class AliRsnFitResult
+//
+// It contains all results from a fit to an invariant mass spectrum
+//
+
+#ifndef ALIRSNFITRESULT
+#define ALIRSNFITRESULT
+
+#include <TMath.h>
+
+class AliRsnFitResult : public TNamed
+{
+  public:
+  
+    enum
+    {
+      kSgAmp,
+      kSgCenter,
+      kSgWidth,
+      kSgIntegral,
+      kBgIntegral,
+      kSg2Bg,
+      kSignificance,
+      kTotalIntegral,
+      kSumIntegral,
+      kSumIntegralError,
+      kChi2,
+      kNDF,
+      kValues
+    };
+    
+    enum ESignalType
+    {
+      kBreitWigner,
+      kGaus
+    };
+    
+    enum EBackgroundType
+    {
+      kSquareRoot,
+      kLinear,
+      kPoly2
+    };
+  
+    AliRsnFitResult(const char *name = "none", ESignalType sType = kGaus, EBackgroundType bType = kPoly2);
+    //AliRsnFitResult(const AliRsnFitResult& copy);
+    //AliRsnFitResult& operator=(cont AliRsnFitResult& copy);
+    
+    Double_t NormSquareRoot(Double_t *x, Double_t *par);
+    Double_t NormPoly2(Double_t *x, Double_t *par);
+    Double_t NormLinear(Double_t *x, Double_t *par);
+    Double_t NormBreitWigner(Double_t *x, Double_t *par);
+    Double_t NormGaus(Double_t *x, Double_t *par);
+    Double_t Signal(Double_t *x, Double_t *par);
+    Double_t Background(Double_t *x, Double_t *par);
+    Double_t Sum(Double_t *x, Double_t *par);
+    
+    Int_t    GetNParBackground();
+    Bool_t   InitFunctions();
+    Bool_t   SingleFit(const char *fitOpts, Double_t mass, Double_t width);
+    
+    void     SetHistogram(TH1F *histogram);
+    void     SetPeakRange(Double_t min, Double_t max) {fPeakRange[0] = TMath::Min(min, max); fPeakRange[1] = TMath::Max(min, max);}
+    void     SetFullRange(Double_t min, Double_t max) {fFullRange[0] = TMath::Min(min, max); fFullRange[1] = TMath::Max(min, max);}
+  
+  private:
+  
+    ESignalType      fSignalType;          // signal function
+    EBackgroundType  fBackgroundType;      // background function
+    
+    Double_t         fFullRange[2];        // full fit range
+    Double_t         fPeakRange[2];        // peak fit range
+    Double_t         fFitResult[kValues];  // array of fit results
+    
+    TH1F            *fHistogram;           // histogram to be fitted
+    TF1             *fSignalFcn;           // function for signal
+    TF1             *fBackgroundFcn;       // function for background
+    TF1             *fSumFcn;              // total function
+  
+    ClassDef(AliRsnFitResult, 1)
+};
+
+#endif
index a5d109c94f5c3569c44b3df64a3e0527d0534e90..05a13f1ae670ece9b89890c328036318ca8a1eaf 100644 (file)
@@ -138,7 +138,7 @@ TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTit
 // even if the bins are equal, since they are defined in this class.
 // Eventually present histoDef's in other slots of array (1, 2) are ignored.
 //
-// This version produces a THnSparseD.
+// This version produces a THnSparseF.
 //
 
   fSize = fAxisList.GetEntries();
@@ -176,13 +176,13 @@ TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTit
   switch (fSize)
   {
     case 1:
-      fH1 = new TH1D(histoName, histoTitle, nbins[0], min[0], max[0]);
+      fH1 = new TH1F(histoName, histoTitle, nbins[0], min[0], max[0]);
       break;
     case 2:
-      fH1 = new TH2D(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
+      fH1 = new TH2F(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
       break;
     case 3:
-      fH1 = new TH3D(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1], nbins[2], min[2], max[2]);
+      fH1 = new TH3F(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1], nbins[2], min[2], max[2]);
       break;
   }
   fH1->Sumw2();
@@ -191,7 +191,7 @@ TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTit
 }
 
 //________________________________________________________________________________________
-THnSparseD* AliRsnFunction::CreateHistogramSparse(const char *histoName, const char *histoTitle)
+THnSparseF* AliRsnFunction::CreateHistogramSparse(const char *histoName, const char *histoTitle)
 {
 //
 // Creates and returns the histogram defined using
@@ -200,7 +200,7 @@ THnSparseD* AliRsnFunction::CreateHistogramSparse(const char *histoName, const c
 // even if the bins are equal, since they are defined in this class.
 // Eventually present histoDef's in other slots of array (1, 2) are ignored.
 //
-// This version produces a THnSparseD.
+// This version produces a THnSparseF.
 //
 
   fSize = fAxisList.GetEntries();
@@ -236,7 +236,7 @@ THnSparseD* AliRsnFunction::CreateHistogramSparse(const char *histoName, const c
   }
 
   // create histogram
-  fHSparse = new THnSparseD(histoName, histoTitle, size, nbins, min, max);
+  fHSparse = new THnSparseF(histoName, histoTitle, size, nbins, min, max);
   fHSparse->Sumw2();
   
   // clean heap
@@ -284,19 +284,19 @@ Bool_t AliRsnFunction::Fill()
     {
       case 1:
         {
-          TH1D *h1 = (TH1D*)fH1;
+          TH1F *h1 = (TH1F*)fH1;
           h1->Fill(values[0]);
         }
         break;
       case 2:
         {
-          TH2D *h2 = (TH2D*)fH1;
+          TH2F *h2 = (TH2F*)fH1;
           h2->Fill(values[0], values[1]);
         }
         break;
       case 3:
         {
-          TH3D *h3 = (TH3D*)fH1;
+          TH3F *h3 = (TH3F*)fH1;
           h3->Fill(values[0], values[1], values[2]);
         }
         break;
@@ -309,7 +309,7 @@ Bool_t AliRsnFunction::Fill()
   {
     // check presence of output histogram
     if (!fHSparse) {
-      AliError("Required a THnSparseD whish is not initialized");
+      AliError("Required a THnSparseF which is not initialized");
       return kFALSE;
     }
     
index d2cf00c12e436f9af935d00058cc43b20a820bcb..fe3b34ed0688a7699cd0f9f474ed419be873ac9b 100644 (file)
@@ -60,7 +60,7 @@ class AliRsnFunction : public TNamed
     void                 AddAxis(AliRsnValue* const axis);
     Int_t                GetNumberOfAxes() {return fAxisList.GetEntries();}
     TH1*                 CreateHistogram(const char *histoName, const char *histoTitle);
-    THnSparseD*          CreateHistogramSparse(const char *histoName, const char *histoTitle);
+    THnSparseF*          CreateHistogramSparse(const char *histoName, const char *histoTitle);
 
     Bool_t               Fill();
 
@@ -75,7 +75,7 @@ class AliRsnFunction : public TNamed
     Bool_t              fUseTH1;      // use TH1 or not?
     Int_t               fSize;        // number of dim of output histogram
     TH1                *fH1;          // output histogram (standard type)
-    THnSparseD         *fHSparse;     // output histogram (sparse type)
+    THnSparseF         *fHSparse;     // output histogram (sparse type)
 
     // ROOT dictionary
     ClassDef(AliRsnFunction, 3)
index 3f7402533967ecd8b0e84a8fc576cedd6efcedf5..77e17455e8d10ad0aa206631f1feeb23f489a3ed 100644 (file)
@@ -18,6 +18,7 @@ ClassImp(AliRsnValue)
 
 //_____________________________________________________________________________
 AliRsnValue::AliRsnValue() :
+  TNamed(),
   fType(kValueTypes),
   fNBins(0),
   fMin(0.0),
@@ -32,7 +33,8 @@ AliRsnValue::AliRsnValue() :
 
 //_____________________________________________________________________________
 AliRsnValue::AliRsnValue
-(EValueType type, Int_t nbins, Double_t min, Double_t max) :
+(const char *name, EValueType type, Int_t nbins, Double_t min, Double_t max) :
+  TNamed(name, ""),
   fType(type),
   fNBins(0),
   fMin(0.0),
@@ -49,7 +51,8 @@ AliRsnValue::AliRsnValue
 
 //_____________________________________________________________________________
 AliRsnValue::AliRsnValue
-(EValueType type, Double_t min, Double_t max, Double_t step) :
+(const char *name, EValueType type, Double_t min, Double_t max, Double_t step) :
+  TNamed(name, ""),
   fType(type),
   fNBins(0),
   fMin(0.0),
@@ -63,6 +66,7 @@ AliRsnValue::AliRsnValue
   SetBins(min, max, step);
 }
 
+/*
 //_____________________________________________________________________________
 const char* AliRsnValue::GetName() const
 {
@@ -87,6 +91,7 @@ const char* AliRsnValue::GetName() const
     default:              return "UNDEF";
   }
 }
+*/
 
 //_____________________________________________________________________________
 TArrayD AliRsnValue::GetArray() const
index 7048177b551cfbf08dab0c318b745c1fa373ada3..f3c9b78e663912ce97082778d6cf7fb8b36c7165 100644 (file)
@@ -9,12 +9,13 @@
 #ifndef ALIRSNVALUE_H
 #define ALIRSNVALUE_H
 
+#include "TNamed.h"
 #include "TArrayD.h"
 
 class AliRsnPairDef;
 class AliRsnMother;
 
-class AliRsnValue : public TObject
+class AliRsnValue : public TNamed
 {
   public:
 
@@ -36,13 +37,13 @@ class AliRsnValue : public TObject
     };
 
     AliRsnValue();
-    AliRsnValue(EValueType type, Int_t n, Double_t min, Double_t max);
-    AliRsnValue(EValueType type, Double_t min, Double_t max, Double_t step);
-    AliRsnValue(const AliRsnValue& copy) : TObject(copy),fType(copy.fType),fNBins(copy.fNBins),fMin(copy.fMin),fMax(copy.fMax),fValue(copy.fValue) {}
-    AliRsnValue& operator=(const AliRsnValue& copy) {fType=copy.fType;fNBins=copy.fNBins;fMin=copy.fMin;fMax=copy.fMax;fValue=copy.fValue;return (*this);}
+    AliRsnValue(const char *name, EValueType type, Int_t n, Double_t min, Double_t max);
+    AliRsnValue(const char *name, EValueType type, Double_t min, Double_t max, Double_t step);
+    AliRsnValue(const AliRsnValue& copy) : TNamed(copy),fType(copy.fType),fNBins(copy.fNBins),fMin(copy.fMin),fMax(copy.fMax),fValue(copy.fValue) {}
+    AliRsnValue& operator=(const AliRsnValue& copy) {SetName(copy.GetName());fType=copy.fType;fNBins=copy.fNBins;fMin=copy.fMin;fMax=copy.fMax;fValue=copy.fValue;return (*this);}
     virtual ~AliRsnValue() { }
     
-    const char* GetName() const;
+    //const char* GetName() const;
     Int_t       GetNBins() const {return fNBins;}
     Double_t    GetMin() const {return fMin;}
     Double_t    GetMax() const {return fMax;}
index af9fb679ee82d2faee8e0b9da59503a1b57abfe3..96815b8e6feb2229f47fc4a4dfe27f4514eedc42 100644 (file)
@@ -11,7 +11,7 @@
 //                   arguments must have to be the task and the 'dataLabel' argument.
 // 
 Bool_t AddAnalysisTaskRsnTest
-(const char *dataLabel, const char *configMacro = "ConfigTaskRsn2010.C")
+(const char *dataLabel, const char *configMacro = "ConfigTaskRsnTest.C")
 {
   // retrieve analysis manager
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
@@ -19,6 +19,7 @@ Bool_t AddAnalysisTaskRsnTest
   // initialize task with all available slots, even if not all of them will be used:
   AliRsnAnalysisSE *task = new AliRsnAnalysisSE("RsnAnalysis");
   task->SetZeroEventPercentWarning(100.0);
+  task->SelectCollisionCandidates();
 
   // load and execute configuration macro
   gROOT->LoadMacro(configMacro);
index 9a61a63c089fbbaf44ac711483e3ec3b27462924..9dad7c9e40c7f962c2e27a4b12135c99304306f2 100644 (file)
@@ -30,13 +30,14 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   Bool_t isData  = strDataLabel.Contains("data");
   Bool_t isPass1 = strDataLabel.Contains("pass1");
   Bool_t isPass2 = strDataLabel.Contains("pass2");
+  Bool_t isPID   = strDataLabel.Contains("PID");
 
   //
   // -- Set cuts for events (applied to all analyses) -----------------------------------------------
   //
   
   // primary vertex range
-  AliRsnCutPrimaryVertex *cutVertex   = new AliRsnCutPrimaryVertex("cutVertex", 3);
+  AliRsnCutPrimaryVertex *cutVertex   = new AliRsnCutPrimaryVertex("cutVertex", 10.0, 0, kFALSE);
   AliRsnCutSet           *cutSetEvent = new AliRsnCutSet("eventCuts", AliRsnCut::kEvent);
   cutSetEvent->AddCut(cutVertex);
   cutSetEvent->SetCutScheme("cutVertex");
@@ -56,10 +57,10 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   //AliRsnPairNtuple      *truePM = new AliRsnPairNtuple("TruePM" , pairDefPM);
   //AliRsnPairNtuple      *pairPP = new AliRsnPairNtuple("PairPP" , pairDefPP);
   //AliRsnPairNtuple      *pairMM = new AliRsnPairNtuple("PairMM" , pairDefMM);
-  AliRsnPairFunctions    *pairPM = new AliRsnPairNtuple("PairPM" , pairDefPM);
-  AliRsnPairFunctions    *truePM = new AliRsnPairNtuple("TruePM" , pairDefPM);
-  AliRsnPairFunctions    *pairPP = new AliRsnPairNtuple("PairPP" , pairDefPP);
-  AliRsnPairFunctions    *pairMM = new AliRsnPairNtuple("PairMM" , pairDefMM);
+  AliRsnPairFunctions    *pairPM = new AliRsnPairFunctions(Form("PairPM%s", (isPID ? "pid" : "nopid")), pairDefPM);
+  AliRsnPairFunctions    *truePM = new AliRsnPairFunctions(Form("TruePM%s", (isPID ? "pid" : "nopid")), pairDefPM);
+  AliRsnPairFunctions    *pairPP = new AliRsnPairFunctions(Form("PairPP%s", (isPID ? "pid" : "nopid")), pairDefPP);
+  AliRsnPairFunctions    *pairMM = new AliRsnPairFunctions(Form("PairMM%s", (isPID ? "pid" : "nopid")), pairDefMM);
 
   //
   // -- Setup cuts ----------------------------------------------------------------------------------
@@ -67,14 +68,14 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   
   // track cut -----------------------
   // --> global cuts for 2010 analysis
-  AliRsnCutESD2010 *cuts2010 = new AliRsnCutESD2010("cuts2010");
+  AliRsnCutESD2010 *cuts2010 = new AliRsnCutESD2010(Form("cutESD2010%s", (isPID ? "pid" : "nopid")));
   // ----> set the flag for sim/data management
   cuts2010->SetMC(isSim);
-  // ----> require to check PID
-  cuts2010->SetCheckITS(1);
-  cuts2010->SetCheckTPC(1);
-  cuts2010->SetCheckTOF(1);
-  // ----> set TPC ranges and calibration
+  // ----> require to check PID or not, depending on the label
+  cuts2010->SetCheckITS(isPID);
+  cuts2010->SetCheckTPC(isPID);
+  cuts2010->SetCheckTOF(isPID);
+  // ----> set TPC rangeisPID and calibration
   cuts2010->SetTPCrange(5.0, 3.0);
   cuts2010->SetTPCpLimit(0.35);
   cuts2010->SetITSband(4.0);
@@ -126,7 +127,7 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   // --> only common cuts for tracks are needed
   AliRsnCutSet *cutSetDaughterCommon = new AliRsnCutSet("commonDaughterCuts", AliRsnCut::kDaughter);
   cutSetDaughterCommon->AddCut(cuts2010);
-  cutSetDaughterCommon->SetCutScheme("cuts2010");
+  cutSetDaughterCommon->SetCutScheme(cuts2010->GetName());
    
   // configure cut managers -------------------
   pairPM->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
@@ -142,9 +143,13 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   //
   
   // function axes
-  AliRsnValue *axisIM = new AliRsnValue(AliRsnValue::kPairInvMass,   4000,  0.9,  4.9);
-  AliRsnValue *axisPt = new AliRsnValue(AliRsnValue::kPairPt,         100,  0.0, 10.0);
-  AliRsnValue *axisY  = new AliRsnValue(AliRsnValue::kPairY,           30, -1.5,  1.5);
+  AliRsnValue *axisIM = new AliRsnValue("IM", AliRsnValue::kPairInvMass,   4000,  0.9,  4.9);
+  AliRsnValue *axisPt = new AliRsnValue("PT", AliRsnValue::kPairPt,         100,  0.0, 10.0);
+  AliRsnValue *axisY[4];
+  axisY[0] = new AliRsnValue("Y1", AliRsnValue::kPairY, 1, -0.5,  0.5);
+  axisY[1] = new AliRsnValue("Y2", AliRsnValue::kPairY, 1, -0.6,  0.6);
+  axisY[2] = new AliRsnValue("Y3", AliRsnValue::kPairY, 1, -0.7,  0.7);
+  axisY[3] = new AliRsnValue("Y4", AliRsnValue::kPairY, 1, -0.8,  0.8);
   
   // the ntuple output requires to get directly the values
   //pairPM->AddValue(axisIM);
@@ -161,17 +166,22 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   //pairMM->AddValue(axisEta);
   
   // functions
-  AliRsnFunction *fcnImPtY = new AliRsnFunction;
-  // --> add axes
-  fcnImPtY->AddAxis(axisIM);
-  fcnImPtY->AddAxis(axisPt);
-  fcnImPtY->AddAxis(axisEta);
+  AliRsnFunction *fcnImPtY[4];
+  for (Int_t i = 0; i < 4; i++)
+  {
+    fcnImPtY[i] = new AliRsnFunction;
+    //fcnImPtY->UseSparse();
+    // --> add axes
+    fcnImPtY[i]->AddAxis(axisIM);
+    fcnImPtY[i]->AddAxis(axisPt);
+    fcnImPtY[i]->AddAxis(axisY[i]);
   
-  // add functions to pairs
-  pairPM->AddFunction(fcnImPtY);
-  truePM->AddFunction(fcnImPtY);
-  pairPP->AddFunction(fcnImPtY);
-  pairMM->AddFunction(fcnImPtY);
+    // add functions to pairs
+    pairPM->AddFunction(fcnImPtY[i]);
+    truePM->AddFunction(fcnImPtY[i]);
+    pairPP->AddFunction(fcnImPtY[i]);
+    pairMM->AddFunction(fcnImPtY[i]);
+  }
   
   //
   // -- Conclusion ----------------------------------------------------------------------------------
index 450e0d166ace6b2c23c0f60b2f9a893c18c95692..231c3000375142edeba2756f3ebe351e10785a33 100644 (file)
@@ -24,19 +24,19 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   
   // interpret the useful information from second argument
   TString strDataLabel(dataLabel);
-  Bool_t isESD   = strDataLabel->Contains("ESD");
-  Bool_t isAOD   = strDataLabel->Contains("AOD");
-  Bool_t isSim   = strDataLabel->Contains("sim");
-  Bool_t isData  = strDataLabel->Contains("data");
-  Bool_t isPass1 = strDataLabel->Contains("pass1");
-  Bool_t isPass2 = strDataLabel->Contains("pass2");
+  Bool_t isESD   = strDataLabel.Contains("ESD");
+  Bool_t isAOD   = strDataLabel.Contains("AOD");
+  Bool_t isSim   = strDataLabel.Contains("sim");
+  Bool_t isData  = strDataLabel.Contains("data");
+  Bool_t isPass1 = strDataLabel.Contains("pass1");
+  Bool_t isPass2 = strDataLabel.Contains("pass2");
 
   //
   // -- Set cuts for events (applied to all analyses) -----------------------------------------------
   //
   
   // primary vertex range
-  AliRsnCutPrimaryVertex *cutVertex   = new AliRsnCutPrimaryVertex("cutVertex", 3);
+  AliRsnCutPrimaryVertex *cutVertex   = new AliRsnCutPrimaryVertex("cutVertex", 10.0, 0, kFALSE);
   AliRsnCutSet           *cutSetEvent = new AliRsnCutSet("eventCuts", AliRsnCut::kEvent);
   cutSetEvent->AddCut(cutVertex);
   cutSetEvent->SetCutScheme("cutVertex");
@@ -47,11 +47,17 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   //
 
   // decay channels
-  AliRsnPairDef         *pairDef = new AliRsnPairDef(AliPID::kKaon, '+', AliPID::kKaon, '-', 333, 1.019455);
+  AliRsnPairDef         *pairDefpp = new AliRsnPairDef(AliPID::kKaon, '+', AliPID::kKaon, '+', 333, 1.019455);
+  AliRsnPairDef         *pairDefpm = new AliRsnPairDef(AliPID::kKaon, '+', AliPID::kKaon, '-', 333, 1.019455);
+  AliRsnPairDef         *pairDefmm = new AliRsnPairDef(AliPID::kKaon, '-', AliPID::kKaon, '-', 333, 1.019455);
 
   // computation objects
-  AliRsnPairFunctions   *pairPMhist = new AliRsnPairFunctions("pairHist", pairDef);
-  AliRsnPairNtuple      *pairPMntp  = new AliRsnPairNtuple   ("pairNtp" , pairDef);
+  AliRsnPairFunctions   *pairPMhist = new AliRsnPairFunctions("pairPMHist", pairDefpm);
+  AliRsnPairFunctions   *pairPPhist = new AliRsnPairFunctions("pairPPHist", pairDefpp);
+  AliRsnPairFunctions   *pairMMhist = new AliRsnPairFunctions("pairMMHist", pairDefmm);
+  AliRsnPairNtuple      *pairPMntp  = new AliRsnPairNtuple   ("pairPMNtp" , pairDefpm);
+  AliRsnPairNtuple      *pairPPntp  = new AliRsnPairNtuple   ("pairPPNtp" , pairDefpp);
+  AliRsnPairNtuple      *pairMMntp  = new AliRsnPairNtuple   ("pairMMNtp" , pairDefmm);
 
   //
   // -- Setup cuts ----------------------------------------------------------------------------------
@@ -63,17 +69,17 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   // ----> set the flag for sim/data management
   cuts2010->SetMC(isSim);
   // ----> require to check PID
-  cuts2010->SetCheckITS(1);
-  cuts2010->SetCheckTPC(1);
-  cuts2010->SetCheckTOF(1);
+  cuts2010->SetCheckITS(kTRUE);
+  cuts2010->SetCheckTPC(kTRUE);
+  cuts2010->SetCheckTOF(kTRUE);
   // ----> set TPC ranges and calibration
   cuts2010->SetTPCrange(5.0, 3.0);
   cuts2010->SetTPCpLimit(0.35);
   cuts2010->SetITSband(4.0);
-  if (isMC) cuts2010->SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
-  else      cuts2010->SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
+  if (isSim) cuts2010->SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
+  else       cuts2010->SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
   // ----> set standard quality cuts for TPC global tracks
-  cuts2010->GetCutsTPC()->SetRequireTPCStandAlone(kTRUE); // require to have the projection at inner TPC wall
+  //cuts2010->GetCutsTPC()->SetRequireTPCStandAlone(kTRUE); // require to have the projection at inner TPC wall
   cuts2010->GetCutsTPC()->SetMinNClustersTPC(70);
   cuts2010->GetCutsTPC()->SetMaxChi2PerClusterTPC(4);
   cuts2010->GetCutsTPC()->SetAcceptKinkDaughters(kFALSE);
@@ -85,8 +91,7 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   cuts2010->GetCutsTPC()->SetDCAToVertex2D(kFALSE); // each DCA is checked separately
   cuts2010->GetCutsTPC()->SetRequireSigmaToVertex(kFALSE);
   // ----> set standard quality cuts for ITS standalone tracks
-  cuts2010->GetCutsITS()->SetRequireITSStandAlone(kTRUE);
-  cuts2010->GetCutsITS()->SetRequireITSPureStandAlone(kFALSE);
+  cuts2010->GetCutsITS()->SetRequireITSStandAlone(kTRUE, kTRUE);
   cuts2010->GetCutsITS()->SetRequireITSRefit(kTRUE);
   cuts2010->GetCutsITS()->SetMinNClustersITS(4);
   cuts2010->GetCutsITS()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
@@ -94,7 +99,6 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   cuts2010->GetCutsITS()->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55"); // DCA pt dependent
   cuts2010->GetCutsITS()->SetMaxDCAToVertexZ(1e6); // disabled
   cuts2010->GetCutsITS()->SetDCAToVertex2D(kFALSE); // each DCA is checked separately
-  cuts2010->GetCutsITS()->SetRequireITSPid(kTRUE);
   // ----> set the configuration for TOF PID checks
   if (isData && (isPass1 || isPass2))
   {
@@ -175,38 +179,42 @@ Bool_t RsnConfigTask(AliRsnAnalysisSE* &task, const char *dataLabel)
   // cut managers
   // define a proper name for each mult bin, to avoid omonyme output histos
   pairPMhist->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
+  pairPPhist->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
+  pairMMhist->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
   //pairPMhist->GetCutManager()->SetMotherCuts        (cutSetMother);
   pairPMntp ->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
+  pairPPntp ->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
+  pairMMntp ->GetCutManager()->SetCommonDaughterCuts(cutSetDaughterCommon);
   //pairPMntp ->GetCutManager()->SetMotherCuts        (cutSetMother);
 
   // function axes
-  AliRsnValue *axisIM   = new AliRsnValue(AliRsnValue::kPairInvMass,   50,  0.9,  1.4);
-  AliRsnValue *axisPt   = new AliRsnValue(AliRsnValue::kPairPt,        50,  0.0, 10.0);
-  AliRsnValue *axisEta  = new AliRsnValue(AliRsnValue::kPairEta,       30, -1.5,  1.5);
+  AliRsnValue *axisIM   = new AliRsnValue("IM", AliRsnValue::kPairInvMass,   50,  0.9,  1.4);
+  AliRsnValue *axisPt   = new AliRsnValue("PT", AliRsnValue::kPairPt,        50,  0.0, 20.0);
   
   pairPMntp->AddValue(axisIM);
   pairPMntp->AddValue(axisPt);
-  pairPMntp->AddValue(axisEta);
+  pairPPntp->AddValue(axisIM);
+  pairPPntp->AddValue(axisPt);
+  pairMMntp->AddValue(axisIM);
+  pairMMntp->AddValue(axisPt);
 
   // functions
-  AliRsnFunction *fcn      = new AliRsnFunction;
   AliRsnFunction *fcnPt    = new AliRsnFunction;
-  AliRsnFunction *fcnPtEta = new AliRsnFunction;
   // --> add axes
-  fcn     ->AddAxis(axisIM);
   fcnPt   ->AddAxis(axisIM);
   fcnPt   ->AddAxis(axisPt);
-  fcnPtEta->AddAxis(axisIM);
-  fcnPtEta->AddAxis(axisPt);
-  fcnPtEta->AddAxis(axisEta);
   
-  pairPMhist->AddFunction(fcn);
   pairPMhist->AddFunction(fcnPt);
-  pairPMhist->AddFunction(fcnPtEta);
+  pairPPhist->AddFunction(fcnPt);
+  pairMMhist->AddFunction(fcnPt);
   
   // add everything to pair manager
   task->GetAnalysisManager()->Add(pairPMhist);
+  task->GetAnalysisManager()->Add(pairPPhist);
+  task->GetAnalysisManager()->Add(pairMMhist);
   task->GetAnalysisManager()->Add(pairPMntp);
+  task->GetAnalysisManager()->Add(pairPPntp);
+  task->GetAnalysisManager()->Add(pairMMntp);
 
   return kTRUE;
 }
diff --git a/PWG2/RESONANCES/macros/test/PluginDataByRun.C b/PWG2/RESONANCES/macros/test/PluginDataByRun.C
new file mode 100644 (file)
index 0000000..61b2623
--- /dev/null
@@ -0,0 +1,177 @@
+// ======================================
+// ===== ALIEN PLUGIN CONFIGURATION =====
+// ======================================
+//
+// This example macro configures an AliAnalysisAlien object
+// to be used to run analysis on AliEn without having to put
+// all the required code and scripts in the alien shell.
+//
+// All the possible configuration parameters are arguments
+// of the macro function, even if most of them have default
+// values which the user will rarely change, and some args
+// do exclude some other ones, depending of user choice.
+//
+// The macro tries to synchronize some output names, using 
+// a unique name ('analysisName') to define all files that
+// describe the output, the analysis macros/executables/JDL.
+//
+// Since the run mode can be more variable than the config
+// it is not set here, but it is required in the run macro
+// which uses the plugin.
+//
+// Considered that the arguments are many, they are explained
+// inside the list of arguments in the macro definition.
+// In ALL cases where a list of strings must be provided, its
+// elements must be separated by a blank character.
+//
+AliAnalysisAlien* PluginDataByRun
+(
+  // common base name for defining those of all created files (run numbers will be appended)
+  const char *analysisName = "AnalysisPhi7TeV",
+
+  // input files to be processed:
+  // -- runList   --> list of run numbers, separated by a minus (e.g.: "23211-22231-2232") --> will name the output directory
+  // -- runPath   --> AliEn path where runs will be searched
+  // -- runPrefix --> string prefix prepended to run numbers
+  // -- pattern   --> name of the file to be used (ESD/AOD)
+  const char *runList     = "117112",
+  const char *runPath     = "/alice/data/2010/LHC10b",
+  const char *runPrefix   = "000",
+  const char *pattern     = "pass2/*/AliESDs.root",
+
+  // working parameters:
+  // -- workDir        --> AliEn working directory (w.r. to AliEn $HOME)
+  // -- outDir         --> AliEn output directory (w.t. to workDir)
+  // -- analysisName   --> common name used for analysis related files (.C, .JDL, exe, log)
+  // -- taskSource     --> name of a task compiled on-fly (without any extension)
+  // -- exePre         --> executable command part before task macro name
+  // -- exePost        --> executable command part after task macro name
+  // -- exeArgs        --> executable arguments to be specified in the JDL
+  // -- outList        --> lit of expected output files (separated by commas)
+  const char *workDir        = "rsnNewTest_v1",
+  const char *outDir         = "run",
+  const char *taskSource     = "",
+  const char *exePre         = "aliroot -q -b",
+  const char *exePost        = "",
+  const char *exeArgs        = "",
+  
+  // additional stuff (libs, includes, code):
+  // -- addInclude  --> additional include paths
+  // -- addLibs     --> additional libraries compiled in ROOT or ALIROOT
+  // -- addPars     --> additional PAR libraries from ALIROOT
+  // -- addExternal --> external tasks to be compiled on-fly
+  const char *addInclude  = "TOF",
+  const char *addLibs     = "",
+  const char *addPars     = "PWG2resonances.par",
+  const char *addExternal = "",
+  
+  // job parameters:
+  // these are all the job parameters which can be inserted in the JDL
+  // their names are the same as the corresponding JDL keywords
+  Int_t       split          = 50,
+  Int_t       maxInitFailed  = 0,
+  Int_t       resubmitThr    = 0,
+  Int_t       TTL            = 22500,
+  const char *inputFormat    = "xml-single",
+  const char *splitMode      = "se",
+  Int_t       price          = 1,
+  const char *jobTag         = "",
+  Int_t       maxMergeFiles  = 50,
+  Int_t       nTestFiles     = 6,
+  
+  // standard package versions
+  const char *apiVersion     = "V1.1x",
+  const char *rootVersion    = "v5-26-00b-6",
+  const char *aliVersion     = "v4-19-22-AN"
+)
+{
+  // append the run number to the reference names used 
+  // for the executeble and for naming the used files
+  Char_t exeName[255];
+  sprintf(exeName, "%s_exe", analysisName);
+  
+  // create plugin object
+  AliAnalysisAlien *plugin = new AliAnalysisAlien;
+  plugin->SetOutputToRunNo(kTRUE);
+  plugin->SetMergeViaJDL();
+  plugin->SetNtestFiles(nTestFiles);
+  plugin->SetNrunsPerMaster(0);
+  
+  // define names of analysis files after 'analysisName'
+  const char *analysisMacro = Form("%s.C"   , analysisName);
+  const char *logFile       = Form("%s.log" , analysisName);
+  const char *jdlFile       = Form("%s.jdl" , analysisName);
+  const char *outFile       = Form("%s.root", analysisName);
+  const char *exeFile       = Form("%s.sh"  , exeName);
+  const char *archiveFile   = Form("%s.zip:stdout,stderr,%s,%s@disk=2", analysisName, logFile, outFile);
+  
+  // package versions
+  plugin->SetAPIVersion(apiVersion);
+  plugin->SetROOTVersion(rootVersion);
+  plugin->SetAliROOTVersion(aliVersion);
+  
+  // work and output directory in AliEn shell
+  plugin->SetGridWorkingDir(workDir);
+  plugin->SetGridOutputDir(outDir);
+  
+  // executable command and arguments
+  // in 'post' args it is added the request to log all output
+  // into a file which is named after the common analysis name
+  plugin->SetExecutableCommand(exePre);
+  plugin->SetExecutableArgs(Form("%s >& %s", exePost, logFile));
+  plugin->SetExecutable(exeFile);
+  if (strlen(exeArgs) > 0) plugin->SetArguments(exeArgs);
+  
+  // declare (if any) the analysis source file
+  // to be compiled runtime separated by blanks
+  // and, in this case, add its header and implementation to additional libs
+  if (strlen(taskSource) > 0)
+  {
+    plugin->SetAnalysisSource(Form("%s.h", taskSource));
+    plugin->SetAdditionalLibs(Form("%s %s.h %s.cxx", addLibs, taskSource, taskSource));
+  }
+  else
+    plugin->SetAdditionalLibs(addLibs);
+  
+  // additional PARs
+  if (strlen(addPars) > 0) plugin->EnablePackage(addPars);
+  
+  // external packages
+  if (strlen(addExternal) > 0) plugin->AddExternalPackage(addExternal);
+    
+  // set names of analysis macro, JDL script, output file and archive
+  // from the commonly used 'outName' argument
+  plugin->SetAnalysisMacro(analysisMacro);
+  plugin->SetJDLName(jdlFile);
+  plugin->SetDefaultOutputs(kFALSE);
+  plugin->SetOutputFiles(outFile);
+  plugin->SetOutputArchive(archiveFile);
+  
+  // add all runs
+  TString sList = runList;
+  list  = sList.Tokenize("-");
+  Int_t i, n = list->GetEntries();
+  plugin->SetRunPrefix(runPrefix);
+  plugin->SetGridDataDir(runPath);
+  plugin->SetDataPattern(pattern);
+  for (i = 0; i < n; i++)
+  {
+    TObjString *os = (TObjString*)list->At(i);
+    plugin->AddRunNumber(os->GetString().Atoi());
+    cout << "Added run number " << os->GetString().Data() << endl;
+  }
+  
+  // JDL parameters
+  plugin->SetSplitMaxInputFileNumber(split);
+  plugin->SetMaxInitFailed(maxInitFailed);
+  plugin->SetMasterResubmitThreshold(resubmitThr);
+  plugin->SetTTL(TTL);
+  plugin->SetPrice(price);
+  plugin->SetInputFormat(inputFormat);
+  plugin->SetSplitMode(splitMode);
+  if (strlen(jobTag) > 0) plugin->SetJobTag(jobTag);
+  plugin->SetMaxMergeFiles(maxMergeFiles);
+    
+  // the end!
+  return plugin;
+}
diff --git a/PWG2/RESONANCES/macros/test/runAlienPlugin.C b/PWG2/RESONANCES/macros/test/runAlienPlugin.C
new file mode 100644 (file)
index 0000000..35aa062
--- /dev/null
@@ -0,0 +1,191 @@
+//
+// This is an example steering macro for running RSN analysis task
+// with the AliEn plugin to launch a multiple analysis.
+//
+// Inputs:
+//   - runList     = list of runs to be processed
+//   - runPath     = path containing the runs
+//   - runMode     = AliEn plugin run mode
+//   - pluginMacro = macro which loads and initializes the plugin
+//   - addTaskName = name of the macro to add the RSN analysis task
+//                   (assumed to have inside it a function named like the file)
+//   - inputSource = name of the file containing all the inputs
+//                   ---> to run on a local collection, the collection file 
+//                        must contain on each line the full path 
+//                        of one input file and it must have the ".txt" extension
+//                   ---> to run on an AliEn collection, the collection file must be an XML
+//                        file collection like those built from the "find -x" method in aliensh.
+//   - dataLabel   = a label which is used to know what kind of data are being read
+//                   (it is propagated to the 'addTask' macro for eventual setting up of something
+//   - outName     = name for the file with RSN package outputs (without ROOT extension)
+//                   in this case it is fundamental to define the names of all plugin objects
+//
+// Notes:
+//   - in case the source is an ESD, and if inputs are a MC production
+//     the MC input handler is created by default
+// 
+//
+// In principle, the user should never modify this macro. 
+//
+void runAlienPlugin
+(
+  const char *runList     = "117112-117116-117099-117220-117048-117109-117060-117054-117065",
+  const char *runPath     = "/alice/data/2010/LHC10b",
+  const char *runMode     = "terminate",
+  const char *pluginMacro = "PluginDataByRun.C",
+  const char *addTaskName = "AddAnalysisTaskRsnTest.C",
+  const char *dataLabel   = "7TeV_pass2_data_ESD",
+  const char *outName     = "rsnTest"
+)
+{
+  
+  // convert the last argument into a BOOL variable
+  TString strDataLabel(dataLabel);
+  Bool_t isESD = strDataLabel.Contains("ESD");
+  Bool_t isAOD = strDataLabel.Contains("AOD");
+  Bool_t isSim = strDataLabel.Contains("sim");   
+  
+  //AliLog::SetGlobalDebugLevel(AliLog::kDebug+2);
+
+  // load compiled libraries (for aliroot session)
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWG2resonances.so");
+  
+  //
+  // === PLUGIN CONFIGURATION =====================================================================
+  //
+  
+  // check token
+  if (!AliAnalysisGrid::CreateToken()) return;
+  
+  // load and execute plugin configuration macro
+  // pass to the macro, as FIRST argument, the common name
+  // which is used for the output, since it must be equal
+  // to the one defined here for the common output (for merging)
+  TString splugin(pluginMacro);
+  gROOT->LoadMacro(pluginMacro);
+  splugin.ReplaceAll(".C", Form("(\"%s\",\"%s\",\"%s\")", outName, runList, runPath));
+  AliAnalysisAlien *plugin = (AliAnalysisAlien*)gROOT->ProcessLine(splugin);
+  
+  // set run mode
+  plugin->SetRunMode(runMode);
+  
+  //
+  // === ANALYSIS MANAGER CONFIGURATION ===========================================================
+  //
+
+  // create analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("taskRsnTest");
+  mgr->SetGridHandler(plugin);
+  mgr->SetCommonFileName(Form("%s.root", outName));
+  
+  // create input handler
+  if (isESD)
+  {
+    AliESDInputHandler *esdHandler = new AliESDInputHandler();
+    mgr->SetInputEventHandler(esdHandler);
+    // if possible, create also MC handler
+    if (isSim)
+    {
+      AliMCEventHandler *mcHandler  = new AliMCEventHandler();
+      mgr->SetMCtruthEventHandler(mcHandler);
+    }
+  }
+  else if (isAOD)
+  {
+    AliAODInputHandler *aodHandler = new AliAODInputHandler();
+    mgr->SetInputEventHandler(aodHandler);
+  }
+  else
+  {
+    ::Error("Required an ESD or AOD input data set");
+    return;
+  }
+  
+  //
+  // === ANALYSIS TASK CREATION AND INCLUSION =====================================================
+  //
+  
+  // add event selection for data
+  gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isSim);
+  
+  // add task macro
+  gROOT->ProcessLine(Form(".x %s(\"%s\")", addTaskName, dataLabel));
+
+  // initialize and start analysis
+  if (!mgr->InitAnalysis()) return;
+  mgr->PrintStatus();
+  mgr->StartAnalysis("grid");
+}
+
+//_________________________________________________________________________________________________
+Bool_t LoadPars(const char *parList, const char *path)
+{
+//
+// Load PAR libraries locally
+// ---
+// Arguments:
+//  - parList = list of PARs without extension, separated by ':'
+//  - path    = path where PARs are stored
+//
+
+  // store position of working directory
+  TString ocwd = gSystem->WorkingDirectory();
+
+  // tokenize list
+  TString     pars(parList);
+  TObjArray  *array = pars.Tokenize(":");
+
+  // loop on list
+  TObjString *ostr;
+  TString     str;
+  Char_t      parName[200], parFile[200];
+  for (Int_t i = 0; i < array->GetEntriesFast(); i++)
+  {
+    ostr = (TObjString*) array->At(i);
+    str = ostr->GetString();
+    sprintf(parName, "%s", str.Data());
+    sprintf(parFile, "%s/%s.par", path, str.Data());
+
+    // check that file exists
+    if (!gSystem->AccessPathName(parFile))
+    {
+      // explode tar-ball and enter it
+      gROOT->ProcessLine(Form(".! tar xzf %s", parFile));
+      gSystem->ChangeDirectory(Form("%s", parName));
+      // checks for BUILD.sh and execute it
+      if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh"))
+      {
+        ::Info("", ">> Building PARs: %s", parName);
+        if (gSystem->Exec("PROOF-INF/BUILD.sh"))
+        {
+          ::Error("LoadPars", Form("BUILD.sh error for '%s'", parFile));
+          gSystem->ChangeDirectory(ocwd);
+          return kFALSE;
+        }
+      }
+      // check and execute SETUP.C
+      if (!gSystem->AccessPathName("PROOF-INF/SETUP.C"))
+      {
+        ::Info("", ">> Setting up PARs: %s", parName);
+        if (gROOT->Macro("PROOF-INF/SETUP.C"))
+        {
+          Error("LoadPars", Form("SETUP.C error for '%s'", parFile));
+          gSystem->ChangeDirectory(ocwd);
+          return kFALSE;
+        }
+      }
+    }
+    else
+    {
+      Error("LoadParsLocal", Form("File '%s' not found", parFile));
+      return kFALSE;
+    }
+  }
+
+  gSystem->ChangeDirectory(ocwd);
+  return kTRUE;
+}
index 9f14e899ee459b843f34016d0fd0061c93104cd9..aacd5637bbf4e7e97a636197bc7a5d62679db2b4 100644 (file)
 void runLocal
 (
   Int_t       nReadFiles  = 1,
-  Int_t       nSkipFiles  = 0,
+  Int_t       nSkipFiles  = 1,
   const char *addTaskName = "AddAnalysisTaskRsnTest.C",
-  const char *inputSource = "test.xml",
-  const char *dataLabel   = "7TeV_pass2_sim_ESD",
+  //const char *inputSource = "000117065.xml",
+  const char *inputSource = "local.txt",
+  const char *dataLabel   = "7TeV_pass2_data_ESD",
   const char *outName     = "rsnTest"
 )
 {
@@ -41,7 +42,7 @@ void runLocal
   Bool_t isAOD = strDataLabel.Contains("AOD");
   Bool_t isSim = strDataLabel.Contains("sim");   
   
-  //AliLog::SetGlobalDebugLevel(AliLog::kDebug+2);
+  //AliLog::SetGlobalDebugLevel(AliLog::kDebug+1);
 
   // check extension of input to distinguish between XML and TXT
   TString sInput(inputSource);
@@ -89,6 +90,10 @@ void runLocal
   gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPhysicsSelection.C");
   AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isSim);
   
+  // add event selection for data
+  gROOT->LoadMacro("AddTaskAnalysisPhi7TeV.C");
+  AddTaskAnalysisPhi7TeV(dataLabel);
+  
   // add task macro
   gROOT->ProcessLine(Form(".x %s(\"%s\")", addTaskName, dataLabel));
 
index 73794aef0109cc78429761e8c323d9776ab372db..1ed1996ef7aff10e4cb825974b7a07eee67be1a6 100644 (file)
@@ -33,6 +33,7 @@ SRCS= RESONANCES/AliRsnDaughter.cxx \
       RESONANCES/AliRsnAnalysisPhi7TeVNoPID.cxx \
       RESONANCES/AliRsnMonitorTrack.cxx \
       RESONANCES/AliRsnAnalysisMonitorTask.cxx \
+      RESONANCES/AliRsnFitResult.cxx \
 #     RESONANCES/AliRsnAnalysisEffSE.cxx \
 
 HDRS= $(SRCS:.cxx=.h)