Analysis task to study the extra SPD clusters.
authorcvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Apr 2012 10:28:07 +0000 (10:28 +0000)
committercvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Apr 2012 10:28:07 +0000 (10:28 +0000)
PWGPP/forward/SPDClustTask/AliSPDClustTask.cxx [new file with mode: 0755]
PWGPP/forward/SPDClustTask/AliSPDClustTask.h [new file with mode: 0755]
PWGPP/forward/SPDClustTask/AnalysisSPDClustTask.C [new file with mode: 0644]
PWGPP/forward/SPDClustTask/FitOutput.root [new file with mode: 0644]
PWGPP/forward/SPDClustTask/runAAFSPDTask.C [new file with mode: 0644]
PWGPP/forward/SPDClustTask/spectraCombine.root [new file with mode: 0644]

diff --git a/PWGPP/forward/SPDClustTask/AliSPDClustTask.cxx b/PWGPP/forward/SPDClustTask/AliSPDClustTask.cxx
new file mode 100755 (executable)
index 0000000..539c108
--- /dev/null
@@ -0,0 +1,487 @@
+/*************************************************************************
+* Copyright(c) 1998-2008, 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.                  * 
+**************************************************************************/
+
+///////////////////////////////////////////////////////////////////////////
+// Class AliSPDClustTask                                            //
+// Analysis task to produce data and MC histos needed for tracklets      //
+// dNdEta extraction in multiple bins in one go                          //
+// Author:  ruben.shahoyan@cern.ch                                       //
+///////////////////////////////////////////////////////////////////////////
+/*
+  Important parameters to set:
+  1) make sure to initialize correct geometry in UserCreateOutputObjects
+  2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand
+...
+*/
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TRandom.h"
+#include "TH1F.h"
+#include "TH2F.h" 
+#include "TList.h"
+#include "TObjArray.h"
+#include "TGeoGlobalMagField.h"
+
+#include "AliAnalysisManager.h"
+
+#include "AliMultiplicity.h"
+#include "AliESDEvent.h"  
+#include "AliESDInputHandler.h"
+#include "AliESDInputHandlerRP.h"
+#include "AliCDBPath.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliCDBStorage.h"
+#include "AliGeomManager.h"
+#include "AliMagF.h"
+#include "AliESDVZERO.h"
+#include "AliESDZDC.h"
+#include "AliRunLoader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "AliStack.h"
+#include "AliGenEventHeader.h"
+#include "AliCentrality.h"
+#include "../ITS/AliITSRecPoint.h"
+#include "../ITS/AliITSgeomTGeo.h"
+#include "../ITS/AliITSMultReconstructor.h" 
+
+#include "AliLog.h"
+
+#include "AliPhysicsSelection.h"
+#include "AliSPDClustTask.h"
+#include "AliITSMultReconstructor.h"
+#include "AliGenEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliESDtrackCuts.h"
+
+ClassImp(AliSPDClustTask)
+
+
+//________________________________________________________________________
+/*//Default constructor
+AliSPDClustTask::AliSPDClustTask(const char *name)
+  : AliAnalysisTaskSE(name),
+*/  
+//________________________________________________________________________
+AliSPDClustTask::AliSPDClustTask(const char *name) 
+  : AliAnalysisTaskSE(name), 
+//
+  fOutput(0), 
+  fUseMC(kFALSE),
+  fHistos(0),
+//
+  fEtaMin(-3.0),
+  fEtaMax(3.0),
+  fZVertexMin(-20),
+  fZVertexMax( 20),
+//
+  fScaleDTBySin2T(kFALSE),
+  fCutOnDThetaX(kFALSE),
+  fNStdDev(1.),
+  fDPhiWindow(0.08),
+  fDThetaWindow(0.025),
+  fDPhiShift(0.0045),
+  fPhiOverlapCut(0.005),
+  fZetaOverlap(0.05),
+  fRemoveOverlaps(kFALSE),
+//
+  fDPhiSCut(0.06),
+  fNStdCut(1.),
+//
+  fMultReco(0),
+  fRPTree(0),
+  fStack(0),
+  fMCEvent(0),
+  //
+  fDontMerge(kFALSE),
+  fhPtIn(0)
+{
+  // Constructor
+
+  DefineOutput(1, TList::Class());
+  //
+  SetScaleDThetaBySin2T();
+  SetNStdDev();
+  SetPhiWindow();
+  SetThetaWindow();
+  SetPhiShift();
+  SetPhiOverlapCut();
+  SetZetaOverlapCut();
+  SetRemoveOverlaps();
+  //
+}
+
+//________________________________________________________________________
+AliSPDClustTask::~AliSPDClustTask()
+{
+  // Destructor
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
+    printf("Deleteing output\n");
+    delete fOutput;
+    fOutput = 0;
+  }
+  //
+  delete fMultReco;
+  delete fHistos;
+  if (fhPtIn) delete fhPtIn;
+  //
+}
+
+//________________________________________________________________________
+void AliSPDClustTask::UserCreateOutputObjects() 
+{
+  //
+  //  OpenFile(1);  fDontMerge = kTRUE;
+  //  
+
+  fOutput = new TList();
+  fOutput->SetOwner(); 
+  //
+  AliCDBManager *man = AliCDBManager::Instance();
+  if (fUseMC) {
+    Bool_t newGeom = kTRUE;
+    man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
+    if (newGeom) {
+      // new geom
+      AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,8);
+      AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
+      if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
+    }
+    else {
+      // old geom
+      AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130845,7);
+      AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
+      if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
+    }
+  }
+  else {
+    man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
+    AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",137366);
+    AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
+    if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
+  }
+  //
+  // Create histograms
+  fHistos = BookHistos();
+  for(Int_t i=0; i<fHistos->GetEntriesFast(); i++) {
+    if(fHistos->At(i)){
+      fOutput->AddLast(fHistos->At(i));
+    }
+  }
+
+  PostData(1, fOutput);
+  //
+}
+
+//________________________________________________________________________
+void AliSPDClustTask::UserExec(Option_t *) 
+{
+  // Main loop
+  //
+  //
+  AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
+  fRPTree = 0;
+  AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler(); 
+  if (!handRP) { AliError("No RP handler"); return; }
+  //
+  fRPTree = handRP->GetTreeR("ITS");
+  if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
+  //
+  AliESDEvent *esd  = handRP->GetEvent();
+  if (!esd) { AliError("No AliESDEvent"); return; }
+  //
+  AliMCEventHandler* eventHandler = 0;
+  fMCEvent = 0;
+  fStack = 0;
+  //
+  if (fUseMC) {
+    eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
+    if (!eventHandler) { AliError("Could not retrieve MC event handler"); return; }
+    fMCEvent = eventHandler->MCEvent();
+    if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
+    fStack = fMCEvent->Stack();
+    if (!fStack) { AliError("Stack not available"); return; }
+  }
+  //
+  const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
+  if (vtxESD->GetNContributors()<1) return;
+  //
+  // ------------------RS: Put here additional event selection if needed ------->>
+  // at the moment, I am just selecting on vertex
+  TString vtxTyp = vtxESD->GetTitle();
+  if (vtxTyp.Contains("vertexer: Z")) {
+    if (vtxESD->GetDispersion()>0.04) return;
+    if (vtxESD->GetZRes()>0.25) return;
+  }
+  //
+  if (vtxESD->GetZ()<fZVertexMin || vtxESD->GetZ()>fZVertexMax) return;
+  //
+  // pile-up rejection
+  if (esd->IsPileupFromSPD(3,0.8)) return;
+  //
+  // ------------------RS: Put here additional event selection if needed -------<<
+  //
+  fVtx[0] = vtxESD->GetX();
+  fVtx[1] = vtxESD->GetY();
+  fVtx[2] = vtxESD->GetZ();
+  //
+  // do we need to initialize the field?
+  AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+  if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
+  //
+  AliStack *stack = NULL;
+  if (fUseMC) {
+    if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
+    stack = fMCEvent->Stack();
+    if (!stack) { AliError("ERROR: Could not get stack"); return; }
+    for (Int_t iPart = 0; iPart < stack->GetNtrack(); ++iPart) {
+      if (!stack->IsPhysicalPrimary(iPart)) continue;
+      TParticle* p = stack->Particle(iPart);
+      if (TMath::Abs(p->Y()) > 0.5) continue;
+      if (TMath::Abs(p->GetPdgCode()) != 211) continue;
+      TH1D *hPt = (TH1D*)fHistos->At(kHPt);
+      hPt->Fill(p->Pt());
+    }
+  }
+  //
+  InitMultReco();
+  fMultReco->Reconstruct(fRPTree, fVtx);
+  //  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
+  FillHistos(stack);
+  //
+  delete fMultReco; 
+  fMultReco = 0;
+  //
+}      
+
+//________________________________________________________________________
+void AliSPDClustTask::Terminate(Option_t *) 
+{
+  Printf("Terminating...");
+}
+
+//_________________________________________________________________________
+void AliSPDClustTask::InitMultReco()
+{
+  // create mult reconstructor
+  if (fMultReco) delete fMultReco;
+  fMultReco = new AliITSMultReconstructor();
+  fMultReco->SetCreateClustersCopy(kTRUE);
+  fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
+  fMultReco->SetNStdDev(fNStdDev);
+  fMultReco->SetPhiWindow( fDPhiWindow );
+  fMultReco->SetThetaWindow( fDThetaWindow );
+  fMultReco->SetPhiShift( fDPhiShift );
+  fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
+  fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
+  fMultReco->SetZetaOverlapCut(fZetaOverlap);
+  fMultReco->SetHistOn(kFALSE); 
+  //
+}
+
+//_________________________________________________________________________
+TObjArray* AliSPDClustTask::BookHistos()
+{
+  // book set of histos
+  //
+  TObjArray* histos = new TObjArray();
+  histos->SetOwner(kFALSE);
+  // this array contains histos stored at slots:
+  // 0-99    : tracklet related histos
+  // 100-199 : SPD1 clusters not used by tracklets
+  // 200-299 : SPD2 clusters not used by tracklets
+  // 300-399 : SPD1 clusters used by tracklets
+  // 400-499 : SPD2 clusters used by tracklets
+  //
+  int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
+  if (nEtaBins<1) nEtaBins = 1;
+  //
+  // just an example: cluster type vs eta
+  for (int iused=0;iused<2;iused++) {
+    for (int spd=0;spd<2;spd++) {
+      TH2F* h = new TH2F(Form("clType_SPD%d_%s",spd,iused ? "used":"free"),
+                        Form("clType SPD%d %s",spd,iused ? "used":"free"),
+                        nEtaBins, fEtaMin,fEtaMax,
+                        15, -0.5, 14.5);
+      histos->AddAtAndExpand(h, kHClusters+iused*200+spd*100 + kClTypevsEta);
+      TH1F* hZ = new TH1F(Form("clZ_SPD%d_%s",spd,iused ? "used":"free"),
+                         Form("clZ SPD%d %s",spd,iused ? "used":"free"),
+                         200, -15, 15);
+      histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZ);
+      TH1F* hEta = new TH1F(Form("clEta_SPD%d_%s",spd,iused ? "used":"free"),
+                           Form("clEta SPD%d %s",spd,iused ? "used":"free"),
+                           nEtaBins, fEtaMin,fEtaMax);
+      histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEta);
+
+      hZ = new TH1F(Form("clZ_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
+                   Form("clZ pions weighted SPD%d %s",spd,iused ? "used":"free"),
+                   200, -15, 15);
+      histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPionsW);
+      hEta = new TH1F(Form("clEta_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
+                     Form("clEta pions weighted SPD%d %s",spd,iused ? "used":"free"),
+                     nEtaBins, fEtaMin,fEtaMax);
+      histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPionsW);
+
+      hZ = new TH1F(Form("clZ_pions_SPD%d_%s",spd,iused ? "used":"free"),
+                   Form("clZ pions SPD%d %s",spd,iused ? "used":"free"),
+                   200, -15, 15);
+      histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPions);
+      hEta = new TH1F(Form("clEta_pions_SPD%d_%s",spd,iused ? "used":"free"),
+                     Form("clEta pions SPD%d %s",spd,iused ? "used":"free"),
+                     nEtaBins, fEtaMin,fEtaMax);
+      histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPions);
+    }
+  }
+  if (fUseMC) {
+    TH1D *hPt = new TH1D("hPt","Pions Pt spectra (MC)",fhPtIn->GetXaxis()->GetNbins(),fhPtIn->GetXaxis()->GetXbins()->GetArray());
+    histos->AddAtAndExpand(hPt, kHPt);
+  }
+  //
+  return histos;
+}
+
+//_________________________________________________________________________
+void AliSPDClustTask::FillHistos(AliStack *stack)
+{
+  // fill info on clusters, separately on associated to tracklets and singles
+  const int kUsedBit = BIT(15);
+  //
+  // 0) --------- just in case: reset kUsedBit of clusters --->>>
+  for (int ilr=0;ilr<2;ilr++) {
+    for (int icl=fMultReco->GetNClustersLayer(ilr);icl--;) {  // loop on clusters of layer
+      AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl);
+      if (!clus) {
+       AliError(Form("Failed to extract cluster %d of layer %d",icl,ilr));
+       return;
+      }
+      clus->ResetBit(kUsedBit);
+    } // loop on clusters of layer
+  } // loop on layers
+  //
+  // 1) --------- mark used clusters by kUsedBit ------------>>>
+  int ntr = fMultReco->GetNTracklets();
+  for (int itr=ntr;itr--;) {
+    Float_t *trc = fMultReco->GetTracklet(itr);
+    AliITSRecPoint *cl0 = fMultReco->GetRecPoint(0,(int)trc[AliITSMultReconstructor::kClID1]); // cluster on SPD1
+    AliITSRecPoint *cl1 = fMultReco->GetRecPoint(1,(int)trc[AliITSMultReconstructor::kClID2]); // cluster on SPD1
+    //
+    if (!cl0 || !cl1) {
+      AliError(Form("Failed to extract clusters associated with tracklet %d",itr));
+      continue;
+    }
+    cl0->SetBit(kUsedBit);
+    cl1->SetBit(kUsedBit);
+  }
+  //
+  // 2) --------- fill needed info for used and unuses clusters
+  TH2F* h2d = 0;
+  TH1F* h1d = 0;
+  //
+  for (int ilr=0;ilr<2;ilr++) { // loop on layers
+    int ncl = fMultReco->GetNClustersLayer(ilr);
+    for (int icl=ncl;icl--;) {  // loop on clusters of layer
+      AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl); // extract the cluster to access clType
+      Float_t* clInfo = fMultReco->GetClusterOfLayer(ilr,icl); // we already have in the MultReco cluster params extracted and
+      // stored in the float array, see AliITSMultReconstructor
+      //
+      float phi = clInfo[AliITSMultReconstructor::kClPh];
+      float eta = -TMath::Log(TMath::Tan(clInfo[AliITSMultReconstructor::kClTh]/2));
+      float z   = clInfo[AliITSMultReconstructor::kClZ];
+      //
+      int used = clus->TestBit(kUsedBit) ? 1:0;
+      //
+      // Annalisa, here you should fill the info on used/unused clusters
+      if (phi > 0 && phi < 1.6) {
+       h2d = (TH2F*)fHistos->At(kHClusters+used*200+ilr*100 + kClTypevsEta);
+       h2d->Fill(eta, clus->GetClusterType());
+       h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZ);
+       h1d->Fill(z);
+       h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEta);
+       h1d->Fill(eta);
+       if (stack) {
+         for(Int_t iLabel = 0; iLabel < 3; ++iLabel) {
+           Int_t clLabel = clus->GetLabel(iLabel);
+           if (clLabel <= 0) continue;
+           Int_t moLabel = FindMotherParticle(stack, clLabel);
+           if (moLabel <= 0) continue;
+           TParticle* p = stack->Particle(moLabel);
+           if (TMath::Abs(p->GetPdgCode()) != 211) continue;
+           Double_t weight = PtWeight(p->Pt());
+           h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPionsW);
+           h1d->Fill(z,weight);
+           h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPionsW);
+           h1d->Fill(eta,weight);
+           h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPions);
+           h1d->Fill(z);
+           h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPions);
+           h1d->Fill(eta);
+           break;
+         }
+       }
+      }
+    } // loop on clusters of layer
+  } // loop on layers
+  //
+}
+
+//________________________________________________________________________
+void AliSPDClustTask::SetInput(const char *filename)
+{
+  // Read pions Pt spectra
+  // from a local input file
+
+  AliInfo(Form("Reading input histograms from %s",filename));
+  if (fUseMC) {
+    TFile *fPt = TFile::Open("spectraCombine.root");
+    if (!fPt) {        AliError("Failed to open input file"); return; }
+    TList *lPt = (TList*)fPt->Get("output");
+    fhPtIn = (TH1D*)lPt->FindObject("PionsPos_MB_combine_sum")->Clone("fhPtIn");
+    fhPtIn->SetDirectory(0);
+    fPt->Close();
+  }
+}
+
+//________________________________________________________________________
+Int_t AliSPDClustTask::FindMotherParticle(AliStack* stack, Int_t i)
+{
+  if(stack->IsPhysicalPrimary(i)) return i;
+  TParticle* p = stack->Particle(i);
+  Int_t imo =  p->GetFirstMother();
+
+  //  printf("imo = %d\n",imo);
+
+  if(imo<=0)
+    {
+      AliInfo("particle is not a PhysPrim and has no mother");
+      return imo;
+    }
+  return FindMotherParticle(stack, imo);
+
+}
+
+//________________________________________________________________________
+Double_t AliSPDClustTask::PtWeight(Double_t pt)
+{
+  if (pt > 0.35)
+    return 1.;
+  else
+    return (1.71-0.71*pt/0.35);
+}
diff --git a/PWGPP/forward/SPDClustTask/AliSPDClustTask.h b/PWGPP/forward/SPDClustTask/AliSPDClustTask.h
new file mode 100755 (executable)
index 0000000..07f18dd
--- /dev/null
@@ -0,0 +1,110 @@
+#ifndef ALISPDCLUSTTASK_H
+#define ALISPDCLUSTTASK_H
+
+///////////////////////////////////////////////////////////////////////////
+// Class AliSPDClustTask                                            //
+// Analysis task to produce data and MC histos needed for tracklets      //
+// dNdEta extraction in multiple bins in one go                          //
+// Author:  ruben.shahoyan@cern.ch                                       //
+///////////////////////////////////////////////////////////////////////////
+
+class TH1F; 
+class TH2F;
+class AliESDEvent;
+class TList;
+
+class AliMCParticle;
+class AliITSMultReconstructor;
+class AliESDTrackCuts;
+
+#include "../ITS/AliITSsegmentationSPD.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliTriggerAnalysis.h" 
+#include <TMath.h>
+
+class AliSPDClustTask : public AliAnalysisTaskSE {
+ public:
+  enum {kHTracklets=0, kHPt = 50, kHClusters=100, kClTypevsEta=0, kClZ=1, kClEta=2, kClZPions=3, kClEtaPions=4, kClZPionsW=5, kClEtaPionsW=6}; // to facilitated access to histos, see BookHistos
+  //
+  AliSPDClustTask(const char *name = "AliSPDClustTask");
+  virtual ~AliSPDClustTask(); 
+  
+  virtual void  UserCreateOutputObjects();
+  virtual void  UserExec(Option_t *option);
+  virtual void  Terminate(Option_t *);
+  void       SetUseMC(Bool_t mc = kFALSE)              {fUseMC = mc;}
+  TObjArray* BookHistos();
+  void       FillHistos(AliStack *stack);
+  // RS
+  void       SetNStdDev(Float_t f=1.)           {fNStdDev = f<1e-5 ? 1e-5:f;}
+  void       SetScaleDThetaBySin2T(Bool_t v=kFALSE) {fScaleDTBySin2T = v;}
+  void       SetCutOnDThetaX(Bool_t v=kFALSE)   {fCutOnDThetaX = v;}
+  void       SetPhiWindow(float w=0.08)         {fDPhiWindow   = w<1e-5 ? 1e-5:w;}
+  void       SetThetaWindow(float w=0.025)      {if (w<0) fCutOnDThetaX=kTRUE; fDThetaWindow = TMath::Abs(w)<1e-5 ? 1e-5:TMath::Abs(w);}
+  void       SetPhiShift(float w=0.0045)        {fDPhiShift = w;}
+  void       SetPhiOverlapCut(float w=0.005)    {fPhiOverlapCut = w;}
+  void       SetZetaOverlapCut(float w=0.05)    {fZetaOverlap = w;}
+  void       SetRemoveOverlaps(Bool_t w=kFALSE) {fRemoveOverlaps = w;}
+  //
+  void       SetDPhiSCut(Float_t c=0.06)        {fDPhiSCut = c;}
+  void       SetNStdCut(Float_t c=1.0)          {fNStdCut = c;}
+  //
+  void       SetEtaCut(Float_t etaCut)          {fEtaMax = TMath::Abs(etaCut); fEtaMin= -fEtaMax;}
+  void       SetEtaMin(Float_t etaMin)          {fEtaMin = etaMin;}
+  void       SetEtaMax(Float_t etaMax)          {fEtaMax = etaMax;}
+  void       SetZVertexMin(Float_t z)           {fZVertexMin = z;}
+  void       SetZVertexMax(Float_t z)           {fZVertexMax = z;}
+  //
+  void       SetInput(const char *filename);
+  Int_t      FindMotherParticle(AliStack* stack, Int_t i);
+  Double_t   PtWeight(Double_t pt);
+  //
+ protected:
+  void       InitMultReco();
+  void       FillClusterInfo();
+  //
+ protected:
+  TList*       fOutput;                   // output list send on output slot 1 
+  //
+  Bool_t       fUseMC; 
+  TObjArray*   fHistos;                   //! histos array
+  Float_t      fVtx[3];                   //! event vertex
+  //
+  // Settings for the reconstruction
+  // tracklet reco settings
+  Float_t      fEtaMin;                    // histos filled only for this eta range
+  Float_t      fEtaMax;                    // histos filled only for this eta range
+  Float_t      fZVertexMin;                // min Z vtx to process
+  Float_t      fZVertexMax;                // max Z vtx to process
+  //
+  Bool_t       fScaleDTBySin2T;            // request dTheta scaling by 1/sin^2(theta)
+  Bool_t       fCutOnDThetaX;              // if true, apart from NStdDev cut apply also the cut on dThetaX
+  Float_t      fNStdDev;                   // cut on weighted distance
+  Float_t      fDPhiWindow;                // max dPhi
+  Float_t      fDThetaWindow;              // max dTheta
+  Float_t      fDPhiShift;                 // mean bend
+  Float_t      fPhiOverlapCut;             // overlaps cut in phi
+  Float_t      fZetaOverlap;               // overlaps cut in Z
+  Bool_t       fRemoveOverlaps;            // request overlaps removal
+  //
+  Float_t      fDPhiSCut;                  // cut on signal dphiS
+  Float_t      fNStdCut;                   // cut on signal weighted distance
+  //
+  AliITSMultReconstructor *fMultReco;              //! mult.reco object
+  TTree*       fRPTree;                    //! tree of recpoints
+  AliStack*    fStack;                     //! MC stack
+  AliMCEvent*  fMCEvent;                   //! MC Event
+  Float_t      fESDVtx[3];                 //  ESD vertex
+  //
+  Bool_t fDontMerge;                       // no merging requested
+  //
+  TH1D*  fhPtIn; // Input histogram containing the pion spectra
+ private:    
+  AliSPDClustTask(const AliSPDClustTask&); // not implemented
+  AliSPDClustTask& operator=(const AliSPDClustTask&); // not implemented 
+  
+  ClassDef(AliSPDClustTask, 1);  
+};
+
+
+#endif
diff --git a/PWGPP/forward/SPDClustTask/AnalysisSPDClustTask.C b/PWGPP/forward/SPDClustTask/AnalysisSPDClustTask.C
new file mode 100644 (file)
index 0000000..05f4ad7
--- /dev/null
@@ -0,0 +1,265 @@
+Bool_t needRecPoints = kTRUE;
+
+void AnalysisSPDClustTask
+(
+ TString dataset="/alice/sim/LHC11d2_000119161",
+ TString outFName = "trbg.root",
+ Int_t   nEvents   = -1,
+ Float_t etaMin     =-1.5,        // min eta range to fill in histos
+ Float_t etaMax     = 1.5,        // max eta range to fill in histos
+ Float_t zMin       = -5,         // process events with Z vertex min
+ Float_t zMax       =  5,         //                     max positions
+ Float_t cutSigNStd  = 1.,        // cut on weighed distance used to extract signal
+ Float_t cutSigDPhiS = -1,        // cut on dPhi-phiBent used to extract signal (if negative -> dphi*sqrt(cutSigNStd)
+ Bool_t  useMC  = kTRUE,          // fill MC info (doRec=kTRUE)
+ // specific parameters for reconstruction
+ Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
+ float  nStdDev     = 1.,         // number of st.dev. for tracklet cut to keep
+ float  dphi        = 0.06,        // dphi window (sigma of tracklet cut)
+ float  dtht        = 0.025,       // dtheta .... (if negative, abs will be used with additional cut on |dthetaX|, apart from w.distance
+ float  phishift    = 0.0045,      // bending shift
+ Bool_t remOvl      = kTRUE,       
+ float  ovlPhiCut   = 0.005, 
+ float  ovlZetaCut  = 0.05,
+ Int_t  nEventsSkip = 0
+ )
+{
+  //  
+  if (cutSigDPhiS<0) cutSigDPhiS = TMath::Sqrt(cutSigNStd)*dphi;
+  //
+  printf("Start Analysis for %s, max %d Events skipping %d, Event Cuts: %.1f<eta<%.1f, %.2f<Zv<%.2f\n",
+        dataset.Data(),nEvents,nEventsSkip,etaMin,etaMax,zMin,zMax);
+  printf("Tracklet cuts: dPhi:%.3f dTheta:%.3f phiShift:%.4f | Keep %.1f NstDev\n"
+        "Scale dTheta: %s | Signal Selection: NstDev:%.1f, dPhiS: %.3f\n", 
+        dphi,dtht,phishift,nStdDev,scaleDTheta ? "ON":"OFF",
+        cutSigNStd,cutSigDPhiS);
+  //
+  if (nEvents<0) nEvents = int(1e9);
+  TString format = GetFormatFromDataSet(dataset);
+  //
+  // ALICE stuff
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) mgr = new AliAnalysisManager("Test train");
+  //
+  InputHandlerSetup(format,useMC);
+  // compile our task
+  gProof->Load("AliSPDClustTask.cxx++");
+  //
+  // load and run AddTask macro
+  //  gROOT->LoadMacro("AddMultTaskTrackletMulti.C");
+  //
+  // create our task
+  
+  AliSPDClustTask *task = new AliSPDClustTask("AliSPDClustTask");
+  // create output container
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist", TList::Class(),AliAnalysisManager::kOutputContainer,outFName);
+  //  coutput1->SetSpecialOutput();
+  // add our task to the manager
+  mgr->AddTask(task);
+
+  // finaly connect input and output
+  mgr->ConnectInput(task, 0,  mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task,1,coutput1);
+  //  mgr->SetSpecialOutputLocation("root://alicers01.cern.ch//tmp/output/");
+  //
+  //
+  task->SetUseMC(useMC);
+  //
+  task->SetEtaMin(etaMin);
+  task->SetEtaMax(etaMax);
+  task->SetZVertexMin(zMin);
+  task->SetZVertexMax(zMax);
+  //
+  task->SetDPhiSCut(cutSigDPhiS);
+  task->SetNStdCut(cutSigNStd);
+  //
+  task->SetScaleDThetaBySin2T(scaleDTheta);
+  task->SetNStdDev(nStdDev);
+  task->SetPhiWindow(dphi);
+  task->SetThetaWindow(dtht);
+  task->SetPhiShift(phishift);
+  task->SetPhiOverlapCut(ovlPhiCut);
+  task->SetZetaOverlapCut(ovlZetaCut);
+  task->SetRemoveOverlaps(remOvl);
+  //
+  task->SetInput("spectraCombine.root");
+  //
+  printf("Requesting physics selection in %s mode\n",useMC ? "MC":"Data");
+  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  //  /*
+  //gROOT->ProcessLine(".L AddTaskPhysicsSelection.C");
+  AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(useMC);
+  task->SelectCollisionCandidates();//AliVEvent::kMB);
+  //
+  //  */
+  // Run analysis
+  mgr->InitAnalysis();
+  // process dataset  
+  mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip); 
+  //
+  TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_"; 
+  evstCmd += outFName;  evstCmd += " \nfi";
+  gSystem->Exec( evstCmd.Data() );
+  
+}
+
+
+TString GetFormatFromDataSet(TString dataset) {
+  
+//   Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
+  TString dsTreeName;
+  if (dataset.Contains("#")) {
+    Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
+    dsTreeName=dataset(dataset.Last('#'),dataset.Length());
+  } else {
+    Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
+    TFileCollection *ds = gProof->GetDataSet(dataset.Data());
+    if (!ds) {
+      Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
+      return "";
+    }
+    dsTreeName = ds->GetDefaultTreeName();
+  }
+
+  if (dsTreeName.Contains("esdTree")) {
+    Info("runAAF.C","ESD input format detected ...");
+    return "ESD";
+  } else if (dsTreeName.Contains("aodTree"))  {
+    Info("runAAF.C","AOD input format detected ...");
+    return "AOD";
+  } else {
+    Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
+    Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
+  }
+  
+  return "";
+}
+
+Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
+{
+  format.ToLower();
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+
+  AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
+
+  if (cin) return;
+
+  if (!format.CompareTo("esd"))
+  {
+    AliESDInputHandler *esdInputHandler = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+    if (!esdInputHandler)
+    {
+      Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
+      if (needRecPoints)
+       esdInputHandler = new AliESDInputHandlerRP();
+      else 
+       esdInputHandler = new AliESDInputHandler();
+      //
+      mgr->SetInputEventHandler(esdInputHandler);
+    }
+    //
+    if (useKine)
+    {
+      AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+
+      if (!mcInputHandler)
+      {
+        Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
+        AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
+        mgr->SetMCtruthEventHandler(mcInputHandler);
+      }
+    }
+
+  }
+  else if (!format.CompareTo("aod"))
+  {
+    AliAODInputHandler *aodInputHandler = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+    if (!aodInputHandler)
+    {
+      Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
+      aodInputHandler = new AliAODInputHandler();
+      mgr->SetInputEventHandler(aodInputHandler);
+    }
+  }
+  else
+  {
+    Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
+    return kFALSE;
+  }
+
+  return kTRUE;
+}
+
+void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
+                    float zMin, float zMax, float zMixBinSz)
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) return;
+  int bufferSize = 1;
+  AliESDInputHandlerRP *esdH = dynamic_cast<AliESDInputHandlerRP*>(mgr->GetInputEventHandler());
+  if (!esdH) return;
+  //
+  AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
+  esdMixH->SetInputHandlerForMixing(esdH);
+  AliMixEventPool *evPool = new AliMixEventPool("Pool");
+  AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
+  AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
+  //  evPool->AddCut(tracklets);
+  evPool->AddCut(zvertex);
+  //evPool->Init();
+  evPool->Print();
+  esdMixH->SetEventPool(evPool);
+  esdH->SetMixingHandler(esdMixH);
+}
+
+void AddPhysicsSelection(Bool_t isMC)
+{
+  // physics selection a la Michele
+  if(!isMC) {
+    //AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
+    //    physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
+    /*
+    physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
+    physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
+    */
+    //
+    //    physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
+    //    physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
+    //
+    // This are needed only to fill the statistics tables
+    //    physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
+    //    physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
+    //    physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
+    //
+    /*
+    physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
+    //
+    physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
+    //
+    */
+    /*
+    //
+    physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
+    //
+    physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
+    physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
+    */
+  } 
+  // if you use the following line, your task only gets the selected events
+  //  task->SelectCollisionCandidates(AliVEvent::kUserDefined);
+  //  task->SelectCollisionCandidates();
+  //
+  //Alternatively, in the UserExec of your task:
+  //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
+  //
+}
diff --git a/PWGPP/forward/SPDClustTask/FitOutput.root b/PWGPP/forward/SPDClustTask/FitOutput.root
new file mode 100644 (file)
index 0000000..9903126
Binary files /dev/null and b/PWGPP/forward/SPDClustTask/FitOutput.root differ
diff --git a/PWGPP/forward/SPDClustTask/runAAFSPDTask.C b/PWGPP/forward/SPDClustTask/runAAFSPDTask.C
new file mode 100644 (file)
index 0000000..47d9d0a
--- /dev/null
@@ -0,0 +1,98 @@
+//
+void runAAFSPDTask(TString dataset="/alice/sim/LHC11d2_000119161",
+                  TString outFName = "spdres.root",
+                  Int_t   nEvents    = -1,//3000,
+                  Float_t etaMin     =-1.5,        // min eta range to fill in histos
+                  Float_t etaMax     = 1.5,        // max eta range to fill in histos
+                  Float_t zMin       = -5,         // process events with Z vertex min
+                  Float_t zMax       =  5,         //                     max positions
+                  //
+                  Float_t cutSigNStd  = 1.,       // cut on weighed distance used to extract signal
+                  Float_t cutSigDPhiS = -1,        // cut on dPhi-phiBent used to extract signal (if negative -> dphi*sqrt(cutSigNStd)
+                  Bool_t  useMC  = kTRUE,          // fill MC info (doRec=kTRUE)
+                  //
+                  // specific parameters for reconstruction
+                  float  phiRot      = 3.14159e+00, // angle for bg. generation with rotation
+                  Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
+                  float  nStdDev     = 1.,         // number of st.dev. for tracklet cut to keep
+                  float  dphi        = 0.06,        // dphi window (sigma of tracklet cut)
+                  float  dtht        = 0.025,       // dtheta .... (if negative, abs will be used with additional cut on |dthetaX|, apart from w.distance
+                  float  phishift    = 0.0045,      // bending shift
+                  Bool_t remOvl      = kTRUE,       
+                  float  ovlPhiCut   = 0.005, 
+                  float  ovlZetaCut  = 0.05,
+                  Int_t  nEventsSkip = 0,
+                  //
+                  TString alirootVer = "VO_ALICE@AliRoot::v5-03-15-AN",
+                  TString rootVer    = "VO_ALICE@ROOT::v5-33-02a",
+                  //
+                  TString proofCluster="cheshkov@alice-caf.cern.ch"
+                  ) 
+{ 
+  //  
+  Bool_t runLocal = kFALSE;//kTRUE; // true only for local test mode
+  if (runLocal) {
+    //    dataset = "/default/shahoian/test_pp";//"/default/shahoian/test";
+    dataset = "default/shahoian/test_sim_lhc11b1a";
+    proofCluster = "";
+    alirootVer = "AliRootProofLite";
+    nEvents = 500;
+  }
+  //
+  if ((!dataset.Contains("sim")) && useMC) {
+    printf("Running with read data dataset, switching OFF useMC\n");
+    useMC = kFALSE;
+  }
+  //
+  printf("Requested: %s %s\n",alirootVer.Data(), rootVer.Data());
+  printf("Output expected in %s\n",outFName.Data());
+  //
+  gEnv->SetValue("XSec.GSI.DelegProxy","2");
+  //
+  TString alirootMode="REC";//"REC";
+  TString extraLibs = "ITSrec:CDB:Geom:"; // not needed in default aliroot mode
+  //extraLibs+= "ANALYSIS:ANALYSISalice";
+  extraLibs+= "ANALYSIS:OADB:ANALYSISalice:EventMixing";
+  TList *list = new TList();
+  // sets $ALIROOT_MODE on each worker to let proof to know to run in special mode
+  list->Add(new TNamed("ALIROOT_MODE"      , alirootMode.Data()));
+  list->Add(new TNamed("ALIROOT_EXTRA_LIBS", extraLibs.Data()));
+  list->Add(new TNamed("ALIROOT_EXTRA_INCLUDES", "ITS:include"));
+  list->Add(new TNamed("ALIROOT_ENABLE_ALIEN","1"));
+  //
+  //REM: same version of AliRoot on client!!!!! Otherwise error!! 
+  TProof::Mgr(proofCluster.Data())->SetROOTVersion(rootVer.Data());
+  //  TProof::Open(proofCluster.Data());//,"workers=10x");
+  TProof::Open(proofCluster.Data(),"workers=1x");
+  if (!gProof) {
+    Error("runAAFMulti.C","Connection to AF failed.");
+    return;
+  }
+  //  gProof->Exec("TObject *o = gEnv->GetTable()->FindObject(\"Proof.UseMergers\");"
+  //          "gEnv->GetTable()->Remove(o);", kTRUE);
+  gProof->SetParameter("PROOF_UseMergers", 0);
+  // Lets enable aliroot + extra libs on proof cluster
+  if (runLocal) gProof->UploadPackage(alirootVer.Data());
+  gProof->EnablePackage(alirootVer.Data(), list);
+  //  gProof->EnablePackage(alirootVer.Data());
+  //
+  if (runLocal) {
+    Int_t numWorkers = gProof->GetParallel();
+    if (numWorkers<1) {printf("No workers\n"); return;}
+    gProof->SetParameter("PROOF_PacketizerStrategy", (Int_t)0);
+    int frac = (Int_t) 5 / numWorkers;
+    if (frac<1) frac = 1;
+    gProof->SetParameter("PROOF_PacketAsAFraction", frac);
+  }
+  //
+  gROOT->LoadMacro(Form("%s/AnalysisSPDClustTask.C",gSystem->pwd()));
+  TStopwatch sw;
+  sw.Start();
+  AnalysisSPDClustTask(dataset,outFName,nEvents,etaMin,etaMax,zMin,zMax,
+                      cutSigNStd,cutSigDPhiS,useMC,
+                      scaleDTheta,nStdDev,dphi,dtht,
+                      phishift,remOvl,ovlPhiCut,ovlZetaCut,nEventsSkip);
+  //
+  sw.Stop();
+  sw.Print();
+}
diff --git a/PWGPP/forward/SPDClustTask/spectraCombine.root b/PWGPP/forward/SPDClustTask/spectraCombine.root
new file mode 100644 (file)
index 0000000..bc33fea
Binary files /dev/null and b/PWGPP/forward/SPDClustTask/spectraCombine.root differ