Fix in HaveCommonParent routine, updated version of AliTrackletTaskUni
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Sep 2011 15:33:33 +0000 (15:33 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Sep 2011 15:33:33 +0000 (15:33 +0000)
PWG0/multVScentPbPb/AddMultTaskTrackletUni.C [new file with mode: 0755]
PWG0/multVScentPbPb/AliTrackletTaskMulti.cxx
PWG0/multVScentPbPb/AliTrackletTaskMultipp.cxx
PWG0/multVScentPbPb/AliTrackletTaskUni.cxx
PWG0/multVScentPbPb/AliTrackletTaskUni.h
PWG0/multVScentPbPb/MyAnalysisMacroUni.C [new file with mode: 0755]
PWG0/multVScentPbPb/runAAFbgpp.C [new file with mode: 0755]

diff --git a/PWG0/multVScentPbPb/AddMultTaskTrackletUni.C b/PWG0/multVScentPbPb/AddMultTaskTrackletUni.C
new file mode 100755 (executable)
index 0000000..8d36f75
--- /dev/null
@@ -0,0 +1,22 @@
+
+AliTrackletTaskUni* AddMultTaskTrackletUni(const char* outName="tracklet.root", TString nomergeDir="")
+{
+  // create manager
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) mgr = new AliAnalysisManager("My test train");
+  // create our task
+  AliTrackletTaskUni *task = new AliTrackletTaskUni("AliTrackletTaskUni");
+  task->SetDontMerge(!nomergeDir.IsNull());
+  // create output container
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist", TList::Class(),AliAnalysisManager::kOutputContainer,outName);
+  if (!nomergeDir.IsNull()) 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);
+  if (!nomergeDir.IsNull()) mgr->SetSpecialOutputLocation(nomergeDir.Data()); //"root://alicers01.cern.ch//tmp/myoutput/");
+  //
+  return task;
+}
index ab22395..6e17d6a 100755 (executable)
@@ -1246,18 +1246,25 @@ Bool_t AliTrackletTaskMulti::HaveCommonParent(const float* clLabs0,const float*
   int npars[2]={0,0};
   const float *labs[2] = {clLabs0,clLabs1};
   int ntr = fStack->GetNtrack();
+  //  printf("\nHaveCommonParent \n");
   for (int il=0;il<2;il++) {
+    
     for (int ilb=0;ilb<3;ilb++) {
       int lbl = (int)labs[il][ilb];
-      if (lbl<0 || lbl>=ntr) continue;
+      if (lbl<2 || lbl>=ntr) continue;
       //
       while (npars[il]<kMaxPar-1) {
-       pars[il][ npars[il]++ ] = lbl;
        TParticle* part = fStack->Particle(lbl);
        if (!part) break;
+       int code = TMath::Abs(part->GetPdgCode());
+       int q = (int)TMath::Abs(part->GetPDG()->Charge());
+       if (code==21 || code<10 || q==1 || q==2 || q==4 ) break;
+       //printf("%d/%d/%d: %4d (%d)%s|",il,ilb,npars[il],lbl,part->GetStatusCode(),part->GetName());
+       pars[il][ npars[il]++ ] = lbl;
        lbl = part->GetFirstMother();
        if (lbl<1 || lbl>=ntr) break;
       }
+      //      printf("\n");
     }
   }
   // compare array of parents
index a324bd3..c9f0eb5 100755 (executable)
@@ -1258,20 +1258,27 @@ Bool_t AliTrackletTaskMulti::HaveCommonParent(const float* clLabs0,const float*
   int npars[2]={0,0};
   const float *labs[2] = {clLabs0,clLabs1};
   int ntr = fStack->GetNtrack();
+  //  printf("\nHaveCommonParent \n");
   for (int il=0;il<2;il++) {
+    
     for (int ilb=0;ilb<3;ilb++) {
       int lbl = (int)labs[il][ilb];
-      if (lbl<0 || lbl>=ntr) continue;
+      if (lbl<2 || lbl>=ntr) continue;
       //
       while (npars[il]<kMaxPar-1) {
-       pars[il][ npars[il]++ ] = lbl;
        TParticle* part = fStack->Particle(lbl);
        if (!part) break;
+       int code = TMath::Abs(part->GetPdgCode());
+       int q = (int)TMath::Abs(part->GetPDG()->Charge());
+       if (code==21 || code<10 || q==1 || q==2 || q==4 ) break;
+       //printf("%d/%d/%d: %4d (%d)%s|",il,ilb,npars[il],lbl,part->GetStatusCode(),part->GetName());
+       pars[il][ npars[il]++ ] = lbl;
        lbl = part->GetFirstMother();
        if (lbl<1 || lbl>=ntr) break;
       }
+      //      printf("\n");
     }
-  }
+  } 
   // compare array of parents
   for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE;
   return kFALSE;
index a316baa..0fdf9ea 100755 (executable)
@@ -39,7 +39,7 @@
 #include "AliESDEvent.h"  
 #include "AliESDInputHandler.h"
 #include "AliESDInputHandlerRP.h"
-#include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h"
+//#include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h"
 #include "AliCDBPath.h"
 #include "AliCDBManager.h"
 #include "AliCDBEntry.h"
@@ -61,7 +61,6 @@
 #include "AliLog.h"
 
 #include "AliPhysicsSelection.h"
-#include "AliESDCentrality.h" 
 #include "AliTrackletTaskUni.h"
 #include "AliITSMultRecBg.h"
 #include "AliGenEventHeader.h"
@@ -205,11 +204,13 @@ AliTrackletTaskUni::AliTrackletTaskUni(const char *name)
   fHistosTrRcblSec(0),
   fHistosCustom(0),
 //
-  fEtaCut(3.0),
+  fEtaMin(-3.0),
+  fEtaMax(3.0),
   fZVertexMin(-20),
   fZVertexMax( 20),
   fMultCutMin(0),
   fMultCutMax(99999),
+  fMCV0Scale(0.7520),
 //
   fScaleDTBySin2T(kFALSE),
   fCutOnDThetaX(kFALSE),
@@ -227,7 +228,8 @@ AliTrackletTaskUni::AliTrackletTaskUni(const char *name)
   fRPTree(0),
   fRPTreeMix(0),
   fStack(0),
-  fMCEvent(0)
+  fMCEvent(0),
+  fDontMerge(kFALSE)    
   /*
   ,
   fTrigger(AliTriggerAnalysis::kAcceptAll),
@@ -284,33 +286,42 @@ AliTrackletTaskUni::~AliTrackletTaskUni()
 void AliTrackletTaskUni::UserCreateOutputObjects() 
 {
   //
+
+  if (fDontMerge) {
+    OpenFile(1);
+    printf("No merging will be done\n");
+  }
+
   fOutput = new TList();
   fOutput->SetOwner(); 
   //
   AliCDBManager *man = AliCDBManager::Instance();
   if (fUseMC) {
+    printf("Loading MC geometry\n");
     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");
+      AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,-1,-1);
       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
-      if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
+      if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,-1,-1)) AliFatal("Failed to misalign geometry");
     }
     else {
       // old geom
-      AliCDBEntry*  obj = man->Get("GRP/Geometry/Data");
+      AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,7,-1);
       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
     }
   }
   else {
-    man->SetDefaultStorage("raw://"); man->SetRun(137045);
-    AliCDBEntry*  obj = man->Get("GRP/Geometry/Data");
+    printf("Loading Raw geometry\n");
+    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",137045,8,-1)) AliFatal("Failed to misalign geometry");
+    if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
   }
   //
+  printf("Geometry loaded\n");  
   // Create histograms
   //---------------------------------------------Standard histos per tracklet type--->>
   UInt_t hPattern = 0xffffffff;
@@ -359,14 +370,59 @@ void AliTrackletTaskUni::UserCreateOutputObjects()
 }
 
 //________________________________________________________________________
+void AliTrackletTaskUni::RegisterStat() 
+{
+  static Bool_t done = kFALSE;
+  if (done) return;
+  TH1* hstat;
+  Printf("Registering used settings");
+  TList *lst = dynamic_cast<TList*>(GetOutputData(1));
+  if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
+    // fill used settings
+    hstat->Fill(kDPhi,fDPhiWindow);
+    hstat->Fill(kDTht,fDThetaWindow);
+    hstat->Fill(kNStd,fNStdDev);
+    hstat->Fill(kPhiShift,fDPhiShift);
+    hstat->Fill(kThtS2,fScaleDTBySin2T);  
+    hstat->Fill(kThtCW,fCutOnDThetaX);  
+    hstat->Fill(kPhiOvl,fPhiOverlapCut);
+    hstat->Fill(kZEtaOvl,fZetaOverlap);
+    hstat->Fill(kNoOvl,fRemoveOverlaps);
+    //
+    hstat->Fill(kPhiRot,fPhiRot);
+    hstat->Fill(kInjScl,fInjScale);
+    hstat->Fill(kEtaMin,fEtaMin);
+    hstat->Fill(kEtaMax,fEtaMax);
+    hstat->Fill(kZVMin,fZVertexMin);
+    hstat->Fill(kZVMax,fZVertexMax);
+    hstat->Fill(kTrcMin,fMultCutMin);    
+    hstat->Fill(kTrcMax,fMultCutMax);    
+    hstat->Fill(kMCV0Scale, fMCV0Scale);
+    //
+    hstat->Fill(kOneUnit,1.);    
+  }
+  else Printf("Did not find stat histo");
+  done = kTRUE;
+  //
+}
+
+
+//________________________________________________________________________
 void AliTrackletTaskUni::UserExec(Option_t *) 
 {
   // Main loop
   //
+  if (fDontMerge) RegisterStat();
+  Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
+
   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
   fRPTree = fRPTreeMix = 0;
-  AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
-  if (!handRP) { printf("No RP handler\n"); return; }
+  AliESDInputHandler *handler = (AliESDInputHandler*)anMan->GetInputEventHandler();
+  AliESDInputHandlerRP *handRP = 0;
+  if (needRecPoints) {
+    handRP = (AliESDInputHandlerRP*)handler;
+    if (!handRP) { printf("No RP handler\n"); return; }
+  }
   AliESDEvent *esd  = handRP->GetEvent();
   if (!esd) { printf("No AliESDEvent\n"); return; }
   //
@@ -393,17 +449,28 @@ void AliTrackletTaskUni::UserExec(Option_t *)
     }
   }
   */
-
+  TH1F* hstat = (TH1F*)fHistosCustom->UncheckedAt(kHStat);
+  //
+  hstat->Fill(kEvTot0); // RS
   const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
   if (vtxESD->GetNContributors()<1) return;
-  if (vtxESD->GetDispersion()>0.04) return;
-  if (vtxESD->GetZRes()>0.25) return;
+  TString vtxTyp = vtxESD->GetTitle();
+  if (vtxTyp.Contains("vertexer: Z")) {
+    if (vtxESD->GetDispersion()>0.04) return;
+    if (vtxESD->GetZRes()>0.25) return;
+  }
+  // pile-up rejection
+  if (esd->IsPileupFromSPD(3,0.8)) {
+    hstat->Fill(kEvTotPlp); // RS
+    return;
+  }
+  //
   const AliMultiplicity* multESD = esd->GetMultiplicity();
+  /*
   const AliESDVertex* vtxESDTPC = esd->GetPrimaryVertexTPC();
   if (vtxESDTPC->GetNContributors()<1 ||
       vtxESDTPC->GetNContributors()<(-10.+0.25*multESD->GetNumberOfITSClusters(0))) return;
-  //
-  TH1F* hstat = (TH1F*)fHistosCustom->UncheckedAt(kHStat);
+  */
   //
   hstat->Fill(kEvTot); // RS
   //
@@ -415,25 +482,13 @@ void AliTrackletTaskUni::UserExec(Option_t *)
   //
   //------------------------------------------------------
   // ZDC cut
+   //------------------------------------------------------
   AliESDZDC *esdZDC = esd->GetESDZDC();
-  // --- ZDC offline trigger ---
-  // Returns if ZDC triggered, based on TDC information
-  Bool_t tdc[32] = {kFALSE};
-  for(Int_t itdc=0; itdc<32; itdc++){
-    for(Int_t i=0; i<4; i++){
-      if (0.025*esdZDC->GetZDCTDCData(itdc, i) != 0){
-       tdc[itdc] = kTRUE;
-      }
-    }
+  float zdcEnergy=0,zemEnergy=0;
+  if (esdZDC) {
+    zdcEnergy = (esdZDC->GetZDCN1Energy() + esdZDC->GetZDCP1Energy() + esdZDC->GetZDCN2Energy()+ esdZDC->GetZDCP2Energy());
+    zemEnergy = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.; 
   }
-  Bool_t zdcNA = tdc[12];
-  Bool_t zdcNC = tdc[10];
-  Bool_t zdcPA = tdc[13];
-  Bool_t zdcPC = tdc[11];
-  //
-  Bool_t zdcC= ((zdcPA) || (zdcNA));
-  Bool_t zdcA= ((zdcPC) || (zdcNC));
-  if (!zdcC || !zdcA) return;
   //-----------------------------------------------------
   Float_t multV0A=0,multV0C=0;
   AliESDVZERO* esdV0 = esd->GetVZEROData();
@@ -442,7 +497,7 @@ void AliTrackletTaskUni::UserExec(Option_t *)
     multV0C = esdV0->GetMTotV0C();
   }
   if (fUseMC) {
-    const double v0Scale = 0.75108;
+    const double v0Scale = fMCV0Scale;
     multV0A *= v0Scale;
     multV0C *= v0Scale;    
   }
@@ -467,7 +522,7 @@ void AliTrackletTaskUni::UserExec(Option_t *)
   //
   //  if((multESD->GetNumberOfTracklets() < fMultCutMin) || (multESD->GetNumberOfTracklets() > fMultCutMax)) return;
   //
-  printf("Multiplicity from ESD:\n");
+  //  printf("Multiplicity from ESD:\n");
   //  multESD->Print();
   //
   AliMCEventHandler* eventHandler = 0;
@@ -483,8 +538,10 @@ void AliTrackletTaskUni::UserExec(Option_t *)
     if (!fStack) { printf("Stack not available\n"); return; }
   }
   //
