PWGCF/FLOW/macros/AddTaskPIDFlowSP.C

   1
   2 class AliAnalysisDataContainer;
   3 class AliFlowTrackCuts;
   4 class AliFlowTrackSimpleCuts;
   5 class AliFlowEventCuts;
   6 class AliFlowEventSimpleCuts;
   7
   8
   9 void AddTaskPIDFlowSP(Int_t triggerSelectionString=AliVEvent::kMB,
  10                                    Float_t etamin=-0.8,
  11                                    Float_t etamax=0.8,
  12                                    Float_t EtaGap=0.2,
  13                                    TString fileNameBase="AnalysisResults",
  14                                    TString uniqueStr="Pion_02",
  15                                    TString Qvector ="Qa",
  16                                    Int_t AODfilterBit = 272,
  17                                    Int_t charge=0,
  18                                    Int_t MinTPCdedx = 10,
  19                                    Int_t ncentrality = 6,
  20                                    Int_t maxITSCls = 7,
  21                                    Double_t maxChi2ITSCls = 37,
  22                                    Bool_t doQA=kTRUE,
  23                                    Bool_t isPID = kTRUE,
  24                                    Bool_t VZERO = kFALSE, // use vzero sp method
  25                                    Bool_t is2011 = kFALSE,
  26                                    Bool_t isAOD = kTRUE,
  27                                    Bool_t UsePIDParContours = kFALSE,
  28                                    AliPID::EParticleType particleType=AliPID::kPion,
  29                                    AliFlowTrackCuts::PIDsource sourcePID=AliFlowTrackCuts::kTOFbayesian) {
  30
  31 // Define a range of the detector to exclude
  32 Bool_t ExcludeRegion = kFALSE;
  33 Double_t excludeEtaMin = -0.;
  34 Double_t excludeEtaMax = 0.;
  35 Double_t excludePhiMin = 0.;
  36 Double_t excludePhiMax = 0.;
  37
  38 //Define the range for eta subevents (for SP method) with TPC
  39 Double_t minA = -0.8;//
  40 Double_t maxA = -0.5*EtaGap;//
  41 Double_t minB = +0.5*EtaGap;//
  42 Double_t maxB = +0.8;//
  43
  44 int centrMin[9] = {0,5,10,20,30,40,50,60,70};
  45 int centrMax[9] = {5,10,20,30,40,50,60,70,80};
  46 const int ncentr = ncentrality;
  47
  48
  49 //---------Data selection---------- ESD only!!!
  50 //kMC, kGlobal, kESD_TPConly, kESD_SPDtracklet
  51 AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kTPCstandalone;
  52 AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kTPCstandalone;
  53
  54 //---------Parameter mixing--------
  55 //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
  56 AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
  57 AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
  58
  59 const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype);  //ESD
  60 const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype); //ESD
  61
  62
  63 //===========================================================================
  64 // EVENTS CUTS:
  65 AliFlowEventCuts* cutsEvent[ncentr];
  66 AliFlowTrackCuts* cutsRP[ncentr];
  67 AliFlowTrackCuts* cutsPOI[ncentr];
  68 TString outputSlotName[ncentr][4];
  69 TString suffixName[ncentr];
  70
  71 for(int icentr=0;icentr<ncentr;icentr++){
  72     cutsEvent[icentr] = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
  73     cutsEvent[icentr]->SetLHC11h(is2011);
  74     cutsEvent[icentr]->SetCentralityPercentileRange(centrMin[icentr],centrMax[icentr]);
  75     cutsEvent[icentr]->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
  76     //  cutsEvent[icentr]->SetRefMultMethod(AliFlowEventCuts::kVZERO);
  77     //cutsEvent[icentr]->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
  78     //cutsEvent[icentr]->SetNContributorsRange(2);
  79     cutsEvent[icentr]->SetPrimaryVertexZrange(-10.,10.);
  80     cutsEvent[icentr]->SetQA(doQA);
  81     cutsEvent[icentr]->SetCutTPCmultiplicityOutliers();
  82
  83
  84     // RP TRACK CUTS:
  85     if(!VZERO){
  86         cutsRP[icentr] = new AliFlowTrackCuts(Form("RP_%d",icentr));
  87         //cutsRP[icentr]->SetParamType(rptype);
  88         //cutsRP[icentr]->SetParamMix(rpmix);
  89         cutsRP[icentr]->SetPtRange(0.2,5.);
  90         cutsRP[icentr]->SetEtaRange(etamin,etamax);
  91         cutsRP[icentr]->SetMinNClustersTPC(70);
  92         cutsRP[icentr]->SetMinChi2PerClusterTPC(0.1);
  93         cutsRP[icentr]->SetMaxChi2PerClusterTPC(4.0);
  94         cutsRP[icentr]->SetMaxDCAToVertexXY(2.4);
  95         cutsRP[icentr]->SetMaxDCAToVertexZ(3.0);
  96         cutsRP[icentr]->SetAcceptKinkDaughters(kFALSE);
  97         cutsRP[icentr]->SetMinimalTPCdedx(MinTPCdedx);
  98         cutsRP[icentr]->SetAODfilterBit(AODfilterBit);
  99     }
 100
 101     if(VZERO) { // use vzero sub analysis
 102         if(!is2011) cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010(); // select vzero tracks
 103         if(is2011)  cutsRP[icentr] = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011(); // select vzero tracks
 104
 105         if(!cutsRP[icentr]) {
 106             cout << " Fatal error: no RP cuts found! " << endl;
 107             return 0x0;
 108         }
 109
 110      }//vzero is not a tracking device it is just a scintillator. so pt range or DCAtoVertex are not set here.
 111     cutsRP[icentr]->SetQA(doQA);
 112
 113
 114     //POIs for SP and QC method
 115     //===========================================================================
 116     AliFlowTrackCuts  *SP_POI[ncentr];
 117     //half window for POIs
 118     //=======================SP POI Cuts
 119     SP_POI[icentr] = DefinePOIcuts(icentr);
 120
 121     SP_POI[icentr]->GetBayesianResponse()->ForceOldDedx(); // for 2010 data to use old TPC PID Response instead of the official one
 122     if(!isAOD){
 123         SP_POI[icentr]->SetMaxSharedITSCluster(maxITSCls);
 124         SP_POI[icentr]->SetMaxChi2perITSCluster(maxChi2ITSCls);
 125         SP_POI[icentr]->SetMinNClustersITS(2);
 126         SP_POI[icentr]->SetRequireITSRefit(kTRUE);
 127         SP_POI[icentr]->SetRequireTPCRefit(kTRUE);
 128         SP_POI[icentr]->SetMaxDCAToVertexXY(0.3);
 129         SP_POI[icentr]->SetMaxDCAToVertexZ(0.3);
 130         SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
 131         SP_POI[icentr]->SetMinimalTPCdedx(10.);
 132     }
 133     //SP_POI[icentr]->SetParamMix(poimix);
 134     SP_POI[icentr]->SetPtRange(0.2,6.);//
 135     SP_POI[icentr]->SetMinNClustersTPC(70);
 136     SP_POI[icentr]->SetMinChi2PerClusterTPC(0.1);
 137     SP_POI[icentr]->SetMaxChi2PerClusterTPC(4.0);
 138
 139
 140
 141     if(!VZERO && Qvector=="Qa"){
 142         SP_POI[icentr]->SetEtaRange( +0.5*EtaGap, etamax );
 143         printf(" > NOTE: Using half TPC (Qa) as POI selection u < \n");
 144
 145     }
 146     if(!VZERO && Qvector=="Qb"){
 147         SP_POI[icentr]->SetEtaRange( etamin,-0.5*EtaGap );
 148         printf(" > NOTE: Using half TPC (Qb) as POI selection u < \n");
 149
 150     }
 151     if(VZERO){
 152         SP_POI[icentr]->SetEtaRange( etamin,etamax );
 153         printf(" > NOTE: Using full TPC as POI selection u < \n");
 154     }
 155     //SP_POI->SetRequireITSRefit(kTRUE);
 156     //SP_POI->SetRequireTPCRefit(kTRUE);
 157     //SP_POI->SetMinNClustersITS(2);
 158     //SP_POI->SetMaxChi2PerClusterITS(1.e+09);
 159     SP_POI[icentr]->SetMaxDCAToVertexXY(2.4);
 160     SP_POI[icentr]->SetMaxDCAToVertexZ(3.0);
 161     //SP_POI->SetDCAToVertex2D(kTRUE);
 162     //SP_POI->SetMaxNsigmaToVertex(1.e+10);
 163     //SP_POI->SetRequireSigmaToVertex(kFALSE);
 164     SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
 165     if(isPID){
 166         SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
 167         SP_POI[icentr]->SetTPCTOFNsigmaPIDCutContours(UsePIDParContours,centrMin[icentr],centrMax[icentr]);
 168     }
 169
 170     if (charge!=0) SP_POI[icentr]->SetCharge(charge);
 171     //SP_POI->SetAllowTOFmismatch(kFALSE);
 172     SP_POI[icentr]->SetRequireStrictTOFTPCagreement(kTRUE);
 173     SP_POI[icentr]->SetMinimalTPCdedx(MinTPCdedx);
 174     if(isAOD) SP_POI[icentr]->SetAODfilterBit(AODfilterBit);
 175     SP_POI[icentr]->SetQA(doQA);
 176     SP_POI[icentr]->SetPriors((centrMin[icentr]+centrMax[icentr])*0.5);
 177
 178
 179
 180
 181     //=====================================================================
 182
 183     if(!VZERO && Qvector=="Qa") suffixName[icentr] = "Qa";
 184     if(!VZERO && Qvector=="Qb") suffixName[icentr] = "Qb";
 185     if(VZERO) suffixName[icentr] = "vzero";
 186     suffixName[icentr] += "-highharmflow";
 187     suffixName[icentr] += Form("%i_", centrMin[icentr]);
 188     suffixName[icentr] += Form("%i_", centrMax[icentr]);
 189     suffixName[icentr] += Form("%.f_", EtaGap*10);
 190
 191     if(isPID){
 192         suffixName[icentr]+=AliFlowTrackCuts::PIDsourceName(sourcePID);
 193         suffixName[icentr]+="_";
 194         suffixName[icentr]+=uniqueStr;
 195         suffixName[icentr]+="_";
 196         //suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
 197     }
 198     else{
 199         suffixName[icentr]+="AllCharged";
 200     }
 201     if (charge<0) suffixName[icentr]+="-";
 202     if (charge>0) suffixName[icentr]+="+";
 203
 204
 205     for(int harmonic=2;harmonic<6;harmonic++){  //for v2,v3,v4 and v5
 206         outputSlotName[icentr][harmonic-2] = "";
 207         outputSlotName[icentr][harmonic-2]+=uniqueStr;
 208         outputSlotName[icentr][harmonic-2]+=Form("_v%i_",harmonic);
 209         outputSlotName[icentr][harmonic-2]+=cutsRP[icentr]->GetName();
 210         outputSlotName[icentr][harmonic-2]+="_";
 211         outputSlotName[icentr][harmonic-2]+=SP_POI[icentr]->GetName();
 212         outputSlotName[icentr][harmonic-2]+=Form("_%i-",centrMin[icentr]);
 213         outputSlotName[icentr][harmonic-2]+=Form("%i_",centrMax[icentr]);
 214
 215
 216         if(isPID){
 217             outputSlotName[icentr][harmonic-2]+=AliFlowTrackCuts::PIDsourceName(sourcePID);//sourcePID
 218             outputSlotName[icentr][harmonic-2]+="_";
 219             outputSlotName[icentr][harmonic-2]+=AliPID::ParticleName(particleType);//particleType
 220         }
 221         else{
 222             outputSlotName[icentr][harmonic-2]+="AllCharged";
 223         }
 224         if (charge<0) outputSlotName[icentr][harmonic-2]+="-";
 225         if (charge>0) outputSlotName[icentr][harmonic-2]+="+";
 226     }
 227 }
 228
 229
 230 TString fileName(fileNameBase);
 231 fileName.Append(".root");
 232
 233
 234
 235 //====================================FLOWPACKAGE TASKS=========================//
 236 AliAnalysisDataContainer *cinput1[ncentr];
 237 AliAnalysisDataContainer *coutputFE[ncentr];
 238 AliAnalysisDataContainer* coutputFEQA[ncentr];
 239 AliAnalysisTaskFlowEvent *taskFE[ncentr];
 240
 241 AliAnalysisDataContainer *flowEvent[ncentr][4];
 242 AliAnalysisTaskFilterFE *tskFilter[ncentr][4];
 243
 244 AliAnalysisDataContainer *coutputSP[ncentr][4];
 245 AliAnalysisTaskScalarProduct *taskSP[ncentr][4];
 246
 247 TString outputQA[ncentr];
 248 TString myNameSP[ncentr][4];
 249 TString slot[ncentr][4];
 250
 251 for (int icentr=0; icentr<ncentr; icentr++) {
 252
 253     // Get the pointer to the existing analysis manager via the static access method.
 254         //==============================================================================
 255     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 256     if (!mgr) {
 257         Error("AddTaskFlowEvent", "No analysis manager to connect to.");
 258         return NULL;
 259     }
 260
 261     // Check the analysis type using the event handlers connected to the analysis
 262     // manager. The availability of MC handler can also be checked here.
 263     //==============================================================================
 264     if (!mgr->GetInputEventHandler()) {
 265         ::Error("AddTaskFlowEvent", "This task requires an input event handler");
 266         return NULL;
 267     }
 268
 269     taskFE[icentr] = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent_%s",suffixName[icentr].Data()),"",doQA);
 270     taskFE[icentr]->SelectCollisionCandidates(triggerSelectionString);
 271
 272   //  if(taskFE[icentr]->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl;
 273
 274     if(!VZERO) taskFE[icentr]->SetSubeventEtaRange(minA, maxA, minB, maxB);
 275     if(VZERO)  taskFE[icentr]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
 276     mgr->AddTask(taskFE[icentr]);
 277
 278     // Pass cuts for RPs and POIs to the task:
 279     taskFE[icentr]->SetCutsEvent(cutsEvent[icentr]);
 280     taskFE[icentr]->SetCutsRP(cutsRP[icentr]);
 281     taskFE[icentr]->SetCutsPOI(SP_POI[icentr]);
 282     if (cutsRP[icentr]->GetParamType()==AliFlowTrackCuts::kVZERO)
 283     {
 284         //TODO: since this is set in a static object all analyses in an analysis train
 285         //will be affected.
 286         taskFE[icentr]->SetHistWeightvsPhiMin(0.);
 287         taskFE[icentr]->SetHistWeightvsPhiMax(200.);
 288     }
 289
 290     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 291
 292     cinput1[icentr] = mgr->GetCommonInputContainer();
 293
 294     coutputFE[icentr] = mgr->CreateContainer(Form("FlowEvent_%s",suffixName[icentr].Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
 295
 296     mgr->ConnectInput(taskFE[icentr],0,cinput1[icentr]);
 297     mgr->ConnectOutput(taskFE[icentr],1,coutputFE[icentr]);
 298     //==========================================================
 299
 300     TString Species = "";
 301     if(isPID) Species += AliPID::ParticleName(particleType);
 302     else      Species += "Allcharged";
 303
 304
 305     for(int harm=2;harm<6;harm++){
 306         myNameSP[icentr][harm-2] = "SP_";
 307         myNameSP[icentr][harm-2] += Qvector;
 308         myNameSP[icentr][harm-2] += Form("_v%i_%s_%.f",harm,outputSlotName[icentr][harm-2].Data(),EtaGap*10);
 309
 310         flowEvent[icentr][harm-2] = mgr->CreateContainer( Form("Filter_%s", myNameSP[icentr][harm-2].Data()),
 311                                                                    AliFlowEventSimple::Class(),
 312                                                                    AliAnalysisManager::kExchangeContainer );
 313
 314         tskFilter[icentr][harm-2] = new AliAnalysisTaskFilterFE( Form("TaskFilter_%s",myNameSP[icentr][harm-2].Data()),cutsRP[icentr], NULL);
 315         if(!VZERO){
 316             tskFilter[icentr][harm-2]->SetSubeventEtaRange(etamin, -.5*EtaGap, +.5*EtaGap, etamax);
 317         }
 318         if(VZERO) tskFilter[icentr][harm-2]->SetSubeventEtaRange(-5,-1.5,+1.5,5);
 319         mgr->AddTask(tskFilter[icentr][harm-2]);
 320         mgr->ConnectInput( tskFilter[icentr][harm-2],0,coutputFE[icentr]);
 321         mgr->ConnectOutput(tskFilter[icentr][harm-2],1,flowEvent[icentr][harm-2]);
 322
 323
 324         taskSP[icentr][harm-2] = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s",outputSlotName[icentr][harm-2].Data()),kFALSE);
 325         taskSP[icentr][harm-2]->SetHarmonic(harm);
 326         taskSP[icentr][harm-2]->SelectCollisionCandidates(triggerSelectionString);
 327         taskSP[icentr][harm-2]->SetRelDiffMsub(1.0);
 328         taskSP[icentr][harm-2]->SetTotalQvector(Qvector);
 329         taskSP[icentr][harm-2]->SetApplyCorrectionForNUA(kTRUE);
 330
 331         TString outputSP = fileName;
 332         outputSP += ":outputSPanalysis";
 333         outputSP+= rptypestr;
 334         slot[icentr][harm-2] = "SP_";
 335         slot[icentr][harm-2] += outputSlotName[icentr][harm-2];
 336         slot[icentr][harm-2] += "_";
 337         slot[icentr][harm-2] += Qvector;
 338         coutputSP[icentr][harm-2] = mgr->CreateContainer(Form("%s_%.f",slot[icentr][harm-2].Data(),EtaGap*10),
 339                                          TList::Class(),AliAnalysisManager::kOutputContainer,outputSP);
 340         mgr->AddTask(taskSP[icentr][harm-2]);
 341         mgr->ConnectInput(taskSP[icentr][harm-2],0,flowEvent[icentr][harm-2]);
 342         mgr->ConnectInput(taskSP[icentr][harm-2],0,coutputFE[icentr]);
 343         mgr->ConnectOutput(taskSP[icentr][harm-2],1,coutputSP[icentr][harm-2]);
 344     }
 345
 346
 347     if (taskFE[icentr]->GetQAOn()) {
 348         outputQA[icentr] = fileName;
 349         outputQA[icentr] += ":QA";
 350         coutputFEQA[icentr] = mgr->CreateContainer(Form("QA_%s",suffixName[icentr].Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputQA[icentr]);
 351         mgr->ConnectOutput(taskFE[icentr],2,coutputFEQA[icentr]);
 352     }
 353
 354
 355 }
 356 }
 357
 358 //===========================================================================
 359
 360 AliFlowEventCuts* DefinecutsEvent(Int_t icentr){
 361     AliFlowEventCuts* cutsEvent = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
 362     return cutsEvent;
 363 }
 364 AliFlowTrackCuts* DefineRPcuts(Int_t icentr){
 365     AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts(Form("RP_%d",icentr));
 366     return cutsRP;
 367 }
 368 AliFlowTrackCuts* DefinePOIcuts(Int_t    icentr){
 369     AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts(Form("POI_%d",icentr));
 370     return cutsPOI;
 371 }
 372
 373
 374
 375