// ------------------------------------------------------------------------- // ----- AsyUnpack source file ----- // ----- Base created 12/04/12 by P. Russotto ----- // ----- Written for the ALADiN ToF-Wall by A. Le Fevre, 11/2012 ----- // ------------------------------------------------------------------------- #include "ATOFUnpack.h" #include "FairLogger.h" #include "TClonesArray.h" #include #include "TSystem.h" #include "TMath.h" using std::cout; using std::endl; inline Double_t PhiWeight(float phi,int Zbound2) { // Weight to apply at a given phi (degrees) in order to recover a flat reaction plane Phi distribution // as function of Zbound2 - ALF, 12/2013 Double_t w=0.; if (iCurrentSystemIndex>=0 && Zbound2>0 && Zbound2<=PhiWeightVsZbound->GetXaxis()->GetXmax() && phi>=-180 && phi<180) w = PhiWeightVsZbound->GetBinContent(PhiWeightVsZbound->GetXaxis()->FindBin(Zbound2), PhiWeightVsZbound->GetYaxis()->FindBin(phi)); //w=1.; // test without correction return w; } ATOFUnpack::ATOFUnpack() :FairTask() { //init calibration } // ------------------------------------------------------------------------- ATOFUnpack::ATOFUnpack(const char* name, Int_t iVerbose) : FairTask(name, iVerbose), fATOF_evt(0), fRootATOFhitEvent(0), fRootATOFhitEventCopy(0), fRootATOFtrackEvent(0), fRootATOFtrackEventCopy(0), fRootATOFcalEvent(0), fRootATOFcalEventCopy(0), fRootATOFglobalEvent(0), fRootATOFglobalEventCopy(0), ffBufferClone(0), fMbsInfo(0), // ALF, 12/2013 fMbsInfoClone(0) // ALF, 12/2013 { } // ------------------------------------------------------------------------- ATOFUnpack::~ATOFUnpack() { } // ------------------------------------------------------------------------- InitStatus ATOFUnpack::Init() { fLogger->Info(MESSAGE_ORIGIN," ATOFUnpack::Init()------------------Start "); corr=false; string dir = gSystem->Getenv("VMCWORKDIR"); string atofcaldir = dir + "/calfiles/ATOF/"; TATOFroot tatofroot; // important - 11/2012 fATOF_evt = new TATOFEvent(); fCAMAC_evt = new TCAMACEvent(); // start detector string atofmapping = atofcaldir + "s394_May2011_twal_2.map"; //if (p_tmap==0) p_tmap = new TATOFparaDsc("p_tmap",new TATOFparMap()); // ATOF mapping fATOF_evt->LoadMapping( atofmapping.c_str() ); // this mapping corrects from an inversion which was already here at the beginning of May 2011 (14th); // otherwise, use s394_May2011_twal_1.map instead - ALF fATOF_evt->LoadCalibrationParameters(atofcaldir.c_str()); // 09/2012 // dans le futur, devra si necessaire etre appele avant chaque run pour adapter les parametres de mapping et d'etalonnage en fonction du numero de run fRootATOFhitEvent = new TRootATOFhitEvent(); // ATOF infos at the hit level fRootATOFtrackEvent = new TRootATOFtrackEvent(); // ATOF infos at the track level (front-rear wall correlations) fRootATOFcalEvent = new TRootATOFcalEvent(); // ATOF calibrated infos fRootATOFglobalEvent = new TRootATOFglobalEvent(); // ATOF calibrated global variables FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) Fatal("Init", "No FairRootManager"); ioman->Register("ATOFHITEVENT", "ATOF HIT EVENT TREE", fRootATOFhitEvent, kFALSE); // kFALSE means that it won't be written on disk ioman->Register("ATOFTRACKEVENT", "ATOF TRACK EVENT TREE", fRootATOFtrackEvent, kFALSE); // (possible to be done after in the merger) ioman->Register("ATOFCALEVENT", "ATOF CAL EVENT TREE", fRootATOFcalEvent, kFALSE); ioman->Register("ATOFGLOBAL", "ATOF_GLOBAL_TREE", fRootATOFglobalEvent, kFALSE); fRootATOFhitEventCopy = new TClonesArray("TRootATOFhitEvent"); fRootATOFtrackEventCopy = new TClonesArray("TRootATOFtrackEvent"); fRootATOFcalEventCopy = new TClonesArray("TRootATOFcalEvent"); fRootATOFglobalEventCopy = new TClonesArray("TRootATOFglobalEvent"); fRootATOFhitEventCopy=(TClonesArray*)fRootATOFhitEvent; fRootATOFtrackEventCopy=(TClonesArray*)fRootATOFtrackEvent; fRootATOFcalEventCopy=(TClonesArray*)fRootATOFcalEvent; fRootATOFglobalEventCopy=(TClonesArray*)fRootATOFglobalEvent; ioman->Register("ATOFHITEVENTCLONE", "ATOF HIT EVENT TREE CLONE", fRootATOFhitEventCopy, kFALSE); ioman->Register("ATOFTRACKEVENTCLONE", "ATOF TRACK EVENT TREE CLONE", fRootATOFtrackEventCopy, kFALSE); ioman->Register("ATOFCALEVENTCLONE", "ATOF CAL EVENT TREE CLONE", fRootATOFcalEventCopy, kFALSE); ioman->Register("ATOFGLOBALCLONE", "ATOF GLOBAL TREE CLONE", fRootATOFglobalEventCopy, kFALSE); ffBufferClone =(TClonesArray*) ioman->GetObject("EVTBUFFERCLONE"); ffBuffer = (TEvtBuffer*)ffBufferClone; if(!ffBufferClone)cout << "ATOF EVT BUFFER NOT FOUND"<GetObject("MBSINFO"); // ALF, 12/2013 fMbsInfo = (TMbsInfo*)fMbsInfoClone; // ALF, 12/2013 // ===== Load phi (vs Zbound) correction weights ===== TFile *finputPhiWeights = new TFile(Form("%s/%s",atofcaldir.c_str(),InputPhiWeight),"read"); // 12/2013 int run = fMbsInfo->GetRunNumber(); // ALF, 12/2013 //int EventNumber = fMbsInfo->GetCamacEventNumber(); iCurrentSystemIndex=-1; for (int is=0;is<3;is++) if (run>=RunMinSystem[is] && run<=RunMaxSystem[is]) iCurrentSystemIndex=is; if (iCurrentSystemIndex==-1) printf(": no current run number defined in the FairRunAna database.\n"); else{ //PhiWeightVsZbound = (TH2D*)finputPhiWeights->Get(Form("hWeight_%s",PhiWeightSystem[iCurrentSystemIndex])); // corrections obtained from raw ATOF phi vs Zbound maps with the condition that CHIMERA and (corrected) MBall are correlated in Phi_r.p., such as the phi distribution is homogeneous for a given Zbound - ALF, 12/2013 PhiWeightVsZbound = (TH2D*)finputPhiWeights->Get(Form("hZbetaWeight_%s",PhiWeightSystem[iCurrentSystemIndex])); // obtained from phi vs Zbound maps, weighted by Z*beta_lab (like the Q-vector) //PhiWeightVsZbound = (TH2D*)finputPhiWeights->Get(Form("hPhirpChiMBWeight_%s",PhiWeightSystem[iCurrentSystemIndex])); // corrections obtained from raw ATOF Phi (reaction plane) vs Zbound maps with the condition that CHIMERA and (corrected) MBall are correlated in Phi_r.p. - ALF, 12/2013 printf("Correction weights for phi (vs Zbound) are loaded for run number %d -> system = %s - ALF\n",run,PhiWeightSystem[iCurrentSystemIndex]); } fLogger->Info(MESSAGE_ORIGIN," ATOFUnpack::Init()------------------End "); cout << " ...return to continue..."<< endl; getchar(); } // ------------------------------------------------------------------------- void ATOFUnpack::SetParTask() { } // ------------------------------------------------------------------------- void ATOFUnpack::Finish() { // if(fIsLastEvent) // fLogger->Info(MESSAGE_ORIGIN," Finish all events in file"); } // ------------------------------------------------------------------------- void ATOFUnpack::Reset() { } // ------------------------------------------------------------------------- void ATOFUnpack::Exec(Option_t* opt) { //_____________________________ATOF subevent________________________________________________ //fATOF_evt->ProcessCurrentEventMBS(SubEventDataPtr,*sebuflenght); fATOF_evt->ProcessCurrentEventMBS(ffBuffer->ATOFBuffer,ffBuffer->ATOF_nrw); // ALF, 02/2013 (don't take long words here to stick to initial version) Int_t i_nhit = ((TATOFdatRaw*)p_datraw->Object())->NHit(); Int_t i_ntrk = 0; Int_t i_hit = 0; fRootATOFhitEvent->nadc = ((TATOFdatRaw*)p_datraw->Object())->NAdc(); fRootATOFhitEvent->ntdc = ((TATOFdatRaw*)p_datraw->Object())->NTdc(); fRootATOFhitEvent->ndrop = ((TATOFdatRaw*)p_datraw->Object())->NDrop(); //printf("nadc=%d ntdc=%d ndrop=%d nhit=%d\n",fRootATOFhitEvent->nadc,fRootATOFhitEvent->ntdc,fRootATOFhitEvent->ndrop,i_nhit); // test // Initialisations of global variables for (Int_t ii = 0; ii < i_nhit; ii++) { fRootATOFhitEvent->slat[i_hit] = 0; fRootATOFhitEvent->stat[i_hit] = 0; fRootATOFhitEvent->adct[i_hit] = 0; fRootATOFhitEvent->adcb[i_hit] = 0; fRootATOFhitEvent->tdct[i_hit] = 0; fRootATOFhitEvent->tdcb[i_hit] = 0; fRootATOFtrackEvent->slatc[i_ntrk] = 0; fRootATOFtrackEvent->statc[i_ntrk] = 0; fRootATOFtrackEvent->adctc[i_ntrk] = 0; fRootATOFtrackEvent->adcbc[i_ntrk] = 0; fRootATOFtrackEvent->tdctc[i_ntrk] = 0; fRootATOFtrackEvent->tdcbc[i_ntrk] = 0; fRootATOFtrackEvent->slatnl[i_ntrk] = 0; fRootATOFtrackEvent->statnl[i_ntrk] = 0; fRootATOFtrackEvent->adctnl[i_ntrk] = 0; fRootATOFtrackEvent->adcbnl[i_ntrk] = 0; fRootATOFtrackEvent->tdctnl[i_ntrk] = 0; fRootATOFtrackEvent->tdcbnl[i_ntrk] = 0; fRootATOFtrackEvent->slatnr[i_ntrk] = 0; fRootATOFtrackEvent->statnr[i_ntrk] = 0; fRootATOFtrackEvent->adctnr[i_ntrk] = 0; fRootATOFtrackEvent->adcbnr[i_ntrk] = 0; fRootATOFtrackEvent->tdctnr[i_ntrk] = 0; fRootATOFtrackEvent->tdcbnr[i_ntrk] = 0; fRootATOFtrackEvent->slatoc[i_ntrk] = 0; fRootATOFtrackEvent->statoc[i_ntrk] = 0; fRootATOFtrackEvent->adctoc[i_ntrk] = 0; fRootATOFtrackEvent->adcboc[i_ntrk] = 0; fRootATOFtrackEvent->tdctoc[i_ntrk] = 0; fRootATOFtrackEvent->tdcboc[i_ntrk] = 0; fRootATOFtrackEvent->slatol[i_ntrk] = 0; fRootATOFtrackEvent->statol[i_ntrk] = 0; fRootATOFtrackEvent->adctol[i_ntrk] = 0; fRootATOFtrackEvent->adcbol[i_ntrk] = 0; fRootATOFtrackEvent->tdctol[i_ntrk] = 0; fRootATOFtrackEvent->tdcbol[i_ntrk] = 0; fRootATOFtrackEvent->slator[i_ntrk] = 0; fRootATOFtrackEvent->stator[i_ntrk] = 0; fRootATOFtrackEvent->adctor[i_ntrk] = 0; fRootATOFtrackEvent->adcbor[i_ntrk] = 0; fRootATOFtrackEvent->tdctor[i_ntrk] = 0; fRootATOFtrackEvent->tdcbor[i_ntrk] = 0; const TATOFrawHit& hit = ((TATOFdatRaw*)p_datraw->Object())->Hit(ii); Int_t i_slat = hit.Slat(); if (i_slat > 0 && i_slat < ((TATOFparMap*)p_tmap->Object())->MaxSlat()) { // condition on slat number // Hit level ============================================= fRootATOFhitEvent->slat[i_hit] = hit.Slat(); fRootATOFhitEvent->stat[i_hit] = hit.Status(); fRootATOFhitEvent->adct[i_hit] = hit.AdctScaled(); fRootATOFhitEvent->adcb[i_hit] = hit.AdcbScaled(); fRootATOFhitEvent->tdct[i_hit] = hit.TdctScaled(); fRootATOFhitEvent->tdcb[i_hit] = hit.TdcbScaled(); // Track Level =========================================== //const TATOFparMapSlat& slatinfo = ((TATOFparMap*)p_tmap->Object())->SlatInfo(i_slat); Int_t i_type = ((TATOFparMap*)p_tmap->Object())->TypeFromSlat(i_slat); // Front (1) or Rear (2) Wall //Int_t i_type = slatinfo.Type(); // empty, replaced by function above - ALF, 01/2013 //printf("ICI 0.0: i_type=%d\n",i_type); if (i_type == 2) { // Rear Wall fRootATOFtrackEvent->slatc[i_ntrk] = i_slat; // rear central slat fRootATOFtrackEvent->statc[i_ntrk] = hit.Status(); fRootATOFtrackEvent->adctc[i_ntrk] = hit.AdctScaled(); fRootATOFtrackEvent->adcbc[i_ntrk] = hit.AdcbScaled(); fRootATOFtrackEvent->tdctc[i_ntrk] = hit.TdctScaled(); fRootATOFtrackEvent->tdcbc[i_ntrk] = hit.TdcbScaled(); Int_t i_snl = ((TATOFparMap*)p_datraw->Object())->NeighborThisLeftFromSlat(i_slat); //Int_t i_snl = slatinfo.NeighborThisLeft(); // empty, replaced by above function - ALF, 01/2013 Int_t i_snr = ((TATOFparMap*)p_datraw->Object())->NeighborThisRightFromSlat(i_slat); //Int_t i_snr = slatinfo.NeighborThisRight(); // empty, replaced by above function - ALF, 01/2013 Int_t i_inl = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_snl); // = -1 when no signal is Int_t i_inr = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_snr); // found in the neighbors //printf("ICI 0.1: i_slat=%d i_snl=%d i_snr=%d i_inl=%d i_inr=%d\n",i_slat,i_snl,i_snr,i_inl,i_inr); if (i_inl >= 0) { // rear next left slat const TATOFrawHit& hitnl = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_inl); fRootATOFtrackEvent->slatnl[i_ntrk] = i_snl; fRootATOFtrackEvent->statnl[i_ntrk] = hitnl.Status(); fRootATOFtrackEvent->adctnl[i_ntrk] = hitnl.AdctScaled(); fRootATOFtrackEvent->adcbnl[i_ntrk] = hitnl.AdcbScaled(); fRootATOFtrackEvent->tdctnl[i_ntrk] = hitnl.TdctScaled(); fRootATOFtrackEvent->tdcbnl[i_ntrk] = hitnl.TdcbScaled(); } else{ fRootATOFtrackEvent->slatnl[i_ntrk] = 0; fRootATOFtrackEvent->statnl[i_ntrk] = 0; fRootATOFtrackEvent->adctnl[i_ntrk] = 0; fRootATOFtrackEvent->adcbnl[i_ntrk] = 0; fRootATOFtrackEvent->tdctnl[i_ntrk] = 0; fRootATOFtrackEvent->tdcbnl[i_ntrk] = 0; } if (i_inr >= 0) { // rear next right slat const TATOFrawHit& hitnr = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_inr); fRootATOFtrackEvent->slatnr[i_ntrk] = i_snr; fRootATOFtrackEvent->statnr[i_ntrk] = hitnr.Status(); fRootATOFtrackEvent->adctnr[i_ntrk] = hitnr.AdctScaled(); fRootATOFtrackEvent->adcbnr[i_ntrk] = hitnr.AdcbScaled(); fRootATOFtrackEvent->tdctnr[i_ntrk] = hitnr.TdctScaled(); fRootATOFtrackEvent->tdcbnr[i_ntrk] = hitnr.TdcbScaled(); } else{ fRootATOFtrackEvent->slatnr[i_ntrk] = 0; fRootATOFtrackEvent->statnr[i_ntrk] = 0; fRootATOFtrackEvent->adctnr[i_ntrk] = 0; fRootATOFtrackEvent->adcbnr[i_ntrk] = 0; fRootATOFtrackEvent->tdctnr[i_ntrk] = 0; fRootATOFtrackEvent->tdcbnr[i_ntrk] = 0; } Int_t i_adcmax = 1; Int_t i_slatmax = 0; for (Int_t i_dslat = 0; i_dslat <=2; i_dslat++) { // look for the strongest hit in front +/- 2 slats Int_t i_sol = ((TATOFparMap*)p_datraw->Object())->NeighborOtherLeftFromSlat(i_slat) - i_dslat; //Int_t i_sol = slatinfo.NeighborOtherLeft() - i_dslat; // empty, replaced by above function - ALF, 01/2013 Int_t i_sor = ((TATOFparMap*)p_datraw->Object())->NeighborOtherRightFromSlat(i_slat) + i_dslat; //Int_t i_sor = slatinfo.NeighborOtherRight() + i_dslat; // empty, replaced by above function - ALF, 01/2013 Int_t i_iol = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_sol); Int_t i_ior = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_sor); if (i_iol >= 0) { const TATOFrawHit& hitol = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_iol); Int_t i_adcsum = int(sqrt(float(hitol.AdctScaled())*float(hitol.AdcbScaled()))); //Int_t i_adcsum = hitol.AdctScaled() + hitol.AdcbScaled(); if (i_adcsum > i_adcmax) { i_adcmax = i_adcsum; i_slatmax = i_sol; } } if (i_ior >= 0) { const TATOFrawHit& hitor = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_ior); Int_t i_adcsum = int(sqrt(float(hitor.AdctScaled())*float(hitor.AdcbScaled()))); //Int_t i_adcsum = hitor.AdctScaled() + hitor.AdcbScaled(); if (i_adcsum > i_adcmax) { i_adcmax = i_adcsum; i_slatmax = i_sor; } } } if (i_slatmax > 0) { // a hit in front region is correlated //const TATOFparMapSlat& slatmaxinfo = ((TATOFparMap*)p_tmap->Object())->SlatInfo(i_slatmax); Int_t i_sol = ((TATOFparMap*)p_tmap->Object())->NeighborThisLeftFromSlat(i_slatmax); //Int_t i_sol = slatmaxinfo.NeighborThisLeft(); Int_t i_sor = ((TATOFparMap*)p_tmap->Object())->NeighborThisRightFromSlat(i_slatmax); //Int_t i_sor = slatmaxinfo.NeighborThisRight(); Int_t i_ioc = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_slatmax); Int_t i_iol = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_sol); Int_t i_ior = ((TATOFdatRaw*)p_datraw->Object())->IndexFromSlat(i_sor); if (i_ioc >= 0) { // front (other) central slat const TATOFrawHit& hitoc = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_ioc); fRootATOFtrackEvent->slatoc[i_ntrk] = i_slatmax; fRootATOFtrackEvent->statoc[i_ntrk] = hitoc.Status(); fRootATOFtrackEvent->adctoc[i_ntrk] = hitoc.AdctScaled(); fRootATOFtrackEvent->adcboc[i_ntrk] = hitoc.AdcbScaled(); fRootATOFtrackEvent->tdctoc[i_ntrk] = hitoc.TdctScaled(); fRootATOFtrackEvent->tdcboc[i_ntrk] = hitoc.TdcbScaled(); } else{ fRootATOFtrackEvent->slatoc[i_ntrk] = 0; fRootATOFtrackEvent->statoc[i_ntrk] = 0; fRootATOFtrackEvent->adctoc[i_ntrk] = 0; fRootATOFtrackEvent->adcboc[i_ntrk] = 0; fRootATOFtrackEvent->tdctoc[i_ntrk] = 0; fRootATOFtrackEvent->tdcboc[i_ntrk] = 0; } if (i_iol >= 0) { // next left other slat const TATOFrawHit& hitol = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_iol); fRootATOFtrackEvent->slatol[i_ntrk] = i_sol; fRootATOFtrackEvent->statol[i_ntrk] = hitol.Status(); fRootATOFtrackEvent->adctol[i_ntrk] = hitol.AdctScaled(); fRootATOFtrackEvent->adcbol[i_ntrk] = hitol.AdcbScaled(); fRootATOFtrackEvent->tdctol[i_ntrk] = hitol.TdctScaled(); fRootATOFtrackEvent->tdcbol[i_ntrk] = hitol.TdcbScaled(); } else{ fRootATOFtrackEvent->slatol[i_ntrk] = 0; fRootATOFtrackEvent->statol[i_ntrk] = 0; fRootATOFtrackEvent->adctol[i_ntrk] = 0; fRootATOFtrackEvent->adcbol[i_ntrk] = 0; fRootATOFtrackEvent->tdctol[i_ntrk] = 0; fRootATOFtrackEvent->tdcbol[i_ntrk] = 0; } if (i_ior >= 0) { // next right other slat const TATOFrawHit& hitor = ((TATOFdatRaw*)p_datraw->Object())->Hit(i_ior); fRootATOFtrackEvent->slator[i_ntrk] = i_sor; fRootATOFtrackEvent->stator[i_ntrk] = hitor.Status(); fRootATOFtrackEvent->adctor[i_ntrk] = hitor.AdctScaled(); fRootATOFtrackEvent->adcbor[i_ntrk] = hitor.AdcbScaled(); fRootATOFtrackEvent->tdctor[i_ntrk] = hitor.TdctScaled(); fRootATOFtrackEvent->tdcbor[i_ntrk] = hitor.TdcbScaled(); } else{ fRootATOFtrackEvent->slator[i_ntrk] = 0; fRootATOFtrackEvent->stator[i_ntrk] = 0; fRootATOFtrackEvent->adctor[i_ntrk] = 0; fRootATOFtrackEvent->adcbor[i_ntrk] = 0; fRootATOFtrackEvent->tdctor[i_ntrk] = 0; fRootATOFtrackEvent->tdcbor[i_ntrk] = 0; } i_ntrk += 1; } } // } // end condition on slat number i_hit++; } // end loop over ToF-Wall hits fRootATOFhitEvent->nhit = i_hit; fRootATOFtrackEvent->ntrack = i_ntrk; //printf("ICI 2: nhit=%d ntrack=%d\n",fRootATOFhitEvent->nhit,fRootATOFtrackEvent->ntrack); ((TATOFdatRaw*)p_datraw->Object())->Clear(); // reset ToF-Wall event informations // does not work in AsyEOSRoot - ALF, 01/2012 //printf("HERE3\n"); // ======== Read the start time from the Stelzer detector ==================================== Float_t tstart=fCAMAC_evt->tstart_sc(ffBuffer->CAMACBuffer,ffBuffer->CAMAC_nrw); // ALF, 02/2013 (don't take long words here to stick to initial version) // =========================================================================================== // =====Tof-Wall event calibration. // ==== Calibrated informations at the track level. // ==== Strategy in 2 steps: // ==== 1- a 1st pass over tracks where only clean rear slats are accepted, recording the front corresponding slats in list 1. // ==== 2- a 2nd pass over rear tracks which are not clean, accepting clean front corresponding slats if not in list 1. // ==== (for the future, can be extracted also at the hit level in order to recover events // ==== where one wall has been hit through a crack => no signal there) //printf("HERE4\n"); fRootATOFcalEvent->npart=0; fRootATOFcalEvent->Zmax1=0; fRootATOFcalEvent->Zmax2=0; fRootATOFcalEvent->Zbound2=0; fRootATOFcalEvent->Ztot=0; fRootATOFcalEvent->nmultihit=0; fRootATOFcalEvent->multH=0; fRootATOFcalEvent->multHe=0; fRootATOFcalEvent->multIMF=0; Bool_t FrontExcluded[100]; // list of front slats which correspond to correct rear slat hits (which information is taken in priority) - index = slat number Bool_t RearOK[200],FrontOK[200]; // list of rear/front slats where the hit information is clean - index = track number int a1,a2,t1,t2,sl,st,imh,cp; float mh,tof,E,y,Z; for (sl=SlatMin;slnpart = fRootATOFtrackEvent->ntrack; //printf("ntracks=%d\n",fRootATOFcalEvent->npart); for (i_ntrk=0;i_ntrkntrack;i_ntrk++){ // 1st loop // Front center slat corresponding to center rear slat sl = fRootATOFtrackEvent->slatoc[i_ntrk]; st = fRootATOFtrackEvent->statoc[i_ntrk]; a1 = fRootATOFtrackEvent->adctoc[i_ntrk]; a2 = fRootATOFtrackEvent->adcboc[i_ntrk]; t1 = fRootATOFtrackEvent->tdctoc[i_ntrk]; t2 = fRootATOFtrackEvent->tdcboc[i_ntrk]; if (!fATOF_evt->Calibration(sl,st,a1,a2,t1,t2,tstart,Z,tof,E,y,imh,mh,cp)) cp=13; //if (cp==13&&sl<48) printf("FRONT: sl=%03d,st=%02d,a1=%04d,a2=%04d,t1=%04d,t2=%04d,tstart=%5.1f,Z=%6.1f,tof=%5.1f,E=%6.0f,imh=%01d,mh=%4.1f,cp=%04d\n",sl,st,a1,a2,t1,t2,tstart,Z,tof,E,imh,mh,cp); fRootATOFcalEvent->slatoc[i_ntrk] = sl; fRootATOFcalEvent->statoc[i_ntrk] = st; fRootATOFcalEvent->Zoc[i_ntrk] = Z; fRootATOFcalEvent->tofoc[i_ntrk] = tof; fRootATOFcalEvent->adcoc[i_ntrk] = E; fRootATOFcalEvent->yoc[i_ntrk] = y; fRootATOFcalEvent->imultihitoc[i_ntrk] = imh; fRootATOFcalEvent->multihitoc[i_ntrk] = mh; fRootATOFcalEvent->cpoc[i_ntrk] = cp; if (fRootATOFcalEvent->statoc[i_ntrk]==15 // 2 TDCs and 2 ADCs are present && fRootATOFcalEvent->cpoc[i_ntrk]==0 // Z id. OK && fRootATOFcalEvent->imultihitoc[i_ntrk]==0){ // no multihit FrontOK[i_ntrk]=kTRUE; } else{ FrontOK[i_ntrk]=kFALSE; } // Rear center slat sl = fRootATOFtrackEvent->slatc[i_ntrk]; st = fRootATOFtrackEvent->statc[i_ntrk]; a1 = fRootATOFtrackEvent->adctc[i_ntrk]; a2 = fRootATOFtrackEvent->adcbc[i_ntrk]; t1 = fRootATOFtrackEvent->tdctc[i_ntrk]; t2 = fRootATOFtrackEvent->tdcbc[i_ntrk]; if (!fATOF_evt->Calibration(sl,st,a1,a2,t1,t2,tstart,Z,tof,E,y,imh,mh,cp)) cp=13; //printf("REAR : sl=%03d,st=%02d,a1=%04d,a2=%04d,t1=%04d,t2=%04d,tstart=%5.1f,Z=%6.1f,tof=%5.1f,E=%6.0f,imh=%01d,mh=%4.1f,cp=%04d\n",sl,st,a1,a2,t1,t2,tstart,Z,tof,E,imh,mh,cp); fRootATOFcalEvent->slatc[i_ntrk] = sl; fRootATOFcalEvent->statc[i_ntrk] = st; fRootATOFcalEvent->Zc[i_ntrk] = Z; fRootATOFcalEvent->tofc[i_ntrk] = tof; fRootATOFcalEvent->adcc[i_ntrk] = E; fRootATOFcalEvent->yc[i_ntrk] = y; fRootATOFcalEvent->imultihitc[i_ntrk] = imh; fRootATOFcalEvent->multihitc[i_ntrk] = mh; fRootATOFcalEvent->cpc[i_ntrk] = cp; float Zfront = fRootATOFcalEvent->Zoc[i_ntrk]; float Zrear = fRootATOFcalEvent->Zc[i_ntrk]; if (fRootATOFcalEvent->statc[i_ntrk]==15 // 2 TDCs and 2 ADCs are present && fRootATOFcalEvent->cpc[i_ntrk]==0 // Z id. OK && fRootATOFcalEvent->imultihitc[i_ntrk]==0){ // no multihit // Checks of weird cases if (FrontOK[i_ntrk] && Zfront>2.5 && Zrear<1.5) // probably a track close to a rear crack -> take front information RearOK[i_ntrk]=kFALSE; else{ FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]=kTRUE; RearOK[i_ntrk]=kTRUE; } } else { RearOK[i_ntrk]=kFALSE; if (fRootATOFcalEvent->statc[i_ntrk]==15 && fRootATOFcalEvent->imultihitc[i_ntrk]==0 && fRootATOFcalEvent->cpc[i_ntrk]<2) { // rear nearly OK if (fRootATOFcalEvent->statoc[i_ntrk]==15 && fRootATOFcalEvent->imultihitoc[i_ntrk]==0 && fRootATOFcalEvent->cpoc[i_ntrk]<2) { // front nearly OK too -> check weird cases if (Zfront>2.5 && Zrear<1.5){ // probably a track close to a rear crack -> take front information RearOK[i_ntrk]=kFALSE; FrontOK[i_ntrk]=kTRUE; } else{ FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]=kTRUE; RearOK[i_ntrk]=kTRUE; } } else{ FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]=kTRUE; RearOK[i_ntrk]=kTRUE; } } } /* if (RearOK[i_ntrk]) printf("-> %d OK",fRootATOFcalEvent->slatc[i_ntrk]); if (FrontOK[i_ntrk]) printf(" %d OK ",fRootATOFcalEvent->slatoc[i_ntrk]); if (FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]) printf(" %d excluded ",fRootATOFcalEvent->slatoc[i_ntrk]); printf("\n"); */ } //printf("HERE5\n"); double Qx=0.,Qy=0.,Q1x=0.,Q1y=0.,Q2x=0.,Q2y=0.; double SumWeights=0.; float ZlongCM=0.,ZtransCM=0.,Zgbeta2longCM=0.,Zgbeta2transCM=0.; float Zlong=0.,Ztrans=0.,Zgbeta2long=0.,Zgbeta2trans=0.; Int_t NATOFaccepted=0; for (i_ntrk=0;i_ntrkntrack;i_ntrk++){ // 2nd loop over tracks // ====== Summary ====== float Zfront = fRootATOFcalEvent->Zoc[i_ntrk]; float Zrear = fRootATOFcalEvent->Zc[i_ntrk]; Bool_t TakeRear=kFALSE; Bool_t TakeFront=kFALSE; Bool_t Undefined=kFALSE; if (RearOK[i_ntrk]) TakeRear=kTRUE; // take rear information else if (FrontOK[i_ntrk] && !FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]) TakeFront=kTRUE; /* if (FrontOK[i_ntrk] // for the central 8-pack, take the FRONT wall && fRootATOFcalEvent->slatc[i_ntrk]>148 && fRootATOFcalEvent->slatc[i_ntrk]<=156){ // in priority, because TDCtop's of TakeRear=kFALSE; // REAR wall have quite often some weird values (up to 32000) TakeFront=kTRUE; // making that cpc=0 (when 0 < TDC < 4096) does not guaranty the absence of } // a "fake" TDC value; seem to be related to a bad TDC card slot - ALF, 01/2013 */ // NOT APPLIED because there are too many hits on the FRONT central 8-pack (from phi dist.) - ALF, 02/2013 if (TakeRear){ fRootATOFcalEvent->iZ[i_ntrk]=int(Zrear+0.5); fRootATOFcalEvent->Z[i_ntrk]=Zrear; fRootATOFcalEvent->tof[i_ntrk]=fRootATOFcalEvent->tofc[i_ntrk]; if (fRootATOFcalEvent->tofc[i_ntrk]!=0) fRootATOFcalEvent->beta[i_ntrk]=DistanceTargetRearWall/fRootATOFcalEvent->tofc[i_ntrk]/30.; else fRootATOFcalEvent->beta[i_ntrk]=-100.; fRootATOFcalEvent->adc[i_ntrk]=fRootATOFcalEvent->adcc[i_ntrk]; fRootATOFcalEvent->multihit[i_ntrk]=fRootATOFcalEvent->multihitc[i_ntrk]; fRootATOFcalEvent->cp[i_ntrk]=fRootATOFcalEvent->cpc[i_ntrk]; fRootATOFcalEvent->sl[i_ntrk]=fRootATOFcalEvent->slatc[i_ntrk]; fRootATOFcalEvent->x[i_ntrk]=-(fRootATOFcalEvent->sl[i_ntrk]-MiddleSlatRear)*2.5; // cm - in ALADiN convention (x>0 pointing to the left) - ALF, 03/2013 fRootATOFcalEvent->y[i_ntrk]=fRootATOFcalEvent->yc[i_ntrk];; // cm FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]=kTRUE; } else if (TakeFront){ fRootATOFcalEvent->iZ[i_ntrk]=int(Zfront+0.5); fRootATOFcalEvent->Z[i_ntrk]=Zfront; fRootATOFcalEvent->tof[i_ntrk]=fRootATOFcalEvent->tofoc[i_ntrk]; if (fRootATOFcalEvent->tofoc[i_ntrk]!=0) fRootATOFcalEvent->beta[i_ntrk]=DistanceTargetFrontWall/fRootATOFcalEvent->tofoc[i_ntrk]/30.; else fRootATOFcalEvent->beta[i_ntrk]=-100.; fRootATOFcalEvent->adc[i_ntrk]=fRootATOFcalEvent->adcoc[i_ntrk]; fRootATOFcalEvent->multihit[i_ntrk]=fRootATOFcalEvent->multihitoc[i_ntrk]; fRootATOFcalEvent->cp[i_ntrk]=fRootATOFcalEvent->cpoc[i_ntrk]; fRootATOFcalEvent->sl[i_ntrk]=fRootATOFcalEvent->slatoc[i_ntrk]; fRootATOFcalEvent->x[i_ntrk]=-(fRootATOFcalEvent->sl[i_ntrk]-MiddleSlatFront)*2.5; // cm - in ALADiN convention (x>0 pointing left) - ALF, 03/2013 fRootATOFcalEvent->y[i_ntrk]=fRootATOFcalEvent->yoc[i_ntrk];; // cm FrontExcluded[fRootATOFcalEvent->slatoc[i_ntrk]]=kTRUE; } else{ // undefined fRootATOFcalEvent->iZ[i_ntrk]=-100; fRootATOFcalEvent->Z[i_ntrk]==-100.; fRootATOFcalEvent->tof[i_ntrk]=-100.; fRootATOFcalEvent->beta[i_ntrk]=-100.; fRootATOFcalEvent->adc[i_ntrk]=-100.; fRootATOFcalEvent->multihit[i_ntrk]=-100.; fRootATOFcalEvent->cp[i_ntrk]=-100; fRootATOFcalEvent->sl[i_ntrk]=-100; fRootATOFcalEvent->x[i_ntrk]=-100.; fRootATOFcalEvent->y[i_ntrk]=-100.; fRootATOFcalEvent->theta[i_ntrk]=-100.; fRootATOFcalEvent->phi[i_ntrk]=-1000.; fRootATOFcalEvent->wphi[i_ntrk]=0.; Undefined=kTRUE; } if (!Undefined){ //if (TakeRear) printf(" ===> Z accepted=%d - rear - sl=%d\n",fRootATOFcalEvent->iZ[i_ntrk],fRootATOFcalEvent->slatc[i_ntrk]); //else if (TakeFront) printf(" ===> Z accepted=%d - front - sl=%d\n",fRootATOFcalEvent->iZ[i_ntrk],fRootATOFcalEvent->slatoc[i_ntrk]); Float_t r=sqrt(pow(fRootATOFcalEvent->x[i_ntrk],2)+pow(fRootATOFcalEvent->y[i_ntrk],2)); Float_t L; (fRootATOFcalEvent->sl[i_ntrk]>100) ? L=DistanceTargetRearWall : L=DistanceTargetFrontWall; fRootATOFcalEvent->theta[i_ntrk]=TMath::ATan(r/L)*TMath::RadToDeg(); if (r!=0.){ fRootATOFcalEvent->phi[i_ntrk]=TMath::ACos(fRootATOFcalEvent->x[i_ntrk]/r)*TMath::RadToDeg(); // for phi, adopt Paolo's convention: if (fRootATOFcalEvent->y[i_ntrk]<0) fRootATOFcalEvent->phi[i_ntrk]=-fRootATOFcalEvent->phi[i_ntrk]; // defined between -180 and 180 deg - ALF, 03/2013 //if (fRootATOFcalEvent->y[i_ntrk]<0) fRootATOFcalEvent->phi[i_ntrk]=360.-fRootATOFcalEvent->phi[i_ntrk]; // here for phi defined >= 0 } else fRootATOFcalEvent->phi[i_ntrk]=-1000.; if (fRootATOFcalEvent->iZ[i_ntrk]>=fRootATOFcalEvent->Zmax1){ fRootATOFcalEvent->Zmax2 = fRootATOFcalEvent->Zmax1; fRootATOFcalEvent->Zmax1 = fRootATOFcalEvent->iZ[i_ntrk]; } else if (fRootATOFcalEvent->iZ[i_ntrk]>fRootATOFcalEvent->Zmax2) fRootATOFcalEvent->Zmax2 = fRootATOFcalEvent->iZ[i_ntrk]; if (fRootATOFcalEvent->iZ[i_ntrk]>1) fRootATOFcalEvent->Zbound2 += fRootATOFcalEvent->iZ[i_ntrk]; fRootATOFcalEvent->Ztot += fRootATOFcalEvent->iZ[i_ntrk]; if (fRootATOFcalEvent->iZ[i_ntrk]==1) fRootATOFcalEvent->multH++; else if (fRootATOFcalEvent->iZ[i_ntrk]==2) fRootATOFcalEvent->multHe++; else if (fRootATOFcalEvent->iZ[i_ntrk]>2) fRootATOFcalEvent->multIMF++; float sinLAB = TMath::Sin(fRootATOFcalEvent->theta[i_ntrk]*TMath::DegToRad()); float cosLAB = TMath::Cos(fRootATOFcalEvent->theta[i_ntrk]*TMath::DegToRad()); float betaLAB=fRootATOFcalEvent->beta[i_ntrk]; float gammaLAB = (betaLAB<1) ? 1./sqrt(1.-pow(betaLAB,2)) : 0.; float betaTlab=betaLAB*sinLAB; float betaLlab=betaLAB*cosLAB; // Lorentz transformation into the centre-of-mass reference frame of the system. float LorentzFactor = 1.-betaSystCM*betaLlab; float gammaSystCM = 1./sqrt(1.-pow(betaSystCM,2)); float betaLcm = (LorentzFactor) ? (betaLlab-betaSystCM)/LorentzFactor : 0.; float betaTcm = (LorentzFactor) ? betaTlab/(gammaSystCM*LorentzFactor) : 0.; float betaCM=sqrt(pow(betaLcm,2)+pow(betaTcm,2)); float gammaCM = (betaCM<1) ? 1./sqrt(1.-pow(betaCM,2)) : 0.; float cosCM = (betaCM) ? betaLcm/betaCM : 0.; float sinCM = (betaCM) ? betaTcm/betaCM : 0.; float yCM = (betaLcm<1) ? 0.5*TMath::Log((1.+betaLcm)/(1.-betaLcm)) : 0.; // rapidity // Build reaction plane - ALF, 02/2013 ==================== int iZ=fRootATOFcalEvent->iZ[i_ntrk]; float phi=fRootATOFcalEvent->phi[i_ntrk]; float theta=fRootATOFcalEvent->theta[i_ntrk]; cp=fRootATOFcalEvent->cp[i_ntrk]; int imultihitc=fRootATOFcalEvent->imultihitc[i_ntrk]; int imultihitoc=fRootATOFcalEvent->imultihitoc[i_ntrk]; Double_t weight=double(iZ)*betaLAB; fRootATOFcalEvent->wphi[i_ntrk]=PhiWeight(phi,fRootATOFcalEvent->Zbound2); // ALF, 12/2013 SumWeights += fRootATOFcalEvent->wphi[i_ntrk]; // ALF, 12/2013 //if (TMath::Abs(slatc-152)<=2 || TMath::Abs(slatoc-52)<=2) weight=weight*CentralSlatWeight; double yrap=0.5*log((1+betaLlab)/(1-betaLlab)); if (phi==-100 || theta==-100) continue; // if (betaLAB>betaSystCM if (yrap>0.448 && betaLAB<1. && cp==0 && (imultihitc+imultihitoc)<=1 && theta>1 && theta<7 && iZ>1){ ZlongCM+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(cosCM,2); ZtransCM+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(sinCM,2); Zgbeta2longCM+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(gammaCM*betaLcm,2); Zgbeta2transCM+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(gammaCM*betaTcm,2); Zlong+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(cosLAB,2); Ztrans+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(sinLAB,2); Zgbeta2long+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(gammaLAB*betaLlab,2); Zgbeta2trans+=float(fRootATOFcalEvent->iZ[i_ntrk])*pow(gammaLAB*betaTlab,2); weight = weight * fRootATOFcalEvent->wphi[i_ntrk]; if (yCM>0.1){ // Danielewitcz-like cut-off Qy+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Sin(phi*TMath::DegToRad()); // phi in ALADiN convention: clockwise, Qx+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Cos(phi*TMath::DegToRad()); // 0 degree pointing to the left (looking downstream) if (NATOFaccepted%2){ // store separately odd-even index Q-vectors, for the determination of the reaction plane accuracy Q1y+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Sin(phi*TMath::DegToRad()); Q1x+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Cos(phi*TMath::DegToRad()); } else{ Q2y+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Sin(phi*TMath::DegToRad()); Q2x+=weight*TMath::Sin(theta*TMath::DegToRad())*TMath::Cos(phi*TMath::DegToRad()); } NATOFaccepted++; // - ALF, 03/2013 } } } else if (fRootATOFcalEvent->imultihitc[i_ntrk]==1 && fRootATOFcalEvent->imultihitoc[i_ntrk]==1) fRootATOFcalEvent->nmultihit++; // a multihit in both walls } // end 2nd loop over tracks // ALF, 12/2013 if (SumWeights>0){ Qx=Qx/SumWeights; Qy=Qy/SumWeights; Q1x=Q1x/SumWeights; Q1y=Q1y/SumWeights; Q2x=Q2x/SumWeights; Q2y=Q2y/SumWeights; } // ALF, 02/2013 double PhiATOF=-1000; double r = TMath::Sqrt(Qx*Qx+Qy*Qy); double r1 = TMath::Sqrt(Q1x*Q1x+Q1y*Q1y); double r2 = TMath::Sqrt(Q2x*Q2x+Q2y*Q2y); if (r>0 && NATOFaccepted>=MinNATOF && SumWeights>0){ PhiATOF=TMath::ACos(Qx/r)*TMath::RadToDeg(); // for Phi reaction plane, adopt Paolo's convention: defined between -180 and 180 deg - ALF, 03/2013 if (Qy<0.) PhiATOF=-PhiATOF; //if (Qy<0.) PhiATOF=360.-PhiATOF; // here for Phi defined >=0 } else { Qx=-1000.; Qy=-1000.; } if (r1==0 || NATOFacceptedphirp=PhiATOF; fRootATOFcalEvent->Nphirp=NATOFaccepted; fRootATOFcalEvent->Qx=Qx; fRootATOFcalEvent->Qy=Qy; fRootATOFcalEvent->Q1x=Q1x; fRootATOFcalEvent->Q1y=Q1y; fRootATOFcalEvent->Q2x=Q2x; fRootATOFcalEvent->Q2y=Q2y; fRootATOFcalEvent->ZLcm=ZlongCM; fRootATOFcalEvent->ZTcm=ZtransCM; fRootATOFcalEvent->Zgbeta2Lcm=Zgbeta2longCM; fRootATOFcalEvent->Zgbeta2Tcm=Zgbeta2transCM; fRootATOFcalEvent->ZL=Zlong; fRootATOFcalEvent->ZT=Ztrans; fRootATOFcalEvent->Zgbeta2L=Zgbeta2long; fRootATOFcalEvent->Zgbeta2T=Zgbeta2trans; fRootATOFglobalEvent->npart=fRootATOFcalEvent->npart; fRootATOFglobalEvent->Zmax1=fRootATOFcalEvent->Zmax1; fRootATOFglobalEvent->Zmax2=fRootATOFcalEvent->Zmax2; fRootATOFglobalEvent->Zbound2=fRootATOFcalEvent->Zbound2; fRootATOFglobalEvent->Ztot=fRootATOFcalEvent->Ztot; fRootATOFglobalEvent->nmultihit=fRootATOFcalEvent->nmultihit; fRootATOFglobalEvent->multH=fRootATOFcalEvent->multH; fRootATOFglobalEvent->multHe=fRootATOFcalEvent->multHe; fRootATOFglobalEvent->multIMF=fRootATOFcalEvent->multIMF; fRootATOFglobalEvent->phirp=fRootATOFcalEvent->phirp; fRootATOFglobalEvent->Nphirp=fRootATOFcalEvent->Nphirp; fRootATOFglobalEvent->Qx=fRootATOFcalEvent->Qx; fRootATOFglobalEvent->Qy=fRootATOFcalEvent->Qy; fRootATOFglobalEvent->Q1x=fRootATOFcalEvent->Q1x; fRootATOFglobalEvent->Q1y=fRootATOFcalEvent->Q1y; fRootATOFglobalEvent->Q2x=fRootATOFcalEvent->Q2x; fRootATOFglobalEvent->Q2y=fRootATOFcalEvent->Q2y; fRootATOFglobalEvent->ZLcm=fRootATOFcalEvent->ZLcm; fRootATOFglobalEvent->ZTcm=fRootATOFcalEvent->ZTcm; fRootATOFglobalEvent->Zgbeta2Lcm=fRootATOFcalEvent->Zgbeta2Lcm; fRootATOFglobalEvent->Zgbeta2Tcm=fRootATOFcalEvent->Zgbeta2Tcm; fRootATOFglobalEvent->ZL=fRootATOFcalEvent->ZL; fRootATOFglobalEvent->ZT=fRootATOFcalEvent->ZT; fRootATOFglobalEvent->Zgbeta2L=fRootATOFcalEvent->Zgbeta2L; fRootATOFglobalEvent->Zgbeta2T=fRootATOFcalEvent->Zgbeta2T; //printf("HERE6\n"); // ================================================================== //printf("============> npart=%d, Zmax1=%d, Zmax2=%d ============\n",fRootATOFcalEvent->npart,fRootATOFcalEvent->Zmax1,fRootATOFcalEvent->Zmax2) } // ------------------------------------------------------------------------- ClassImp(ATOFUnpack)