-  fRPTree = handRP->GetTreeR("ITS");
-  if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
+  if (needRecPoints) {
+    fRPTree = handRP->GetTreeR("ITS");
+    if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
+  }
   //
   //
   // registed Ntracklets and ZVertex of the event
@@ -498,13 +555,18 @@ void AliTrackletTaskUni::UserExec(Option_t *)
   ((TH1F*)fHistosCustom->UncheckedAt(kHV0))->Fill(multV0A+multV0C);
   //
   // normal reconstruction
+  static int prnRec = 10;
+  static int prnInj = 10;
+  //
   hstat->Fill(kEvProcData);
   if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done
     InitMultReco();
     fMultReco->Run(fRPTree, vtxf);
-    printf("Multiplicity Reconstructed:\n");
     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
-    if (mlt) mlt->Print();
+    if (mlt && (--prnRec)>0) {
+      printf("Multiplicity Reconstructed: %d\n",prnRec);
+      mlt->Print();
+    }
     if (GetDoNormalReco()) FillHistos(kData,mlt);
     //
   }
@@ -514,10 +576,13 @@ void AliTrackletTaskUni::UserExec(Option_t *)
   if (GetDoInjection()) {
     if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above
     fMultReco->SetRecType(AliITSMultRecBg::kBgInj);
-    fMultReco->Run(fRPTree, vtxf);
-    printf("Multiplicity from Injection:\n");
+    fMultReco->Run(fRPTree, vtxf);    
     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
-    if (mlt) mlt->Print();
+    if (mlt && (--prnInj)>0) {
+      printf("Multiplicity from Injection: %d\n",prnInj);
+      mlt->Print();
+    }
+    // if (mlt) mlt->Print();
     hstat->Fill(kEvProcInj);
     FillHistos(kBgInj,mlt);
   }
