PWGCF/FLOW/macros/AddTaskPIDFlowSP.C

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