PWGHF/hfe/macros/AddTaskFlowTPCEMCalQCSP.C

   1 ///////////////////////////////////////////////////////////////////
   2 //                                                               //
   3 // AddTaskFlowTPCEMCalQCSP macro                                 //
   4 // Author: Andrea Dubla, Utrecht University, 2012                //
   5 //                                                               //
   6 ///////////////////////////////////////////////////////////////////
   7 class AliAnalysisDataContainer;
   8 class AliFlowTrackCuts;
   9 class AliFlowTrackSimpleCuts;
  10 class AliFlowEventCuts;
  11 class AliFlowEventSimpleCuts;
  12 class AliAnalysisDataContainer;
  13 class AliHFEextraCuts;
  14
  15 AliAnalysisTaskFlowTPCEMCalQCSP*  AddTaskFlowTPCEMCalQCSP(
  16                                                           TString uniqueID = "",
  17                                                           Float_t centrMin ,
  18                                                           Float_t centrMax ,
  19                                                           Double_t InvmassCut,
  20                                                           Double_t pTCut,
  21                                                           Int_t Trigger,
  22                                                           Bool_t multCorrcut,
  23                                                           Double_t minTPC,
  24                                                           Double_t maxTPC,
  25                                                           Double_t minEovP,
  26                                                           Double_t maxEovP,
  27                                                           Double_t minM20,
  28                                                           Double_t maxM20,
  29                                                           Double_t minM02,
  30                                                           Double_t maxM02,
  31                                                           Double_t Dispersion,
  32                                                           Int_t minTPCCluster,
  33                                                           AliHFEextraCuts::ITSPixel_t pixel,
  34                                                           Bool_t NUA = kTRUE,
  35                                                           Bool_t PhotonicElectronDCA = kFALSE,
  36                                                           Int_t TPCClusterforAsso = 80,
  37                                                           Bool_t AssoITSref = kTRUE,
  38                                                           Bool_t purity = kTRUE,
  39                                                           Bool_t SideBandsFlow = kFALSE,
  40                                                           Bool_t Phi_minus_psi = kFALSE,
  41                                                           const char *Cent = "V0M",
  42                                                           Bool_t QC = kTRUE, // use qc2 and qc4
  43                                                           Bool_t SP_TPC = kTRUE, //use tpc sp method
  44                                                           Bool_t VZERO_SP = kFALSE, // use vzero sp method
  45                                                           Bool_t BaseH = kFALSE, // base histo
  46                                                           Int_t harmonic = 2,
  47                                                           Bool_t shrinkSP = kTRUE,
  48                                                           Bool_t debug = kFALSE,
  49                                                           Int_t RPFilterBit = 1,
  50                                                           Bool_t op_ang = kFALSE,
  51                                                           Double_t op_angle_cut=3.,
  52                                                           TString histoflatname = "alien:///alice/cern.ch/user/a/adubla/CentrDistrBins005.root"
  53                                                           )
  54
  55 {
  56
  57
  58
  59     if(debug) cout << " === Adding Task ElectFlow === " << endl;
  60     TString fileName = AliAnalysisManager::GetCommonFileName();
  61     fileName += ":ElectroID_";
  62     fileName += uniqueID;
  63     if(debug) cout << "    --> Reconstruction data container: " << fileName << endl;
  64     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  65     if (!mgr) {
  66         if(debug) cout << " Fatal error: no analysis manager found! " << endl;
  67         return 0x0;
  68     }
  69     if (!mgr->GetInputEventHandler()) {
  70         if(debug) cout << " Fatal error: no imput event handler found!" << endl;
  71         return 0x0;
  72     }
  73
  74     //create a task
  75     AliAnalysisTaskFlowTPCEMCalQCSP *taskHFE = ConfigHFEemcalMod(kFALSE, minTPCCluster, pixel);    //kTRUE if MC
  76
  77     if(debug) cout << " === AliAnalysisElectronFlow === " << taskHFE << endl;
  78     if(!taskHFE) {
  79         if(debug) cout << " --> Unexpected error occurred: NO TASK WAS CREATED! (could be a library problem!) " << endl;
  80         return 0x0;
  81     }
  82     taskHFE->SetTrigger(Trigger);
  83
  84     if(Trigger==0 || Trigger==4){
  85         TFile *fFlat=TFile::Open(histoflatname.Data());
  86         TCanvas *c=fFlat->Get("cintegral");
  87         TH1F *hfl=(TH1F*)c->FindObject("hint");
  88         taskHFE->SetHistoForCentralityFlattening(hfl,centrMin,centrMax,0.,0);
  89     }
  90
  91     // Set centrality percentiles and method V0M, FMD, TRK, TKL, CL0, CL1, V0MvsFMD, TKLvsV0M, ZEMvsZDC
  92     taskHFE->SetCentralityParameters(centrMin, centrMax, Cent);
  93     taskHFE->SetInvariantMassCut(InvmassCut);
  94     taskHFE->SetIDCuts(minTPC, maxTPC, minEovP, maxEovP, minM20, maxM20, minM02, maxM02, Dispersion);
  95     taskHFE->SetFlowSideBands(SideBandsFlow);
  96     taskHFE->Setphiminuspsi(Phi_minus_psi);
  97     taskHFE->SetPurity(purity);
  98     taskHFE->SetpTCuttrack(pTCut);
  99     taskHFE->SelectPhotonicElectronMethod(PhotonicElectronDCA);
 100     taskHFE->SetOpeningAngleflag(op_ang);
 101     taskHFE->SetOpeningAngleCut(op_angle_cut);
 102     taskHFE->SetAssoTPCCluster(TPCClusterforAsso);
 103     taskHFE->SetAssoITSRefit(AssoITSref);
 104     taskHFE->SetMultCorrelationCut(multCorrcut);
 105     //set RP cuts for flow package analysis
 106     cutsRP = new AliFlowTrackCuts(Form("RFPcuts%s",uniqueID));
 107     if(!cutsRP) {
 108         if(debug) cout << " Fatal error: no RP cuts found, could be a library problem! " << endl;
 109         return 0x0;
 110     }
 111
 112     if(!VZERO_SP) {
 113         AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
 114         cutsRP->SetParamType(rptype);
 115         cutsRP->SetAODfilterBit(RPFilterBit);
 116         cutsRP->SetPtRange(0.2, 5.0);
 117         cutsRP->SetEtaRange(-0.7, 0.7);
 118         cutsRP->SetMinNClustersTPC(70);
 119         cutsRP->SetMinChi2PerClusterTPC(0.1);
 120         cutsRP->SetMaxChi2PerClusterTPC(4.0);
 121         cutsRP->SetRequireTPCRefit(kTRUE);
 122         cutsRP->SetMaxDCAToVertexXY(0.3);
 123         cutsRP->SetMaxDCAToVertexZ(0.3);
 124         cutsRP->SetAcceptKinkDaughters(kFALSE);
 125         cutsRP->SetMinimalTPCdedx(10.);
 126         if(debug) cout << "    --> kGlobal RP's " << cutsRP << endl;
 127     }
 128     if(VZERO_SP) { // use vzero sub analysis
 129         cutsRP = cutsRP->GetStandardVZEROOnlyTrackCuts(); // select vzero tracks
 130         SP_TPC = kFALSE; // disable other methods
 131         QC = kFALSE;
 132         if(debug) cout << "    --> VZERO RP's " << cutsRP << endl;
 133     }
 134
 135     AliFlowTrackSimpleCuts *POIfilterLeft = new AliFlowTrackSimpleCuts();
 136     AliFlowTrackSimpleCuts *POIfilterRight = new AliFlowTrackSimpleCuts();
 137     if(VZERO_SP || SP_TPC){
 138         POIfilterLeft->SetEtaMin(-0.7);
 139         POIfilterLeft->SetEtaMax(0.0);
 140         POIfilterLeft->SetMassMin(263731); POIfilterLeft->SetMassMax(263733);
 141
 142         POIfilterRight->SetEtaMin(0.0);
 143         POIfilterRight->SetEtaMax(0.7);
 144         POIfilterRight->SetMassMin(263731); POIfilterRight->SetMassMax(263733);
 145     }
 146
 147
 148     AliFlowTrackSimpleCuts *POIfilterQC = new AliFlowTrackSimpleCuts();
 149     if(QC){
 150         POIfilterQC->SetEtaMin(-0.7);
 151         POIfilterQC->SetEtaMax(0.7);
 152         POIfilterQC->SetMassMin(263731); POIfilterQC->SetMassMax(263733);
 153
 154     }
 155
 156     if(SideBandsFlow){
 157
 158         AliFlowTrackSimpleCuts *POIfilterLeftH = new AliFlowTrackSimpleCuts();
 159         AliFlowTrackSimpleCuts *POIfilterRightH = new AliFlowTrackSimpleCuts();
 160         if(SP_TPC){
 161             POIfilterLeftH->SetEtaMin(-0.7);
 162             POIfilterLeftH->SetEtaMax(0.0);
 163             POIfilterLeftH->SetMassMin(2636); POIfilterLeftH->SetMassMax(2638);
 164
 165             POIfilterRightH->SetEtaMin(0.0);
 166             POIfilterRightH->SetEtaMax(0.7);
 167             POIfilterRightH->SetMassMin(2636); POIfilterRightH->SetMassMax(2638);
 168         }
 169
 170
 171         AliFlowTrackSimpleCuts *POIfilterQCH = new AliFlowTrackSimpleCuts();
 172         if(QC){
 173             POIfilterQCH->SetEtaMin(-0.7);
 174             POIfilterQCH->SetEtaMax(0.7);
 175             POIfilterQCH->SetMassMin(2636); POIfilterQCH->SetMassMax(2638);
 176
 177         }
 178
 179     }
 180
 181     taskHFE->SetRPCuts(cutsRP);
 182
 183
 184     AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(Form("ccontainer0_%s",uniqueID.Data()),TList::Class(),AliAnalysisManager::kOutputContainer,fileName);
 185
 186     mgr->ConnectInput(taskHFE,0,mgr->GetCommonInputContainer());
 187     mgr->ConnectOutput(taskHFE,1,coutput3);
 188
 189
 190     if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl;
 191     AliAnalysisDataContainer *flowEvent = mgr->CreateContainer(Form("FlowContainer_%s",uniqueID.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
 192     mgr->ConnectOutput(taskHFE, 2, flowEvent);
 193     if(debug) cout << "    --> Created IO containers " << flowEvent << endl;
 194
 195     if(SideBandsFlow){
 196         if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl;
 197         AliAnalysisDataContainer *flowEventCont = mgr->CreateContainer(Form("FlowContainer_Cont_%s",uniqueID.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
 198         mgr->ConnectOutput(taskHFE, 3, flowEventCont);
 199         if(debug) cout << "    --> Created IO containers " << flowEventCont << endl;
 200     }
 201
 202
 203     mgr->AddTask(taskHFE);
 204
 205     if (QC) {  // add qc tasks
 206         AddQCmethod(Form("QCTPCin_%s",uniqueID.Data()), harmonic, flowEvent,  debug ,uniqueID, -0.7, -0.0, 0.0, 0.7,false,POIfilterQC, NUA, shrinkSP, BaseH);
 207         if(debug) cout << "    --> Hanging QC task ...succes! "<< endl;
 208     }
 209     if (SP_TPC) {  // add sp subevent tasks
 210         AddSPmethod(Form("SPTPCQa_in_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qa", harmonic, flowEvent, false, shrinkSP, debug,uniqueID, false, POIfilterRight, NUA,BaseH);
 211         if(debug) cout << "    --> Hanging SP Qa task ... succes!" << endl;
 212         AddSPmethod(Form("SPTPCQb_in_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qb", harmonic, flowEvent,  false, shrinkSP, debug,uniqueID, false, POIfilterLeft, NUA,BaseH);
 213         if(debug) cout << "    --> Hanging SP Qb task ... succes!"<< endl;
 214     }
 215     if (VZERO_SP) {  // add sp subevent tasks
 216         AddSPmethod(Form("SPVZEROQa_in_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qa", harmonic, flowEvent, false, shrinkSP, debug,uniqueID, true, POIfilterRight, NUA,BaseH);
 217         if(debug) cout << "    --> Hanging SP Qa task ... succes!" << endl;
 218         AddSPmethod(Form("SPVZEROQb_in_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qb", harmonic, flowEvent,  false, shrinkSP, debug,uniqueID, true, POIfilterLeft, NUA,BaseH);
 219         if(debug) cout << "    --> Hanging SP Qb task ... succes!"<< endl;
 220     }
 221
 222     //=========================================Flow event for elctronContamination==============================================================================================
 223     if(SideBandsFlow){
 224         if (QC) {  // add qc tasks
 225             AddQCmethod(Form("QCTPCCont_%s",uniqueID.Data()), harmonic, flowEventCont,  debug ,uniqueID, -0.7, -0.0, 0.0, 0.7,false,POIfilterQCH, NUA, shrinkSP,BaseH);
 226             if(debug) cout << "    --> Hanging QC task ...succes! "<< endl;
 227         }
 228         if (SP_TPC) {  // add sp subevent tasks
 229             AddSPmethod(Form("SPTPCQa_Cont_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qa", harmonic, flowEventCont, false, shrinkSP, debug,uniqueID, false, POIfilterRightH, NUA,BaseH);
 230             if(debug) cout << "    --> Hanging SP Qa task ... succes!" << endl;
 231             AddSPmethod(Form("SPTPCQb_Cont_%s", uniqueID.Data()), -0.7, -.0, .0, +0.7, "Qb", harmonic, flowEventCont,  false, shrinkSP, debug,uniqueID, false, POIfilterLeftH, NUA,BaseH);
 232             if(debug) cout << "    --> Hanging SP Qb task ... succes!"<< endl;
 233         }
 234     }
 235     //==========================================================================================================================================================================
 236
 237
 238     return taskHFE;
 239
 240 }
 241
 242 //_____________________________________________________________________________
 243
 244
 245 //_____________________________________________________________________________
 246 void AddSPmethod(char *name, double minEtaA, double maxEtaA, double minEtaB, double maxEtaB, char *Qvector, int harmonic, AliAnalysisDataContainer *flowEvent, bool bEP, bool shrink = false, bool debug, TString uniqueID, Bool_t VZERO_SP = kFALSE,  AliFlowTrackSimpleCuts* POIfilter, Bool_t kNUA, Bool_t kBaseH)
 247 {
 248     // add sp task and invm filter tasks
 249     if(debug) (bEP) ? cout << " ****** Reveived request for EP task ****** " << endl : cout << " ******* Switching to SP task ******* " << endl;
 250     TString fileName = AliAnalysisManager::GetCommonFileName();
 251     (bEP) ? fileName+=":EP" : fileName+=":SP";
 252     //          if(etagap) {
 253     //            fileName+="_SUBEVENTS";
 254     //          if(debug) cout << "    --> Setting up subevent analysis <-- " << endl;
 255     //    }
 256     if(debug) cout << "    --> fileName " << fileName << endl;
 257     TString myFolder = fileName;
 258     if(debug) cout << "    --> myFolder " << myFolder << endl;
 259     TString myNameSP;
 260     (bEP) ? myNameSP = Form("%sEPv%d%s", name, harmonic, Qvector): myNameSP = Form("%sSPv%d%s", name, harmonic, Qvector);
 261     if(debug) cout << " myNameSP " << myNameSP << endl;
 262     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 263     AliAnalysisDataContainer *flowEventOut = mgr->CreateContainer(Form("Filter_%s",myNameSP.Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
 264     AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myNameSP.Data()), NULL, POIfilter);
 265     tskFilter->SetSubeventEtaRange(minEtaA, maxEtaA, minEtaB, maxEtaB);
 266     if(VZERO_SP) tskFilter->SetSubeventEtaRange(-10, 0, 0, 10);
 267     mgr->AddTask(tskFilter);
 268     mgr->ConnectInput(tskFilter, 0, flowEvent);
 269     mgr->ConnectOutput(tskFilter, 1, flowEventOut);
 270     AliAnalysisDataContainer *outSP = mgr->CreateContainer(myNameSP.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
 271     AliAnalysisTaskScalarProduct *tskSP = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s", myNameSP.Data()), kFALSE);
 272     tskSP->SetApplyCorrectionForNUA(kNUA);
 273     tskSP->SetHarmonic(harmonic);
 274     tskSP->SetTotalQvector(Qvector);
 275     if (bEP) tskSP->SetBehaveAsEP();
 276     if (shrink)tskSP->SetBookOnlyBasicCCH(kBaseH);
 277     mgr->AddTask(tskSP);
 278     mgr->ConnectInput(tskSP, 0, flowEventOut);
 279     mgr->ConnectOutput(tskSP, 1, outSP);
 280 }
 281 //_____________________________________________________________________________
 282 void AddQCmethod(char *name, int harmonic, AliAnalysisDataContainer *flowEvent, Bool_t debug, TString uniqueID,double minEtaA, double maxEtaA, double minEtaB, double maxEtaB,Bool_t VZERO_SP = kFALSE,  AliFlowTrackSimpleCuts* POIfilter, Bool_t kNUA, bool shrink = false, Bool_t kBaseH)
 283 {
 284     // add qc task and invm filter tasks
 285     if(debug) cout << " ****** Received request for QC v" << harmonic << " task " << name << ", IO ****** " << flowEvent << endl;
 286     TString fileName = AliAnalysisManager::GetCommonFileName();
 287     fileName+=":QC";
 288     if(debug) cout << "    --> Common filename: " << fileName << endl;
 289     TString myFolder = Form("v%d", harmonic);
 290     if(debug) cout << "    --> myFolder: " << myFolder << endl;
 291     TString myName = Form("%s", name);
 292     if(debug) cout << "    --> myName: " << myName << endl;
 293     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 294     AliAnalysisDataContainer *flowEventOut = mgr->CreateContainer(Form("Filter_%s", myName.Data()), AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
 295     AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myName.Data()), NULL, POIfilter);
 296     tskFilter->SetSubeventEtaRange(minEtaA, maxEtaA, minEtaB, maxEtaB);
 297     //    if(VZERO_SP) tskFilter->SetSubeventEtaRange(-10, 0, 0, 10);
 298     mgr->AddTask(tskFilter);
 299     mgr->ConnectInput(tskFilter, 0, flowEvent);
 300     mgr->ConnectOutput(tskFilter, 1, flowEventOut);
 301
 302     AliAnalysisDataContainer *outQC = mgr->CreateContainer(myName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
 303     AliAnalysisTaskQCumulants *tskQC = new AliAnalysisTaskQCumulants(Form("TaskQCumulants_%s", myName.Data()), kFALSE);
 304     tskQC->SetApplyCorrectionForNUA(kNUA);
 305     tskQC->SetHarmonic(harmonic);
 306     if (shrink)tskQC->SetBookOnlyBasicCCH(kBaseH);
 307     mgr->AddTask(tskQC);
 308     mgr->ConnectInput(tskQC, 0, flowEventOut);
 309     mgr->ConnectOutput(tskQC, 1, outQC);
 310 }
 311 //_____________________________________________________________________________
 312
 313
 314 //_____________________________________________________________________________
 315
 316 AliAnalysisTaskFlowTPCEMCalQCSP* ConfigHFEemcalMod(Bool_t useMC,Int_t minTPCCulster,AliHFEextraCuts::ITSPixel_t pixel){
 317     //
 318     // HFE standard task configuration
 319     //
 320
 321     Bool_t kAnalyseTaggedTracks = kTRUE;
 322
 323     AliHFEcuts *hfecuts = new AliHFEcuts("hfeCutsEMCAL","HFE Standard Cuts");  //TODO....change the cuts values to PbPb
 324     //  hfecuts->CreateStandardCuts();
 325     hfecuts->SetMinNClustersTPC(minTPCCulster);
 326     hfecuts->SetMinNClustersITS(3);
 327     hfecuts->SetMinNTrackletsTRD(0);
 328     hfecuts->SetMinRatioTPCclusters(0.6);
 329
 330     //   hfecuts->SetEtaRange(-0.9,0.9);
 331     //   hfecuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 332     hfecuts->SetRequireITSPixel();
 333     hfecuts->SetCutITSpixel(pixel);//kAny
 334     hfecuts->SetMaxChi2perClusterITS(-1);
 335     hfecuts->SetMaxChi2perClusterTPC(3.5);
 336     hfecuts->SetCheckITSLayerStatus(kFALSE); // shud be put back
 337     //  hfecuts->UnsetVertexRequirement();
 338     hfecuts->SetVertexRange(10.);
 339     hfecuts->SetRequireSigmaToVertex();
 340     //hfecuts->SetSigmaToVertex(10);
 341     hfecuts->SetTOFPIDStep(kFALSE);
 342     //  hfecuts->SetQAOn();
 343     hfecuts->SetPtRange(0, 30);
 344
 345     AliAnalysisTaskFlowTPCEMCalQCSP *task = new AliAnalysisTaskFlowTPCEMCalQCSP("HFE_Flow_TPCEMCal");
 346     printf("task ------------------------ %p\n ", task);
 347
 348
 349     task->SetHFECuts(hfecuts);
 350
 351     //   task->SetInvariantMassCut(0.05);
 352     //  task->SetRejectKinkMother(kTRUE);
 353     //  task->SetRemovePileUp(kTRUE);
 354
 355     // Define PID
 356     AliHFEpid *pid = task->GetPID();
 357     if(useMC) pid->SetHasMCData(kTRUE);
 358     pid->AddDetector("TPC", 0);
 359     pid->AddDetector("EMCAL", 1);
 360
 361     printf("*************************************\n");
 362     printf("Configuring standard Task:\n");
 363     //  task->PrintStatus();
 364     pid->PrintStatus();
 365     printf("*************************************\n");
 366     return task;
 367
 368
 369 }
 370
 371 //_____________________________________________________________________________