@@ -528,14 +593,13 @@ void AliTrackletTaskUni::UserExec(Option_t *)
     fMultReco->SetRecType(AliITSMultRecBg::kBgRot);
     fMultReco->SetPhiRotationAngle(fPhiRot);
     fMultReco->Run(fRPTree, vtxf);
-    printf("Multiplicity from Rotation:\n");
     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
-    if (mlt) mlt->Print();
     hstat->Fill(kEvProcRot);
     FillHistos(kBgRot,mlt);
   }
   //
   if (GetDoMixing()) {
+    /*
     AliMixEventInputHandler* handToMix = (AliMixEventInputHandler*)handRP->MixingHandler();
     if (!handToMix) { printf("No Mixing handler\n"); return; }
     handToMix->GetEntry();
@@ -560,6 +624,8 @@ void AliTrackletTaskUni::UserExec(Option_t *)
     if (mlt) mlt->Print();
     hstat->Fill(kEvProcMix);
     FillHistos(kBgMix,mlt);
+    */
+    AliFatal("Mixing is outphased");
     //
   }
   // =============================================================================<<<
@@ -572,34 +638,7 @@ void AliTrackletTaskUni::UserExec(Option_t *)
 void AliTrackletTaskUni::Terminate(Option_t *) 
 {
   Printf("Terminating...");
-  TH1* hstat;
-  TList *lst = dynamic_cast<TList*>(GetOutputData(1));
-  printf("Term: %p %p %p\n",fOutput,lst,fHistosCustom);
-  if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
-    Info("Terminate","Registering used settings");
-    // fill used settings
-    hstat->Fill(kDPhi,fDPhiWindow);
-    hstat->Fill(kDTht,fDThetaWindow);
-    hstat->Fill(kNStd,fNStdDev);
-    hstat->Fill(kPhiShift,fDPhiShift);
-    hstat->Fill(kThtS2,fScaleDTBySin2T);  
-    hstat->Fill(kThtCW,fCutOnDThetaX);  
-    hstat->Fill(kPhiOvl,fPhiOverlapCut);
-    hstat->Fill(kZEtaOvl,fZetaOverlap);
-    hstat->Fill(kNoOvl,fRemoveOverlaps);
-    //
-    hstat->Fill(kPhiRot,fPhiRot);
-    hstat->Fill(kInjScl,fInjScale);
-    hstat->Fill(kEtaCut,fEtaCut);
-    hstat->Fill(kZVMin,fZVertexMin);
-    hstat->Fill(kZVMax,fZVertexMax);
-    hstat->Fill(kTrcMin,fMultCutMin);    
-    hstat->Fill(kTrcMax,fMultCutMax);    
-    //
-    hstat->Fill(kOneUnit,1.);    
-  }
-  //
-  
+  RegisterStat();
   //  AliAnalysisTaskSE::Terminate();
 }
 
