PWG0/multVScentPbPb/MyAnalysisMacroTrackletMulti.C

   1 void MyAnalysisMacroTrackletMulti
   2 (TString dataset="/alice/sim/LHC10f8c_130844",
   3  TString outFName = "trbg.root",
   4  Int_t   nEvents   = -1,
   5  Float_t etaCut     = 0.5,        // max |eta| range to fill in histos
   6  Float_t zMin       = -7,         // process events with Z vertex min
   7  Float_t zMax       =  7,         //                     max positions
   8  Int_t   useCentVar = 0,          // centrality variable to use: 0=V0, 1=SPD2corr
   9  //
  10  Float_t cutSigNStd  = 1.5,        // cut on weighed distance used to extract signal
  11  Float_t cutSigDPhiS = -1,        // cut on dPhi-phiBent used to extract signal (if negative -> dphi*sqrt(cutSigNStd)
  12  Bool_t  useMC  = kTRUE,          // fill MC info (doRec=kTRUE)
  13  Float_t scaleMCV0 = 0.7520,      // rescale MC V0 to match the data
  14  //
  15  Bool_t doRec  = kTRUE,           // fill data histos from new reco
  16  Bool_t doInj  = kTRUE,           // create Inj. bg
  17  Bool_t doRot  = kFALSE,          // create Rot. bg
  18  Bool_t doMix  = kFALSE,//kTRUE,  // create Mix. bg
  19  //
  20  // specific parameters for reconstruction
  21  float  phiRot      = 3.14159e+00, // angle for bg. generation with rotation
  22  float  injScale    = 1.,//0.7,    // inject injScale*Ncl(Lr1/Lr2) hits
  23  Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
  24  float  nStdDev     = 25.,         // number of st.dev. for tracklet cut to keep
  25  float  dphi        = 0.06,        // dphi window (sigma of tracklet cut)
  26  float  dtht        = 0.025,       // dtheta .... (if negative, abs will be used with additional cut on |dthetaX|, apart from w.distance
  27  float  phishift    = 0.0045,      // bending shift
  28  Bool_t remOvl      = kTRUE,
  29  float  ovlPhiCut   = 0.005,
  30  float  ovlZetaCut  = 0.05,
  31  Int_t  nEventsSkip = 0,
  32  //----------------------- Ntracklets selection parameters important for mixing, to be tuned
  33  Float_t ntMin      =   1,         // process events with ESDmult
  34  Float_t ntMax      = 20000,       // within this range
  35  Float_t ntMixBinSz = 20000,       // ESDMult bin size for mixing
  36  //----------------------- Zv selection parameters important for mixing, to be tuned
  37  Float_t zMixBinSz  =  14,       //0.1,  // Zv. bin for mixing
  38  //---------------------------------------------------------------------------------
  39  //
  40  Bool_t checkReconstructables = kFALSE//kTRUE, // fill histos for reconstructable (needs useMC and doRec)
  41  //
  42  )
  43 {
  44   //
  45   if (cutSigDPhiS<0) cutSigDPhiS = TMath::Sqrt(cutSigNStd)*dphi;
  46   //
  47   printf("Start Analysis for %s, max %d Events skipping %d, Event Cuts: |eta|<%.1f, %.2f<Zv<%.2f\n",
  48          dataset.Data(),nEvents,nEventsSkip,etaCut,zMin,zMax);
  49   printf("Centrality variable: %d\n",useCentVar);
  50   printf("Tracklet cuts: dPhi:%.3f dTheta:%.3f phiShift:%.4f | Keep %.1f NstDev\n"
  51          "Scale dTheta: %s | Signal Selection: NstDev:%.1f, dPhiS: %.3f\n",
  52          dphi,dtht,phishift,nStdDev,scaleDTheta ? "ON":"OFF",
  53          cutSigNStd,cutSigDPhiS);
  54   //
  55   printf("UseMC: %s. V0 scale: %.4f\n",useMC ? "ON":"OFF",scaleMCV0);
  56   printf("Operations:  \n"
  57          "Reco:%s (RemOvl:%s phi:%.3f zeta:%.3f)\n"
  58          "Inj:%s (scale: %.2f)\n"
  59          "Rot:%s (phi: %.4f)\n"
  60          "Mix:%s (Nt:%.1f:%.1f:%.1f Zv:%.1f:%.1f:%.2f)\n",
  61          doRec ? "ON":"OFF",remOvl? "ON":"OFF",ovlPhiCut,ovlZetaCut,
  62          doInj ? "ON":"OFF",injScale,
  63          doRot ? "ON":"OFF",phiRot,
  64          doMix ? "ON":"OFF",ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
  65   //
  66   if (nEvents<0) nEvents = int(1e9);
  67   TString format = GetFormatFromDataSet(dataset);
  68   //
  69   // ALICE stuff
  70   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  71   if (!mgr) mgr = new AliAnalysisManager("Test train");
  72   //
  73   InputHandlerSetup(format,useMC);
  74   if (doMix) MixHandlerSetup(ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
  75   // compile our task
  76   gProof->Load("AliITSMultRecBg.cxx++");
  77   gProof->Load("AliTrackletTaskMulti.cxx++");
  78   //
  79   // load and run AddTask macro
  80   gROOT->LoadMacro("AddMultTaskTrackletMulti.C");
  81   //
  82   // create our task
  83   AliTrackletTaskMulti *mltTask = AddMultTaskTrackletMulti(outFName.Data());
  84   //
  85   mltTask->SetUseCentralityVar(useCentVar);
  86   mltTask->SetDoNormalReco(doRec);
  87   mltTask->SetDoInjection(doInj);
  88   mltTask->SetDoRotation(doRot);
  89   mltTask->SetDoMixing(doMix);
  90   //
  91   mltTask->SetUseMC(useMC);
  92   mltTask->SetCheckReconstructables(checkReconstructables);
  93   //
  94   mltTask->SetEtaCut(etaCut);
  95   mltTask->SetZVertexMin(zMin);
  96   mltTask->SetZVertexMax(zMax);
  97   //
  98   mltTask->SetDPhiSCut(cutSigDPhiS);
  99   mltTask->SetNStdCut(cutSigNStd);
 100   mltTask->SetScaleMCV0(scaleMCV0);
 101   //
 102   mltTask->SetScaleDThetaBySin2T(scaleDTheta);
 103   mltTask->SetNStdDev(nStdDev);
 104   mltTask->SetPhiWindow(dphi);
 105   mltTask->SetThetaWindow(dtht);
 106   mltTask->SetPhiShift(phishift);
 107   mltTask->SetPhiOverlapCut(ovlPhiCut);
 108   mltTask->SetZetaOverlapCut(ovlZetaCut);
 109   mltTask->SetPhiRot(phiRot);
 110   mltTask->SetInjScale(injScale);
 111   mltTask->SetRemoveOverlaps(remOvl);
 112   //
 113   printf("new Task: %p\n",mltTask);
 114   //
 115   AddPhysicsSelection(useMC);
 116   mltTask->SelectCollisionCandidates(useMC ? (AliVEvent::kMB) : (AliVEvent::kUserDefined) );
 117   //mltTask->SelectCollisionCandidates();//AliVEvent::kMB);
 118   //
 119   // Run analysis
 120   mgr->InitAnalysis();
 121   // process dataset
 122   mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip);
 123   //
 124   TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_";
 125   evstCmd += outFName;  evstCmd += " \nfi";
 126   gSystem->Exec( evstCmd.Data() );
 127
 128 }
 129
 130
 131 TString GetFormatFromDataSet(TString dataset) {
 132
 133 //   Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
 134   TString dsTreeName;
 135   if (dataset.Contains("#")) {
 136     Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
 137     dsTreeName=dataset(dataset.Last('#'),dataset.Length());
 138   } else {
 139     Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
 140     TFileCollection *ds = gProof->GetDataSet(dataset.Data());
 141     if (!ds) {
 142       Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
 143       return "";
 144     }
 145     dsTreeName = ds->GetDefaultTreeName();
 146   }
 147
 148   if (dsTreeName.Contains("esdTree")) {
 149     Info("runAAF.C","ESD input format detected ...");
 150     return "ESD";
 151   } else if (dsTreeName.Contains("aodTree"))  {
 152     Info("runAAF.C","AOD input format detected ...");
 153     return "AOD";
 154   } else {
 155     Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
 156     Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
 157   }
 158
 159   return "";
 160 }
 161
 162 Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
 163 {
 164   format.ToLower();
 165
 166   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 167
 168   AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
 169
 170   if (cin) return;
 171
 172   if (!format.CompareTo("esd"))
 173   {
 174     AliESDInputHandler *esdInputHandler = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 175
 176     if (!esdInputHandler)
 177     {
 178       Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
 179       esdInputHandler = new AliESDInputHandlerRP();
 180       mgr->SetInputEventHandler(esdInputHandler);
 181     }
 182
 183     if (useKine)
 184     {
 185       AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 186
 187       if (!mcInputHandler)
 188       {
 189         Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
 190         AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
 191         mgr->SetMCtruthEventHandler(mcInputHandler);
 192       }
 193     }
 194
 195   }
 196   else if (!format.CompareTo("aod"))
 197   {
 198     AliAODInputHandler *aodInputHandler = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 199
 200     if (!aodInputHandler)
 201     {
 202       Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
 203       aodInputHandler = new AliAODInputHandler();
 204       mgr->SetInputEventHandler(aodInputHandler);
 205     }
 206   }
 207   else
 208   {
 209     Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
 210     return kFALSE;
 211   }
 212
 213   return kTRUE;
 214 }
 215
 216 void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
 217                      float zMin, float zMax, float zMixBinSz)
 218 {
 219   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 220   if (!mgr) return;
 221   int bufferSize = 1;
 222   AliESDInputHandlerRP *esdH = dynamic_cast<AliESDInputHandlerRP*>(mgr->GetInputEventHandler());
 223   if (!esdH) return;
 224   //
 225   AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
 226   esdMixH->SetInputHandlerForMixing(esdH);
 227   AliMixEventPool *evPool = new AliMixEventPool("MyPool");
 228   AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
 229   AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
 230   //  evPool->AddCut(tracklets);
 231   evPool->AddCut(zvertex);
 232   //evPool->Init();
 233   evPool->Print();
 234   esdMixH->SetEventPool(evPool);
 235   esdH->SetMixingHandler(esdMixH);
 236 }
 237
 238 void AddPhysicsSelection(Bool_t isMC)
 239 {
 240   // physics selection a la Michele
 241   printf("Requesting physics selection in %s mode\n",isMC ? "MC":"Data");
 242   gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
 243   //isMC is true when processing monte carlo, the second 0 disables the cluster vs tracklets
 244   AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(isMC,0);
 245   if(!isMC) {
 246     AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
 247     physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
 248     /*
 249     physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
 250     physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
 251     */
 252     //
 253     physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
 254     physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
 255     //
 256     // This are needed only to fill the statistics tables
 257     physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
 258     physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
 259     physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
 260     //
 261     /*
 262     physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
 263     physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
 264     physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
 265     //
 266     physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
 267     physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
 268     physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
 269     //
 270     */
 271     /*
 272     //
 273     physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
 274     physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
 275     physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
 276     //
 277     physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
 278     physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
 279     physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
 280     */
 281   }
 282   // if you use the following line, your task only gets the selected events
 283   //  task->SelectCollisionCandidates(AliVEvent::kUserDefined);
 284   //
 285   //Alternatively, in the UserExec of your task:
 286   //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
 287   //
 288 }