@@ -631,9 +670,13 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
   TH1F* hstat;
   //
   // ------------ job parameters, statistics ------------------------------>>>
-  hstat = new TH1F("hStat","Run statistics",kNStatBins,0.5,kNStatBins-0.5);
+  hstat = new TH1F("hStat","Run statistics",kNStatBins,0.5,kNStatBins+0.5);
   //
+  for (int ib=1;ib<=kNStatBins;ib++) hstat->GetXaxis()->SetBinLabel(ib,"-"); // dummy label
+  hstat->GetXaxis()->SetBinLabel(kEvTot0, "Ev.Tot0");
   hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
+  hstat->GetXaxis()->SetBinLabel(kEvTotPlp, "Ev.Tot Plp");
+  //
   hstat->GetXaxis()->SetBinLabel(kEvProcData,"Ev.ProcData");
   hstat->GetXaxis()->SetBinLabel(kEvProcInj,"Ev.ProcInj");
   hstat->GetXaxis()->SetBinLabel(kEvProcRot,"Ev.ProcRot");
@@ -650,7 +693,8 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
   //
   hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}");
   hstat->GetXaxis()->SetBinLabel(kInjScl,"inj");
-  hstat->GetXaxis()->SetBinLabel(kEtaCut,"#eta cut");
+  hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
+  hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
   hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut");
   hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut");
   hstat->GetXaxis()->SetBinLabel(kTrcMin,"Mult_{min} cut");
@@ -734,7 +778,7 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
   AddHisto(histos,hnV0,kHV0);
   //
   //----------------------------------------------------------------------
-  int nEtaBinsS = int(2*fEtaCut/0.1);
+  int nEtaBinsS = int((fEtaMax-fEtaMin)/0.1);
   if (nEtaBinsS<1) nEtaBins = 1;
   //
   int nZVBinsS = int(fZVertexMax-fZVertexMin);
@@ -742,7 +786,7 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
 
   if (fUseMC) {
     // Z vertex vs Eta distribution for primaries
-    TH2F* hzvetap = new  TH2F("zvEtaPrimMC","Z vertex vs eta PrimMC",nEtaBinsS,-fEtaCut,fEtaCut,nZVBinsS,fZVertexMin,fZVertexMax);
+    TH2F* hzvetap = new  TH2F("zvEtaPrimMC","Z vertex vs eta PrimMC",nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
     hzvetap->GetXaxis()->SetTitle("#eta");
     hzvetap->GetYaxis()->SetTitle("Zvertex");
     AddHisto(histos,hzvetap,kHZVEtaPrimMC);
@@ -750,6 +794,7 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
   //
   if (GetDoMixing()) {
     //
+    /*
     // Difference in Z vertex for mixed events
     TH1F* hzdiff = new TH1F("MixSPDVertexDiff","SPD #Delta Z Vertex distribution ",100,-5,5);
     hzdiff->GetXaxis()->SetTitle("#Delta Z Vertex [cm]");
@@ -760,6 +805,7 @@ TObjArray* AliTrackletTaskUni::BookCustomHistos()
     TH1F* hntdiff = new TH1F("MixNTrackletsDiff"," SPD tracklets Diff ",2000,-3000,3000);
     hntdiff->GetXaxis()->SetTitle("# tracklet diff");
     AddHisto(histos,hntdiff,kHNTrMixDiff);
+    */
   }
   // 
   // --------------------------------------------------
@@ -793,9 +839,9 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
   //
   const int kNDPhiBins = 100;
   const int kNDThtBins = 100;
-  int nDistBins = int(fNStdDev)*2;
+  int nDistBins = int(fNStdDev)*4;
 
-  int nEtaBins = int(2*fEtaCut/0.1);
+  int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
   if (nEtaBins<1) nEtaBins = 1;
   //
   int nZVBins = int(fZVertexMax-fZVertexMin);
@@ -813,8 +859,8 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
     sprintf(buffn,"%s_DistZvEta",pref);
     sprintf(bufft,"(%s)Weighted Dist.(#Delta) vs Zv vs Eta",pref);
     int nbnEZD[3]    = {nEtaBins,nZVBins,nDistBins};
-    double xmnEZD[3] = {-fEtaCut,fZVertexMin,0};
-    double xmxEZD[3] = { fEtaCut,fZVertexMax,fNStdDev};
+    double xmnEZD[3] = { fEtaMin,fZVertexMin,0};
+    double xmxEZD[3] = { fEtaMax,fZVertexMax,fNStdDev};
     hsp = new THnSparseF(buffn,bufft,3,nbnEZD,xmnEZD,xmxEZD);
     hsp->GetAxis(0)->SetTitle("#eta");
     hsp->GetAxis(1)->SetTitle("Zv");
@@ -829,8 +875,8 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
     sprintf(buffn,"%s_DistZvDPhiS",pref);
     sprintf(bufft,"(%s) #Delta#varphi-#delta_{#varphi} vs Zv vs Eta",pref);
     int nbnEZP[3]    = {nEtaBins,nZVBins, int(dphir*2/0.005)};
-    double xmnEZP[3] = {-fEtaCut,fZVertexMin,-dphir};
-    double xmxEZP[3] = { fEtaCut,fZVertexMax, dphir};
+    double xmnEZP[3] = { fEtaMin,fZVertexMin,-dphir};
+    double xmxEZP[3] = { fEtaMax,fZVertexMax, dphir};
     hsp = new THnSparseF(buffn,bufft,3,nbnEZP,xmnEZP,xmxEZP);
     hsp->GetAxis(0)->SetTitle("#eta");
     hsp->GetAxis(1)->SetTitle("Zv");
@@ -841,7 +887,7 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
   if (selHistos & (0x1<<kHEtaZvCut) ) {
     sprintf(buffn,"%s_ZvEtaCutT",pref);
     sprintf(bufft,"(%s) Zv vs Eta with tracklet cut",pref);
-    h2 = new TH2F(buffn,bufft,nEtaBins,-fEtaCut,fEtaCut, nZVBins, fZVertexMin,fZVertexMax);
+    h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
     h2->GetXaxis()->SetTitle("#eta");
     h2->GetYaxis()->SetTitle("Zv");
     AddHisto(histos,h2,kHEtaZvCut);
@@ -869,7 +915,7 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
   if (selHistos & (0x1<<kHEtaDPhiS) ) {
     sprintf(buffn,"%s_EtaDPhiS",pref);
     sprintf(bufft,"(%s) #Delta#varphi-#delta_{#varphi} vs #eta",pref);
-    h2 = new TH2F(buffn,bufft,nEtaBins, -fEtaCut,fEtaCut,kNDPhiBins,-dphir,dphir);
+    h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,kNDPhiBins,-dphir,dphir);
     h2->GetXaxis()->SetTitle("#eta");
     h2->GetYaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]");
     AddHisto(histos,h2,kHEtaDPhiS);
@@ -878,7 +924,7 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
   if (selHistos & (0x1<<kHEtaDThetaX) ) {
     sprintf(buffn,"%s_EtaDThetaX",pref);
     sprintf(bufft,"(%s) #Delta#theta%s vs #eta",pref,fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
-    h2 = new TH2F(buffn,bufft,nEtaBins, -fEtaCut,fEtaCut,kNDThtBins,-dthtr,dthtr);
+    h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,kNDThtBins,-dthtr,dthtr);
     h2->GetXaxis()->SetTitle("#eta");
     sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
     h2->GetYaxis()->SetTitle(bufft);
@@ -888,7 +934,7 @@ TObjArray* AliTrackletTaskUni::BookHistosSet(const char* pref, UInt_t selHistos)
   if (selHistos & (0x1<<kHEtaDist) ) {
     sprintf(buffn,"%s_EtaDist",pref);
     sprintf(bufft,"(%s) Weighted Dist.(#Delta) vs #eta",pref);
-    h2 = new TH2F(buffn,bufft,nEtaBins, -fEtaCut,fEtaCut,nDistBins,0,fNStdDev);
+    h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,nDistBins,0,fNStdDev);
     h2->GetXaxis()->SetTitle("#eta");
     sprintf(bufft,"#Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
            "[#Delta#theta%s/#sigma#theta]^{2}",fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
@@ -1069,7 +1115,7 @@ void AliTrackletTaskUni::FillHistosSet(TObjArray* histos, double /*phi*/,double
   double dphiS  = TMath::Abs(dphi) - fDPhiShift;
   if (dphi<0) dphiS = -dphiS;
   double eta    = -TMath::Log(TMath::Tan(theta/2));
-  if (TMath::Abs(eta)>fEtaCut) return;
+  if (eta<fEtaMin || eta>fEtaMax) return;
   //
   double dThetaX = dtheta;
   if (fScaleDTBySin2T) {
@@ -1089,6 +1135,9 @@ void AliTrackletTaskUni::FillHistosSet(TObjArray* histos, double /*phi*/,double
     ((THnSparseF*)histos->UncheckedAt(kHEtaZvDPhiS))->Fill(ezp);
   }
   //
+  if (histos->UncheckedAt(kHDPhiDTheta)) 
+    ((TH2F*)histos->UncheckedAt(kHDPhiDTheta))->Fill(dphi,dtheta);
+  //
   if (histos->UncheckedAt(kHDPhiSDThetaX)) 
     ((TH2F*)histos->UncheckedAt(kHDPhiSDThetaX))->Fill(dphiS,dThetaX);
   //
@@ -1171,18 +1220,25 @@ Bool_t AliTrackletTaskUni::HaveCommonParent(const float* clLabs0,const float* cl
   int npars[2]={0,0};
   const float *labs[2] = {clLabs0,clLabs1};
   int ntr = fStack->GetNtrack();
+  //  printf("\nHaveCommonParent \n");
   for (int il=0;il<2;il++) {
+    
     for (int ilb=0;ilb<3;ilb++) {
       int lbl = (int)labs[il][ilb];
-      if (lbl<0 || lbl>=ntr) continue;
+      if (lbl<2 || lbl>=ntr) continue;
       //
       while (npars[il]<kMaxPar-1) {
-       pars[il][ npars[il]++ ] = lbl;
        TParticle* part = fStack->Particle(lbl);
        if (!part) break;
+       int code = TMath::Abs(part->GetPdgCode());
+       int q = (int)TMath::Abs(part->GetPDG()->Charge());
+       if (code==21 || code<10 || q==1 || q==2 || q==4 ) break;
+       //printf("%d/%d/%d: %4d (%d)%s|",il,ilb,npars[il],lbl,part->GetStatusCode(),part->GetName());
+       pars[il][ npars[il]++ ] = lbl;
        lbl = part->GetFirstMother();
        if (lbl<1 || lbl>=ntr) break;
       }
+      //      printf("\n");
     }
   }
   // compare array of parents
index 29cbe44..f30ca6c 100755 (executable)
@@ -67,7 +67,9 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
 
   // bins for saved parameters
   enum {kDummyBin,
-       kEvTot,       // events read
+       kEvTot0,      // events read
+       kEvTot,       // events read after vertex quality selection
+       kEvTotPlp,    // events with pile-up
        kEvProcData,  // events with data mult.object (ESD or reco)
        kEvProcInj,   // events Injected
        kEvProcRot,   // events Rotated
@@ -85,12 +87,15 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
        //
        kPhiRot,      // rotation phi
        kInjScl,      // injection scaling
-       kEtaCut,      // eta cut
+       kEtaMin,      // eta cut
+       kEtaMax,      // eta cut
        kZVMin,       // min ZVertex to process
        kZVMax,       // max ZVertex to process
        kTrcMin,      // min mult to process
        kTrcMax,      // max mult to process
        //
+       kMCV0Scale,   // scaling value for V0 in MC     
+       //
        kOneUnit=49,  // just 1 to track mergings
        kNWorkers=50, // n workers
        kNStatBins
@@ -103,6 +108,7 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
   virtual void  UserCreateOutputObjects();
   virtual void  UserExec(Option_t *option);
   virtual void  Terminate(Option_t *);
+  void          RegisterStat();
 
   void       SetUseMC(Bool_t mc = kFALSE)              {fUseMC = mc;}
   void       SetCheckReconstructables(Bool_t c=kFALSE) {fCheckReconstructables = c;}
@@ -122,8 +128,11 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
   void       SetPhiRot(float w=0)               {fPhiRot = w;}
   void       SetInjScale(Float_t s=1.)          {fInjScale = s>0? s:1.;}
   void       SetRemoveOverlaps(Bool_t w=kFALSE) {fRemoveOverlaps = w;}
+  void       SetScaleMCV0(Float_t s=1.0)        {fMCV0Scale = s;}  
   //
-  void       SetEtaCut(Float_t eta)             {fEtaCut = eta;}
+  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       SetMultCutMin(Int_t n=0)           {fMultCutMin = n;}
@@ -139,6 +148,7 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
   void       SetDoRotation(Bool_t v=kTRUE)      {fDoRotation = v;}
   void       SetDoMixing(Bool_t v=kTRUE)        {fDoMixing = v;}
   //
+  void       SetDontMerge(Bool_t v=kTRUE)       {fDontMerge = v;}
   /*
   void       SetTrigger(AliTriggerAnalysis::Trigger trigger)  { fTrigger = trigger; }
   void       SetMCCentralityBin(MCCentralityBin mccentrbin)   { fMCCentralityBin = mccentrbin;}
@@ -183,11 +193,13 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
   //
   // Settings for the reconstruction
   // tracklet reco settings
-  Float_t      fEtaCut;                    // histos filled only for this eta range
+  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
   Int_t        fMultCutMin;                // min mult in ESD to process?
   Int_t        fMultCutMax;                // max mult in ESD to process?
+  Float_t      fMCV0Scale;                 // scaling factor for V0 in MC
   //
   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
@@ -215,6 +227,7 @@ class AliTrackletTaskUni : public AliAnalysisTaskSE {
   Float_t      fCentrUpLim;                // to select centrality bin on data
   TString      fCentrEst;                  // to select centrality estimator
   */
+  Bool_t fDontMerge;                       // no merging requested
   static const char*  fgkPDGNames[];                //!pdg names
   static const Int_t  fgkPDGCodes[];                //!pdg codes
   //
diff --git a/PWG0/multVScentPbPb/MyAnalysisMacroUni.C b/PWG0/multVScentPbPb/MyAnalysisMacroUni.C
new file mode 100755 (executable)
index 0000000..428a0c5
--- /dev/null
@@ -0,0 +1,269 @@
+Bool_t needRecPoints = kFALSE;
+
+void MyAnalysisMacroUni
+(
+ TString dataset="/alice/sim/LHC10f8c_130844",
+ TString outFName = "trbg.root",
+ TString noMergeDir = "",
+ Int_t   nEvents   = -1,
+ Bool_t  useMC     = kTRUE,          // fill MC info (doRec=kTRUE)
+ Float_t etaMin     =-0.5,        // min eta range to fill in histos
+ Float_t etaMax     = 0.5,        // max eta range to fill in histos
+ Float_t zMin       = -7,         // process events with Z vertex min
+ Float_t zMax       =  7,         //                     max positions
+ Int_t  ntMin       = 1,
+ Int_t  ntMax       = 999999,
+ float  injScale    = 1.,     // inject injScale*Ncl(Lr1/Lr2) hits
+ Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
+ float  nStdDev     = 25.,         // 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,
+ Float_t scaleMCV0  = 0.8,     // rescale MC V0 to match data
+ Bool_t checkReconstructables = kFALSE//kTRUE, // fill histos for reconstructable (needs useMC and doRec) 
+ )
+{
+  //  
+  //
+  needRecPoints = kTRUE; //doRec || doInj || doRot || doMix;
+    //
+  printf("Start Analysis for %s, max %d Events, Event Cuts: %.1f<eta<%.1f, %.2f<Zv<%.2f\n",
+        dataset.Data(),nEvents,etaMin,etaMax,zMin,zMax);
+  printf("Tracklet cuts: dPhi:%.3f dTheta:%.3f phiShift:%.4f | Keep %.1f NstDev\n"
+        "Scale dTheta: %s\n", 
+        dphi,dtht,phishift,nStdDev,scaleDTheta ? "ON":"OFF");
+  //
+  printf("UseMC: %s. V0 scale: %.4f\n",useMC ? "ON":"OFF",scaleMCV0);
+  //
+  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("AliITSMultRecBg.cxx++");
+  gProof->Load("AliTrackletTaskUni.cxx++");
+  //
+  // load and run AddTask macro
+  gROOT->LoadMacro("AddMultTaskTrackletUni.C");
+  //
+  // create our task
+  AliTrackletTaskUni *mltTask = AddMultTaskTrackletUni(outFName.Data(),noMergeDir);
+  
+  //
+  mltTask->SetDoNormalReco(kTRUE);
+  mltTask->SetDoInjection(kTRUE);
+  mltTask->SetDoRotation(kFALSE);
+  mltTask->SetDoMixing(kFALSE);  
+  //
+  mltTask->SetUseMC(useMC);
+  mltTask->SetCheckReconstructables(checkReconstructables);
+  //
+  mltTask->SetEtaMin(etaMin);
+  mltTask->SetEtaMax(etaMax);
+  mltTask->SetZVertexMin(zMin);
+  mltTask->SetZVertexMax(zMax);
+  //
+  mltTask->SetMultCutMin(ntMin);
+  mltTask->SetMultCutMax(ntMax);
+  //
+  //  mltTask->SetNStdCut(cutSigNStd);
+  mltTask->SetScaleMCV0(scaleMCV0);
+  //
+  mltTask->SetScaleDThetaBySin2T(scaleDTheta);
+  mltTask->SetNStdDev(nStdDev);
+  mltTask->SetPhiWindow(dphi);
+  mltTask->SetThetaWindow(dtht);
+  mltTask->SetPhiShift(phishift);
+  mltTask->SetPhiOverlapCut(ovlPhiCut);
+  mltTask->SetZetaOverlapCut(ovlZetaCut);
+  mltTask->SetInjScale(injScale);
+  mltTask->SetRemoveOverlaps(remOvl);
+  //
+  printf("new Task: %p\n",mltTask);
+  //
+  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,0);
+  mltTask->SelectCollisionCandidates();//AliVEvent::kMB);
+  //
+  //  */
+  // Run analysis
+  mgr->InitAnalysis();
+  // process dataset  
+  mgr->StartAnalysis("proof", dataset.Data(), nEvents, 0); 
+  //
+  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("MyPool");
+  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/PWG0/multVScentPbPb/runAAFbgpp.C b/PWG0/multVScentPbPb/runAAFbgpp.C
new file mode 100755 (executable)
index 0000000..030a135
--- /dev/null
@@ -0,0 +1,96 @@
+//
+void runAAFbgpp(TString dataset="/alice/sim/LHC11d3_000146806", //"/alice/sim/LHC10f8c_130844",
+               TString outFName = "LHC11d3_000146806_v1.root",
+               Int_t   nEvents    = -1,//3000,
+               TString noMergeDir = "root://alicers01.cern.ch//tmp/myoutput/", // "" // if non-zero string, no merging is done
+               Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
+               float  nStdDev     = 25.,         // number of st.dev. for tracklet cut to keep
+               float  dphi        = 0.08,        // 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  injScale    = 1.,//0.7,    // inject injScale*Ncl(Lr1/Lr2) hits
+               //
+               Bool_t useMC  = kTRUE,           // fill MC info (doRec=kTRUE)
+               Float_t etaMin     = -2.,        // min eta range to fill in histos
+               Float_t etaMax     =  2.,        // max eta range to fill in histos
+               Float_t zMin       = -17,         // process events with Z vertex min
+               Float_t zMax       =  17,         //                     max positions
+               Float_t scaleMCV0  = 0.8,     // rescale MC V0 to match data
+               Float_t ntMin      =   1,         // process events with ESDmult 
+               Float_t ntMax      = 20000,       // within this range
+               //
+               Bool_t checkReconstructables = kFALSE,//kTRUE, // fill histos for reconstructable (needs useMC and doRec) 
+               // 
+               //---------------------------------------------------------------------------------
+               float  phiRot      = 3.14159e+00, // angle for bg. generation with rotation
+               float  ovlPhiCut   = 0.005, 
+               float  ovlZetaCut  = 0.05,
+               TString alirootVer = "VO_ALICE@AliRoot::v4-21-33-AN",
+               TString rootVer    = "default",//"VO_ALICE@ROOT::v5-27-06b",
+               //
+               TString proofCluster="shahoian@alice-caf.cern.ch"
+               ) 
+{ 
+  //  
+  Bool_t runLocal = kTRUE; // true only for local test mode
+  if (runLocal) {
+    dataset = "/default/shahoian/tstsim_LHC11d3_146806";
+    //dataset = "/default/shahoian/test";
+    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";
+  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("runAAFbgpp.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);
+  //
+  gROOT->LoadMacro("MyAnalysisMacroUni.C");
+
+  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);
+  }
+  MyAnalysisMacroUni(dataset,outFName,noMergeDir,nEvents,useMC,
+                    etaMin,etaMax,zMin,zMax,ntMin,ntMax,
+                    injScale,scaleDTheta,nStdDev,dphi,dtht,
+                    phishift,remOvl,ovlPhiCut,ovlZetaCut,scaleMCV0,checkReconstructables);
+  //
+}